---
title: "Simulation du processus de Lotka-Volterra stochastique"
author: S. Robin
date: septembre 2022
output: html_document
# output: pdf_document
# geometry: margin=1.5cm
# classoption: a4paper
---

Le modèle de Lotka-Volterra est un modèle déterministe décrivant l'évolution conjointe d'une population de proies (de taille $n_1(t)$) et d'une population de prédateurs (de taille $n_2(t)$) en supposant le système différentiel
$$
\left\{ \begin{array}{rcl} 
  n_1(t) & = & + a n_1(t) - b n_1(t) n_2(t), \\
  n_2(t) & = & - c n_2(t) + d n_1(t) n_2(t).
\end{array} \right.
$$

On définit un processus Markovien analogue décrivant l'évolution d'effectifs aléatoires $N_1(t)$ et $N_2(t)$. On suppose que ces effectifs sont chacun des processus de naissances et morts (non-indépendants).
On note $\lambda_1(n_1, n_2)$ le taux de naissance des proies et $\mu_1(n_1, n_2)$ leur taux de mort (idem $\lambda_2(n_1, n_2)$ et $\mu_2(n_1, n_2)$ pour les prédateurs) et, par analogie avec le modèle de Lotka-Volterra déterministe, on pose
\begin{align*}
  \lambda_1(n_1, n_2) & = a \; n_1, & 
  \mu_1(n_1, n_2) & = b \; n_1 \; n_2, \\
  \mu_2(n_1, n_2) & = c \; n_2, &
  \lambda_2(n_1, n_2) & = d \; n_1 \; n_2.
\end{align*}
avec
$$
a = 1, \qquad b = 1/100, \qquad c = 1, \qquad d = 1/100.
$$

```{r lvParm, echo=TRUE}
a <- 1; b <- 1/100; c <- 1; d <- 1/100
# a <- 1; b <- 1/100; c <- 1/2; d <- 5/1000
```

On appellera "état au temps $t$" le couple $(N_1(t), N_2(t)) \in \mathbb{N}^2$.

1 - Déterminer la liste des états vers les quels le processus peut, partant de l'état $(n_1, n_2)$,  transiter suite à un évenement et à quel taux se produit cette transition.

\begin{align*}
  (n_1, n_2) & \rightarrow (n_1+1, n_2) : \quad \lambda_1(n_1, n_2), \\ 
  (n_1, n_2) & \rightarrow (n_1-1, n_2) : \quad \mu_1(n_1, n_2), \\
  (n_1, n_2) & \rightarrow (n_1, n_2+1) : \quad \lambda_2(n_1, n_2), \\
  (n_1, n_2) & \rightarrow (n_1, n_2-1) : \quad \mu_2(n_1, n_2).
\end{align*}

```{r lvTaux, echo=TRUE}
mouvement <- matrix(c(1, 0, -1, 0, 0, 1, 0, -1), 4, 2, byrow=TRUE)
cat('Mouvements possibles :')
print(mouvement)
Taux <- function(etat){c(a*etat[1], b*etat[1]*etat[2], d*etat[1]*etat[2], c*etat[2])}
```

2 - Donner la loi du temps de séjour dans l'état $(n_1, n_2)$ et la probabilités de transiter vers chacun des état voisins.

$$
T_{n_1, n_2} \sim \mathcal{E}(\lambda_1 + \mu_1 + \lambda_2 + \mu_2)
$$
\begin{align*}
  (n_1, n_2) & \rightarrow (n_1+1, n_2) : \quad \lambda_1/(\lambda_1 + \mu_1 + \lambda_2 + \mu_2), \\ 
  (n_1, n_2) & \rightarrow (n_1-1, n_2) : \quad \mu_1/(\lambda_1 + \mu_1 + \lambda_2 + \mu_2), \\
  (n_1, n_2) & \rightarrow (n_1, n_2+1) : \quad \lambda_2/(\lambda_1 + \mu_1 + \lambda_2 + \mu_2), \\
  (n_1, n_2) & \rightarrow (n_1, n_2-1) : \quad \mu_2/(\lambda_1 + \mu_1 + \lambda_2 + \mu_2).
\end{align*}
```{r lvProbaTrans, echo=TRUE}
ProbTrans <- function(etat){Taux(etat)/sum(Taux(etat))}
```

3 - Simuler une suite d'états visités par le processus $\{(N_1(t), N_2(t))\}_{t \geq 0}$ partant de l'état $(n_1 = 200, n_2 = 50)$ et pour une succession de 1000 événements.

```{r lvSimulEtat, echo=TRUE}
n1init <- 200; n2init <- 50
nEvenements <- 50000
population <- matrix(NA, nEvenements, 2); population[1, ] <- c(n1init, n2init)
for(i in 2:nEvenements){
  population[i, ] <- population[i-1, ] + 
    mouvement[sample(1:4, 1, prob=ProbTrans(population[i-1, ])), ]
  population[i, which(population[i, ] < 0)] <- 0 
  if(sum(population[i, ]) == 0){break}
}
nEvenements <- i; population <- population[1:nEvenements, ]
head(population); tail(population)
plot(population, type='l', main=nEvenements, xlab='n1', ylab='n2')
points(n1init, n2init, col=2)
points(population[nEvenements, 1], population[nEvenements, 2], col=4)
```

3 - Simuler la suite des temps de séjours dans les états successivement visités.

```{r lvSimulTemps, echo=TRUE}
duree <- rexp(nEvenements-1)
for(i in 1:(nEvenements-1)){duree[i] <- duree[i] / sum(Taux(population[i, ]))}
temps <- c(0, cumsum(duree))
par(mfrow=c(1, 1), mex=.6)
plot(temps, population[, 1], col=2, type='l')
lines(temps, population[, 2], col=4)
```