---
title: "Simulation de processus de Poisson, de naissances et de morts"
author: S. Robin
date: septembre 2022
output: html_document
# output: pdf_document
# geometry: margin=1.5cm
# classoption: a4paper
---

# Processus de Poisson

On considère un processus de $\{N(t)\}_{t \geq 0}$ d'intensité $\lambda = 1/10$.
```{r poisPArm, echo=TRUE}
lambda <- 1/10
```

1 - Simuler une trajectoire de $N(t)$ comportant $n = 100$ sauts.

```{r poisNfixe, echo=TRUE}
n <- 1e2
duree <- rexp(n, lambda)
temps <- c(0, cumsum(duree))
population <- 0:n
plot(temps, population, type='s')
abline(0, lambda, col=2)
```

On peut en tracer plusieurs réalisations et comparer les trajectoires avec les quantiles (ici de niveau $.025$ et $.975$) de la loi de Poisson de paramètre $\lambda t$.

```{r poisNfixeRepet, echo=TRUE}
tGrid <- seq(0, round(n/lambda), length.out=100)
plot(tGrid, qpois(.975, lambda*tGrid), type='l', lty=2, col=2, xlab='temps', ylab='population')
lines(tGrid, lambda*tGrid, col=2)
lines(tGrid, qpois(.025, lambda*tGrid), col=2, lty=2)
for(i in 1:20){
  temps <- c(0, cumsum(rexp(n, lambda)))
  lines(temps, population, type='s', col=8)
}
```

2 - Simuler une trajectoire de $N(t)$ sur l'intervalle de temps $[0, 1000]$.

```{r poisTfixe, echo=TRUE}
tMax <- 1000
N <- rpois(1, lambda*tMax)
temps <- c(0, sort(tMax*runif(N)))
population <- 0:N
plot(temps, population, type='s')
abline(0, lambda, col=2)
```

# Processus de naissance

On considère un processus de naissances pure $\{N(t)\}_{t \geq 0}$ de taux de naissance linéaire : 
$$
\lambda(n) = \lambda \; n, \qquad \text{avec } \lambda = 1/10.
$$

```{r birthParm, echo=TRUE}
lambda <- 1/10
TauxNaissance <- function(x){lambda*x} # linéaire
# TauxNaissance <- function(x){lambda*x^2} # quadratique
```

1 - Simuler une trajectoire de $N(t)$ partant de un individu jusqu'à ce que la population atteigne la taille $n=100$.

```{r birthNfixe, echo=TRUE}
n0 <- 1; n <- 1e2
population <- n0:n
duree <- rexp(n-n0)
duree <- duree / TauxNaissance(population[-(n-n0)])
temps <- c(0, cumsum(duree))
plot(temps, population, type='s')
lines(temps, exp(lambda*temps), col=2)
```

```{r birthAlternative, echo=TRUE}
# ProcNaissance <- function(n, lambda){
#   duree <- rexp(n)
#   population <- 1:(n+1)
#   duree <- duree / TauxNaissance(population[1:n])
#   temps <- c(0, cumsum(duree))
#   return(list(temps=temps, population=population))
# }
# for(i in 1:10){
#   simul <- ProcNaissance(n, lambda)
#   points(c(0, simul$temps), c(simul$population, (n+1)), type='s', col=8)
# }
```

2 - Même question avec un taux de naissance quadratique
$$
\lambda(n) = \lambda \; n^2.
$$

```{r birthParmQuad, echo=TRUE}
TauxNaissance <- function(x){lambda*x^2} # quadratique
duree <- rexp(n-n0)
duree <- duree / TauxNaissance(population[-(n-n0)])
temps <- c(0, cumsum(duree))
plot(temps, population, type='s')
lines(temps, exp(lambda*temps), col=2)
```

# Processus de mort

On considère un processus de morts pure $\{N(t)\}_{t \geq 0}$ de taux de mort 
$$
\mu(n) = \mu \; n, \qquad \text{avec } \mu = 1/20.
$$

```{r deathParm, echo=TRUE}
mu <- 1/20
TauxMort <- function(x){mu*x}
```

1 - Simuler une trajectoire de $N(t)$ partant de $n = 100$ individus, jusqu'à extinction de la population.

```{r deathNinit, echo=TRUE}
populationInit <- 1e2
duree <- rexp(populationInit)
population <- populationInit:0
duree <- duree / TauxNaissance(population[1:populationInit])
temps <- c(0, cumsum(duree))
plot(temps, population, type='s')
lines(temps, populationInit*exp(-mu*temps), col=2)
```

# Processus de naissance et mort

On considère un processus de naissances et morts $\{N(t)\}_{t \geq 0}$ de taux de naissances $\lambda(n)$ et $\mu(n)$ linéaires
\begin{align*}
  \lambda(n) & = \lambda \; n & \text{avec} \quad \lambda & = 1/10, \\
  \text{et} \quad  \mu(n) & = \mu \; n & \text{avec} \quad \mu & = 2/25.
\end{align*}

```{r bdParm, echo=TRUE}
lambda <- 1/10; mu <- 2/25
TauxNaissance <- function(x){lambda*x}
TauxMort <- function(x){mu*x}
```

1 - Simuler une trajectoire du processus $N(t)$ partant d'une population initiale de taille $n(0) = 10$ et connaissant 1000 événements (naissances ou morts).

```{r bdSimul, echo=TRUE}
n <- 1e3; populationInit <- 1e1
population <- rep(0, (n+1)); population[1] = populationInit
temps <- rep(NA, (n+1)); temps[1] = 0
for(i in 1:n){
  if(population[i]>0){
    temps[i+1] <- temps[i] + rexp(1)/(TauxNaissance(population[i]) + TauxMort(population[i]))
    if(runif(1) < TauxNaissance(population[i]) / (TauxNaissance(population[i]) + TauxMort(population[i]))){
      population[i+1] <- population[i]+1 # naissance
      }else{
        population[i+1] <- population[i]-1 # mort
      }
  }
}  
plot(temps, population, type='s')
abline(h=0)
```

```{r bdAlternative, echo=TRUE}
# # Algo de Gillespie
# population <- rep(0, (n+1)); population[1] = populationInit
# for(i in 1:n){
#   if(population[i]>0){
#     if(runif(1) < TauxNaissance(population[i]) / (TauxNaissance(population[i]) + TauxMort(population[i]))){
#       population[i+1] <- population[i]+1 # naissance
#     }else{
#       population[i+1] <- population[i]-1 # mort
#     }
#   }
# }
# # Dates des événements
# duree <- rexp(n)
# duree <- duree / (TauxNaissance(population[1:n]) + TauxMort(population[1:n]))
# temps <- c(0, cumsum(duree))
```