Petit cours d'épidémiologie mathématique
Simulation des systèmes stochastiques

Julien Arino

Department of Mathematics & Data Science Nexus
University of Manitoba*

Centre canadien de modélisation des maladies (CCDM/CCMM)
NSERC-PHAC EID Modelling Consortium (CANMOD, MfPH, OMNI/RÉUNIS)

* The University of Manitoba campuses are located on original lands of Anishinaabeg, Cree, Oji-Cree, Dakota and Dene peoples, and on the homeland of the Métis Nation.

Plan de ce cours

  • Simulation et analyse numérique des CMTD
  • Simulation des CMTC
  • Paralléliser son code

Remarques / Ressources

Ceci est un cours orienté sur l'utilisation: je touche à peine aux algorithmes, je me concentre sur leur utilisation

Le code est disponible dans ce répo Github dans le répertoire CODE

Simulation et analyse numérique des CMTD

Le SIS DTMC

On va considérer un modèle SIS en population totale constante . le nombre de susceptibles, le nombre d'infectieux. Les transitions suivantes sont possibles entre et :

  • Nouvelle infection, avec probabilité
  • Guérison, avec probabilité
  • Pas de transition: complément des deux probabilités précédentes,

Puisque , on ne considère que les infectés

Le SIS en marche aléatoire (version simplifiée)

Transitions:

  • Nouvelle infection
  • Guérison
  • Pas de transition

On peut donc considérer une marche aléatoire absorbante sur

Attention à !!!

Pour , on va calculer et (pour , on sait que toute la masse est sur puisque est absorbant). On veut

puisque l'on doit avoir des probabilités. n'est pas un problème, puisqu'on rattrape les choses avec

Donc on doit avoir

Calcul de la borne sup de

On a:

  • (pourvu que ..)
  • , avec

Matrice de transition

donnée par

# Construction de la matrice de transition
T = mat.or.vec(nr = (Pop+1), nc = (Pop+1))
for (row in 2:Pop) {
  I = row-1
  mv_right = gamma*I*Delta_t # Guérisons
  mv_left = beta*I*(Pop-I)*Delta_t # Infections
  T[row,(row-1)] = mv_right
  T[row,(row+1)] = mv_left
}
# Dernier rang, mouvement vers la gauche uniquement
T[(Pop+1),Pop] = gamma*(Pop)*Delta_t
# On vérifie que les valeurs ne sont pas trop grandes
if (max(rowSums(T))>1) {
  T = T/max(rowSums(T))
}
diag(T) = 1-rowSums(T)

Simulation d'une CMTD

library(DTMCPack)
sol = MultDTMC(nchains = 20, tmat = T, io = IC, n = t_f)

donne 20 réalisations d'une CMTD avec matrice de transition T, conditions initiales IC (un vecteur de probabilité d'être dans les differents états de ) et temps final t_f

Voir le code

Aller plus loin

DTMCPack est bien pour obtenir des réalisation d'une CMTD, mais pour l'étudier plus en détail, markovchain est bien plus complet

library(markovchain)
mcSIS <- new("markovchain", 
             states = sprintf("I_%d", 0:Pop),
             transitionMatrix = T,
             name = "SIS")

De façon intéressante, markovchain redéfinit l'étrange notation de Rpour les produits matriciels où "* est Hadamard, %*% est l'habituel", si bien que pour un objet défini comme markovchain, mcSIS*mcSIS est le produit matriciel usuel

