# Processus stochastique : modèle SIR
###### Thomas Besognet, Clément Thion

Voici une représentation du modèle SIR qu'on modélise :

$$
\begin{array}{ccccc}
&& \boxed{\text{S}} & \overrightarrow{\beta_i /n} & \boxed{\text{I}} & \overrightarrow{\gamma} & \boxed{\text{R}} \\
\end{array}
$$

#### Question 1

On a $X_0=(100, 0)$. Par conséquent, il y a trois possibilités pour $X_1$ :

$$
\begin{split}
& X_1=(99, 1) \text{ : un individu tombe malade (A)} \\
& X_1=(100, -1) \text{ : un individu malade est retiré (B)} \\
& X_1=(100, 0) \text{ : la situation reste la même (C)} \\
\end{split}
$$

On va calculer la probabilité de passage pour chaque valeur de $X_1$. On note $P_i$ la probabilité de la possibilité $i$ de $X_1$. Par définition :

$$
\begin{split}
& P_A = \frac{\beta s_i}{N}=0 \text{ car } i=0 \\
& P_b = \gamma_i = 0 \text{ car } i=0 \\
& P_C = 1-(\frac{\beta s_i}{N} + \gamma_i) = 1 \\
\end{split}
$$

Donc forcément $X_1=X_0=(100, 0)$. Les possibilités pour $X_2$ sont alors les mêmes que pour $X_1$ avec les mêmes probabilités.
Par récurrence immédaite : $X_1=X_2 ... X_n=(100, 0)\ \forall n\in \mathbb{N}$. C'est biologiquement logique car si il n'y a pas de malade, la maladie ne peut pas se transmettre.

#### Question 2

En regardant chaque possibilité et en retirant les doublons, on trouve que :

