Empêcher un modèle stochastique SSA Gillespie de s'exécuter Négatif

Question

J'ai produit un modèle stochastique d'infection (ver parasite), en utilisant un SSA Gillespie. Le modèle a utilisé le package "GillespieSSA" (https://cran.r-project.org/web/packages/GillespieSSA/index.html). En résumé, le code modélise une population de compartiments discrets. Le mouvement entre les compartiments dépend des équations de vitesse définies par l'utilisateur. L'algorithme SSA agit pour calculer le nombre d'événements produits par chaque équation de taux pour un pas de temps donné (tau) et met à jour la population en conséquence, le processus se répète jusqu'à un moment donné. Le problème est, le nombre d'événements est supposé Poisson distribué (Poisson (taux [i] * tau)), produit donc une erreur lorsque le taux est négatif, y compris lorsque les nombres de population deviennent négatifs.

# Parameter Values 
sir.parms <- c(deltaHinfinity=0.00299, CHi=0.00586, deltaH0=0.0854, aH=0.5, 
       muH=0.02, SigmaW=0.1, SigmaM =0.8, SigmaL=104, phi=1.15, f = 0.6674, 
       deltaVo=0.0166, CVo=0.0205, alphaVo=0.5968, beta=52, mbeta=7300 ,muV=52, g=0.0096, N=100) 
# Inital Population Values 
sir.x0 <- c(W=20,M=10,L=0.02) 
# Rate Equations 
sir.a <- c("((deltaH0+deltaHinfinity*CHi*mbeta*L)/(1+CHi*mbeta*L))*mbeta*L*N" 
      ,"SigmaW*W*N", "muH*W*N", "((1/2)*phi*f)*W*N", "SigmaM*M*N", "muH*M*N", 
      "(deltaVo/(1+CVo*M))*beta*M*N", "SigmaL*L*N", "muV*L*N", "alphaVo*M*L*N", "(aH/g)*L*N") 
# Population change for even 
sir.nu <- matrix(c(+0.01,0,0, 
        -0.01,0,0, 
        -0.01,0,0, 
        0,+0.01,0, 
        0,-0.01,0, 
        0,-0.01,0, 
        0,0,+0.01/230, 
        0,0,-0.01/230, 
        0,0,-0.01/230, 
        0,0,-0.01/230, 
        0,0,-0.01/32),nrow=3,ncol=11,byrow=FALSE) 
runs <- 10 
set.seed(1) 

# Data Frame of output 
sir.out <- data.frame(time=numeric(),W=numeric(),M=numeric(),L=numeric()) 
# Multiple runs and combining data and SSA methods 
for(i in 1:runs){ 
    sim <- ssa(sir.x0,sir.a,sir.nu,sir.parms, method="ETL", tau=1/12, tf=140, simName="SIR") 
    sim.out <- data.frame(time=sim$data[,1],W=sim$data[,2],M=sim$data[,3],L=sim$data[,4]) 

    sim.out$run <- i 
    sir.out <- rbind(sir.out,sim.out) 
} 

Ainsi, les taux sont calculés et le modèle met à jour les valeurs de la population pour chaque pas de temps, avec le magasin de données dans une trame de données, puis attachés ensemble avec des pistes précédentes. Cependant, lorsque les niveaux de la population deviennent très faibles, des événements peuvent survenir de sorte que le nombre d'événements qui se produisent en réduisant une population est plus grand que le nombre dans le compartiment. Une méthode consiste à rendre le pas de temps très petit, mais cela augmente considérablement la durée de la simulation.

Ma question existe-t-il un moyen d'augmenter le code de sorte que, lorsque les données sont créées/calculées à chaque pas de temps, les valeurs des nombres de population négatifs sont converties en 0?

J'ai essayé de travailler sur ce problème, mais je ne peux que trouver des méthodes qui modifient les valeurs une fois la simulation terminée, les valeurs négatives causant toujours des problèmes dans les exécutions elles-mêmes. E.g.
si (sir.out $ L < 0) sir.out $ L == 0

Toute aide serait appréciée

Answer 1

Je crois que le problème est la méthode que vous définissez ("ETL") dans la région visée fonction. La méthode ETL produira éventuellement des nombres négatifs. Vous pouvez essayer la méthode « OTL », basé sur sélection de la taille de l'étape efficace pour la méthode de simulation tau-saut - dans lequel il y a quelques paramètres que vous pouvez modifier, mais la commande de base est la suivante:

ssa(sir.x0,sir.a,sir.nu,sir.parms, method="OTL", tf=140, simName="SIR") 

Ou la méthode directe, qui ne produira pas de nombre négatif que ce soit:

ssa(sir.x0,sir.a,sir.nu,sir.parms, method="D", tf=140, simName="SIR")