> summary(mcSIS)
SIS  Markov chain that is composed by: 
Closed classes: 
I_0 
Recurrent classes: 
{I_0}
Transient classes: 
{I_1,I_2,I_3,I_4,I_5,I_6,I_7,I_8,I_9,I_10,I_11,I_12,I_13,I_14,I_15,
I_16,I_17,I_18,I_19,I_20,I_21,I_22,I_23,I_24,I_25,I_26,I_27,I_28,
I_29,I_30,I_31,I_32,I_33,I_34,I_35,I_36,I_37,I_38,I_39,I_40,I_41,
I_42,I_43,I_44,I_45,I_46,I_47,I_48,I_49,I_50,I_51,I_52,I_53,I_54,
I_55,I_56,I_57,I_58,I_59,I_60,I_61,I_62,I_63,I_64,I_65,I_66,I_67,
I_68,I_69,I_70,I_71,I_72,I_73,I_74,I_75,I_76,I_77,I_78,I_79,I_80,
I_81,I_82,I_83,I_84,I_85,I_86,I_87,I_88,I_89,I_90,I_91,I_92,I_93,
I_94,I_95,I_96,I_97,I_98,I_99,I_100}
The Markov chain is not irreducible 
The absorbing states are: I_0
Function Role
absorbingStates absorbing states of the transition matrix, if any
steadyStates the vector(s) of steady state(s) in matrix form
meanFirstPassageTime matrix or vector of mean first passage times
meanRecurrenceTime vector of mean number of steps to return to each recurrent state
hittingProbabilities matrix of hitting probabilities for a Markov chain
meanAbsorptionTime expected number of steps for a transient state to be absorbed by any recurrent class
absorptionProbabilities transient states being absorbed by each recurrent state
canonicForm canonic form of transition matrix
period the period of an irreducible DTMC
summary DTMC summary

Simulation des CMTC

  • Algorithme de Gillespie
  • Utilisation d'une librairie
  • Limitation de Gillespie et tau-leap

Algorithme de Gillespie

Poids Transition Effet
, nouvelle infection d'un susceptible
, guérison d'un infectieux

On utilise et on omet la dynamique de

Algorithme de Gillespie (modèle SIS avec seulement I)

posons et
tant que {}

  • tirer de
  • tirer de
  • trouver t.q.
  • selon {}
    • 1: nouvelle infection,
    • 2: fin de la période infectieuse,

Utilisation d'une librairie

  • Vous pouvez implémenter l'algorithme de Gillespie à la main, c'est un bon exercice. On le fait plus loin, d'ailleurs
  • Mais dans R, il y a plusieurs librairies permettant de faire ça facilement
  • Ces librairies ont aussi l'avantage qu'elles implémentent le tau-leap (voir plus loin dans ce cours)

Simulation d'une CMTC avec une librairie

library(GillespieSSA2)
IC <- c(S = (Pop-I_0), I = I_0)
params <- c(gamma = gamma, beta = beta)
reactions <- list(
  reaction("beta*S*I", c(S=-1,I=+1), "new_infection"),
  reaction("gamma*I", c(S=+1,I=-1), "recovery")
)
set.seed(NULL)
sol <- ssa(
    initial_state = IC,
    reactions = reactions,
    params = params,
    method = ssa_exact(),
    final_time = t_f,
)
plot(sol$time, sol$state[,"I"], type = "l",
     xlab = "Time (days)", ylab = "Number infectious")

Limitation de Gillespie et tau-leap

Parfois, simuler une CMTC peut être compliqué

  • Le temps inter-évènements est distribué exponentiellement
  • Étape critique de l'algorithme de Gillespie:
    • poids de tous les évènements possibles (propensité)
    • Tirer au hasard d'une
  • Donc le temps inter-évènement si devient très grand (Cours 05: exponentielle de paramètre a un temps de séjour moyen )
  • Ceci peut causer un "arrêt" de la simulation

Exemple: un processus de naissance et mort

  • Individus naissent au taux per capita
  • Individus meurrent au taux per capita
  • Faisons un Gillespie classique "à la main"

Algorithme de Gillespie (modèle naissance & mort)

posons et
tant que {}

  • tirer de
  • tirer from
  • trouver t.q.
  • selon {}
    • 1: naissance,
    • 2: mort, 
b = 0.01   # Taux de naissance
d = 0.01   # Taux de mortalité
t_0 = 0    # Temps initial
N_0 = 100  # Population initiale

# Vecteurs pour le temps et l'état, initialisés avec la condition initiale
t = t_0
N = N_0

t_f = 1000  # Temps final

# On garde le temps et l'état actuel en mémoire. On pourrait juste regarder 
# la fin des vecteurs t et N, mais ça demande plus d'opérations
t_curr = t_0
N_curr = N_0