$$
X_1 \in \left\{ \begin{array}{rcl}
(94, 6) \\
(95, 4) \\
(95, 5)
\end{array}\right.
X_2 \in \left\{ \begin{array}{rcl}
(93, 7) \\
(94, 5) \\
(94, 6) \\
(95, 3) \\
(95, 4) \\
(95, 5) \\
\end{array}\right.
X_3 \in \left\{ \begin{array}{rcl}
(92, 8) \\
(93, 6) \\
(93, 7) \\
(94, 4) \\
(94, 5) \\
(95, 2) \\
(95, 3) \\
(95, 4) \\
(95, 5) \\
(95, 6) \\
\end{array}\right.
$$

Pour $X_0=(95,5)$, on a l'impression que, si on appelle $V$ l'ensemble des valeurs possibles de $X_i$, on a $Card(V_i)=\sum_{i=0}^n (i+1)$ (*simple conjecture*).

#### Question 3

La chaîne n'est pas irréductible car il existe des états absorbants. Par exemple, l'une des possibilités pour $X_5$ est $X_5=(95,0)$. Or il s'agit d'un cas sans malade, comme à la question 1. A partir de ce cas, la situation est stabilisée et plus rien ne se passe : on ne peut plus atteindre les autres états.

#### Question 4
##### a)

On sait que dans le cas d'une loie uniforme sur $[0,1]$, $P(X\in [a, b])=b-a$. Donc si $U$ suit une loie uniforme sur $[0,1]$, on a trois cas possibles :

$$
\begin{split}
\text{A: }P(U\in[0, \frac{\beta s_i}{n}]) & = \frac{\beta s_i}{n} \\
& = P(X_{n+1}=(s-i, i+1)|X_n=(s,i)) \\
& = P(X_{suivant}=(s-1, i+1)) \\
\end{split}
$$

Correspond au cas où une personne susceptible devient inféctée par la maladie.
$$
\begin{split}
\text{B: }P(U\in[\frac{\beta s_i}{n}, \frac{\beta s_i}{n}+\gamma_i]) & = \frac{\beta s_i}{n}+\gamma_i-\frac{\beta s_i}{n} = \gamma_i \\
& = P(X_{n+1}=(s, i-1)|X_n=(s,i)) \\
& = P(X_{suivant}=(s, i-1)) \\
\end{split}
$$
Correspond au cas où une personne cesse d'être inféctée (passe dans les retirés, *i.e* meurt ou devient immunisée).

$$
\begin{split}
\text{C: }P(U\in ]\frac{\beta s_i}{n}+\gamma_i, 1]) & = 1-\frac{\beta s_i}{n}-\gamma_i\\
& = P(X_{n+1}=(s, i)|X_n=(s,i)) \\
& = P(X_{suivant}=(s, i)) \\
\end{split}
$$
Correspond au cas où il ne se passe rien (pas de nouvel infecté ni de nouveau retiré).

##### Questions 4.b et 4.c

Le code utilisé pour obtenir les graph suivant est disponible en annexe. On propose deux simulations pour $I_0=10$ et $I_0=20$.


On simule le processus pour $i=10$ et $i=20$. Les évolutions ainsi observées sont cohérentes ;
- $I$ le nombre d'infecté commence par augmenter (la maladie se répend parmis les susceptibles), puis diminue (*il y a davantage de personnes qui sont retirées que de personnes inféctées*), jusqu'à atteindre 0 (*plus d'évolution à partir de là).
- $S$ le nombre de susceptible diminue au cours du temps jusqu'à atteindre un pallier où il reste constant, correspondant au moment où il n'y a plus d'infecté.
- $R$ le nombre de retiré augmente à mesure que le nombre de susceptible diminue, jusqu'à atteindre lui aussi un palier constant au moment où il n'y a plus d'infection dans la population.

On sent bien que, plus $I_0$ est grand, plus le niveau des palliers constants vont être extrêmes (*Plus le nombre d'inféctée de départ est élevé, moins les susceptibles auront de chance d'échapper à une inffection, i.e plus $R$ sera haut et $S$ bas au moment où $I$ vaudra 0*).

#### Question 5

##### 5.a
Pour la rédaction des pseudo codes on réutilise la syntaxe employée dans l'ennoncé du sujet.

Fonction temps(i, MAX)
    """Donne le premier instant n où I=0
    """
    I[0]=i
    S[0]=N-i
    X[0]=(S[0],I[0])
    n=0
    tant que I[n] > 0 et que n<MAX: 
        U=random()
        augmenter la taille de I S X d'un item 
        Si U<= beta* S[n]*I[n]/N
            S[n+1]=S[n]-1
            I[n+1]=I[n]+1
        Sinon 
            Si U<= beta* S[n]*I[n]/N+gamma*I[n]
                S[n+1]=S[n]
                I[n+1]=I[n]-1
            Sinon
                S[n+1]=S[n]
                I[n+1]=I[n]
            Fin Si
        Fin Si
        X[n+1]=(S[n+1], I[n+1])
    Fin tant que
    retourner(n)
Fin fonction

##### 5.b

La loi forte des grands nombres énonce que $\displaystyle \lim_{n \to \infty} \frac{1}{n} \sum_{i=1}^{n}X_i \xrightarrow{PS} E(\bar{X_n})$, avec $\bar{X_n}$ la moyenne empirique de $(X_n)_{n\in \mathbb{N}}$, et $(X_n)_{n\in \mathbb{N}}$ une suite de variables aléatoires indépendantes et identiquement distribuées.

On a $T_i^{(l)}=\inf(n\geq 0, I_n^{(l)}=0)$. Par construction de $T_i^{(l)}$, $(T_i^{(l)})_{n\in [o,h],\ l\in [1,M]}$ est une suite de variables :
- *indépendantes* : à aucun moment $T_i^{(l)}$ ne dépend d'une valeur de $T_i^{(l-j)} \forall j$ ;
$T_i^{(l)}$ ne dépend que des paramètres du modèle $I_0, N, \beta, \gamma$
- *identiquement distribuée* : $I_0, N, \beta, \gamma$ sont maintenus constants pour les $M$ simulations.

(*On peut aussi remarquer qu'on nous dit dans l'énoncé que $(I_n)_{n\in[0,h]}$ est une suite de variable aléatoire i.i.d, donc $(T_i^{(l)})_{n\in [o,h],\ l\in [1,M]}$ est aussi une suite de v.a i.i.d puisque construit directement à partir de $(I_n)$*)

Par conséquent, on peut appliquer la loi forte des grands nombres à la suite de variables aléatoires $(T_n)$. On obtient que:

$\displaystyle \lim_{M \to +\infty} \frac{1}{M} \sum_{i=1}^{M}T_i^{(l)} \xrightarrow{PS} E(\bar{T_M}^{(l)})=E_{(N-i, i)}(T)$

##### 5.c

Une application pratique de la loi forte des grands nombres est la méthode de Monté Carlo : pour estimer $E_{(N-i,i)}(T)$, on prend comme estimateur du temps d'atteinte la moyenne empirique des $M$ simulations. C'est à dire que $\frac{1}{M} \sum_{i=1}^{M}\bar{T_i}^{(l)} \simeq E_{(N-i, i)}(T)$ pour $M$ grand (*Conséquence du point 5.b*).

Ce pseudo code est assez différent dans la forme à sa version encodée. En effet on a préféré séparer la partie Montecarlo de la partie simulation pour différentes valeurs de i. Pour autant voilà ce qui se passe en terme logique une fois ramené à une unique fonction. On n'a pas détaillé la partie calcule intervalle de confiance ni affichage puisque ce n'est pas la question.

Fonction temps_esp_montecarlo(i_min, i_max, M)
    """évolution, pour différentes valeurs de i,
        de l'espérance du temps d'atteinte E(T)
    """
    T_esti=[0,...,0] (autant que le nombre de valeurs de i testées)
    I_range=[i_min, i_min+1, ..., i_max]
    pour i0 allant de i_min à i_max
        pour M fois
            T[i0] = T[i0]+temps(i0, h)
        T[i0] = T[i0]/M
        fin pour
    fin pour
    
    calcul intervalle de confiance
    affichage
    

Voilà le résultat une fois encodé :


On a réalisé la même simulation mais pour $M=10000$ (*l'intervalle de confiance n'est plus visible*) :


##### 5.d

On trouve $i=5$ pour que le temps d'atteinte de $I_n=0$ soit le plus élevé.
Ce résultat nous parait biologiquement cohérent : si on avait un grand nombre d'infectés au début, l'épidémie irait forcément très vite et le virus se propagerait. A l'inverse , si on a trop peu d'infectés au départ, il y a un trop grand risque que ces deux personnes guérissent avant d'avoir contaminé quelqu'un et donc que l'épidémie ne dure vraiment pas longtemps. Il faut donc le minimum de contaminés nécessaire au départ pour que la maladie se répande.

##### 5.e

On note $(T_i)$ la suite des $E_{(N-i, i)}(T)$, $i\in[0,N]$ obtenu, pour chaque $i$, par méthode de Montecarlo.

$(T_i)_{n\in [o,h],\ i\in[0,N]}$ est une suite de variables aléatoires dépendant des paramètres $I_0, N, \beta, \gamma, M$.

Les $T_i$ sont identiquement distribuées $\forall i\in[0,N]$ car, $\forall i\in[0,N], N, \beta, \gamma, M$ sont constants.

On postule sans trop de risque que les $T_i=\inf(n\geq 0, I_0=i, I_n=0)$ sont indépendants $\forall i\in[0,N]$, puisque à aucun moment, $\forall i\neq j$, la valeur de $T_{j}$ intervient dans le calcul de $T_{i}$.

On peut donc appliquer le théorème central limite, qui nous donne que $(T_i)_ {n\in [o,h]}, i\in[0,N]$ peut être approximée par une loie normale $\mathcal{N}(N\mu,N\sigma^2)$, avec $\mu$ et $\sigma^2$ respectivement l'espérance et la variance de la loi des $(T_i)_{n\in [0,h]}$. On approxime $N\mu$ et $N\sigma^2$ par leurs estimateurs sans biais empiriques.

On peut ainsi calculer les intervalles de confiance correspondants. On se donne un seuil de confiance à $0.95$. On a alors $\forall i\in[0, N]$:
$$
E(T_i\in \Big[T_i \pm 1.96 \sqrt{\frac{\sigma^2}{M}}\Big]
$$
avec $n\in [0,h]$, $M$ le nombre de simulation faite dans Montecarlo, et $1.96=\alpha, P(|X|<\alpha)=0.95, X\sim \mathcal{N}(0,1)$.

On applique les calcules dans la fonction *temps_esp_montecarlo* et on les fait apparaître sur les graphiques précédents.