r - 防止 Gillespie SSA 随机模型运行为负值

Question

我使用 Gillespie SSA 制作了一个随机感染模型(寄生虫)。该模型使用“GillespieSSA”包( https://cran.r-project.org/web/packages/GillespieSSA/index.html )。简而言之，该代码模拟了一组离散的隔间。隔室之间的移动取决于用户定义的速率方程。 SSA 算法用于计算给定时间步 (tau) 内每个速率方程产生的事件数，并相应地更新总体，过程重复到给定时间点。问题是，事件数量假设为泊松分布(Poisson(rate[i]*tau))，因此当比率为负时(包括人口数量变为负数时)会产生错误。

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

因此，计算比率，模型更新每个时间步的总体值，数据存储在数据框中，然后与之前的运行连接在一起。然而，当人口水平变得非常低时，可能会发生事件，使得减少人口的事件数量大于隔间中的数量。一种方法是使时间步非常小，但这大大增加了模拟的长度。

我的问题是有没有一种方法可以增强代码，以便在每个时间步创建/计算数据时，任何负数的人口数值都会转换为 0？

我已经尝试解决这个问题，但似乎只能想出在模拟完成后改变值的方法，负值仍然会导致运行本身出现问题。例如
if (sir.out$L < 0) Sir.out$L == 0

如有任何帮助，我们将不胜感激

Answer 1

我相信问题出在您在 ssa 函数中设置的方法(“ETL”)。 ETL方法最终会产生负数。您可以尝试“OTL”方法，基于 tau 跳跃模拟方法的有效步长选择 - 其中还有一些参数可以调整，但基本命令是:

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

或者直接方法，它不会产生任何负数:

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