while (t_curr<=t_f) {
  xi_t = (b+d)*N_curr
  # Pour éviter les problèmes lorsque N=0, on termine si c'est le cas
  if (N_curr == 0) {
    break
  }
  tau_t = rexp(1, rate = xi_t)
  t_curr = t_curr+tau_t
  v = c(b*N_curr, xi_t)/xi_t
  zeta_t = runif(n = 1)
  pos = findInterval(zeta_t, v)+1
  switch(pos,
         { 
           # Naissance
           N_curr = N_curr+1
         },
         {
           # Mort
           N_curr = N_curr-1
         })
  N = c(N, N_curr)
  t = c(t, t_curr)
}

Choisir au hasard une valeur d'une loi exponentielle

Si vous n'avez pas accès à un générateur de nombre aléatoires exponentiels (pas de problème en R), vous pouvez utiliser une loi uniforme sur (ça, ça existe toujoura)

On veut de , i.e., a la densité de probabilité

On utilise l'inverse de la densité cumulative

qui à ses valeurs dans . Tirons de et résolvons pour

Ce dernier ne s'est pas bien passé

  • Je voulais 1000 unités de temps (jours?)
  • J'ai interrompu lorsque parce que j'ai perdu patience
  • (Transparent : la simul s'est arrêtée parce que la population s'est éteinte, je ne l'ai pas arreêté!)
  • Au temps d'arrêt
    • le temps allait très lentement
> tail(diff(t))
[1] 1.357242e-04 1.291839e-04 5.926044e-05 7.344020e-05 1.401148e-04 4.423529e-04

Le tau-leap au secours !

  • Je ne vais pas rentrer dans les détails
  • Méthode d'approximation (pas comme Gillespie classique, qui est exact)
  • En gros: considère des "groupes" d'évènements plutôt que des évènements individuels
  • Bonne nouvelle: GillespieSSA2 (qu'on a vu avant) et adaptivetau implémentent le tau-leap

Paralléliser son code

Parallélisation

  • Pour faire beaucoup de simulations (voir de multiples réalisations), c'est une bonne idée de paralléliser son code
  • Les réalisations d'une chaîne de Markov (ainsi que bien d'autres problèmes dans ce petit cours) sont ce que l'on appelle embarrassingly parallel: les réalisations sont indépendantes, donc pourvu d'avoir assez de CPU et assez de RAM, on peut les faire tourner en même temps sans se poser de questions
  • En R, on utilisera parLapply pour ce genre de simulations
  • Implémentation facile, mais encore plus facile avec un exemple :)

On écrit une fonction, e.g., run_one_sim qui .. fait une simulation, puis on l'appelle avec parLapply après avoir "instancié" un cluster

no_cores <- detectCores()-1
cl <- makeCluster(no_cores)
clusterEvalQ(cl,{
  library(GillespieSSA2)
})
clusterExport(cl,
              c("params", "run_one_sim"),
              envir = .GlobalEnv)
SIMS = parLapply(cl = cl, 
                 X = 1:params$number_sims, 
                 fun =  function(x) run_one_sim(params))
stopCluster(cl)

Voir simulate_CTMC_parallel.R sur Github

Bénéfices de la parallélisation

À titre d'exemple, pour montrer à quel point on peut gagner en temps

On fait tourner le code parallèle pour 100 simuls entre tictoc::tic() et tictoc::toc(), donnant 66.958 sec elapsed, puis la version séquentielle

tictoc::tic()
SIMS = lapply(X = 1:params$number_sims, 
              FUN =  function(x) run_one_sim(params))
tictoc::toc()

qui donne 318.141 sec elapsed sur une machine Intel Corei9-8950HK 6C/12T @ 2.90GHz (4.75 plus rapide). Sur une machine plus orientée calcul: 12.067 sec elapsed versus 258.985 sec elapsed sur un processeur AMD Ryzen Threadripper 3970X 32C/64T @ 4.15GHz (21.46 plus rapide !)

Remarque: interpolation des solutions

  • Dans la figure sur le transparent précédent, on présente la moyenne et la moyenne conditionnée sur la non extinction
  • C'est un peu plus compliqué qu'il n'y paraît: chaque réalisation a un ensemble d'instants de saut différent, puisque ces temps sont tirés au hasard
  • Pour calculer une moyenne, on doit donc interpoler les solutions
  • Dans le code, c'est ce que fait la fonction approx, que l'on utilise sur la solution générée par la fonction ssa