Verhindern, dass ein stochastisches Gillespie-SSA-Modell von negativ läuft

Question

Ich habe ein stochastisches Infektionsmodell (parasitischer Wurm) unter Verwendung einer Gillespie SSA hergestellt. Das Modell verwendete das Paket "GillespieSSA" (https://cran.r-project.org/web/packages/GillespieSSA/index.html). Kurz gesagt modelliert der Code eine Population von diskreten Kompartimenten. Die Bewegung zwischen den Abteilen hängt von benutzerdefinierten Ratengleichungen ab. Der SSA-Algorithmus berechnet die Anzahl der Ereignisse, die von jeder Ratengleichung für einen gegebenen Zeitschritt (tau) erzeugt werden, und aktualisiert die Population entsprechend, wobei der Prozess bis zu einem gegebenen Zeitpunkt wiederholt wird. Das Problem ist, dass die Anzahl der Ereignisse Poisson-verteilt angenommen wird (Poisson (rate [i] * tau)), so dass ein Fehler auftritt, wenn die Rate negativ ist, einschließlich wenn die Populationszahlen negativ werden.

# 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) 
} 

Somit Raten berechnet und das Modell aktualisiert die Bevölkerungswerte für jeden Zeitschritt, mit dem Datenspeicher in einem Datenrahmen, dann zusammen befestigt mit vorhergehenden Läufen. Wenn die Bevölkerungszahlen jedoch sehr niedrig sind, können Ereignisse auftreten, so dass die Anzahl der Ereignisse, die bei einer Bevölkerungsreduktion auftreten, größer ist als die Anzahl im Abteil. Eine Methode besteht darin, den Zeitschritt sehr klein zu machen, dies verlängert jedoch die Länge der Simulation sehr stark.

Meine Frage ist, gibt es eine Möglichkeit, den Code so zu erweitern, dass, wenn die Daten zu jedem Zeitschritt erstellt/berechnet werden, alle Werte von negativen Populationszahlen in 0 umgewandelt werden?

Ich habe versucht, an diesem Problem zu arbeiten, aber scheint nur in der Lage zu sein, Methoden zu entwickeln, die die Werte ändern, sobald die Simulation abgeschlossen ist, wobei die negativen Werte immer noch Probleme in den Läufen verursachen. Zum Beispiel
if (sir.out $ L < 0) sir.out $ L == 0

Jede Hilfe würde geschätzt

Answer 1

Ich glaube, das Problem ist die Methode Sie ("ETL") in der ssa gesetzt Funktion. Die ETL-Methode wird schließlich negative Zahlen erzeugen. Sie können die „OTL“ Methode versuchen, basierend auf Effiziente Schrittgröße Auswahl für die tau-springend Simulationsverfahren - in denen es noch ein paar Parameter, die Sie optimieren können, aber die grundlegende Befehl lautet:

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

Oder die direkte Methode, die überhaupt keine negative Zahl produziert:

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