防止 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？

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

任何帮助将不胜感激

Answer 1

我认为问题在于您在 ssa 函数中设置的方法 ("ETL")。 ETL 方法最终会产生负数。您可以尝试 "OTL" 方法，基于 Efficient step size selection for the tau-leaping simulation method- 其中还有一些您可以调整的参数，但是基本命令是：

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")

防止 Gillespie SSA 随机模型出现运行负值

Preventing a Gillespie SSA Stochastic Model From Running Negative

modeling

r

stochastic-process

stochastic

防止 Gillespie SSA 随机模型出现 运行 负值

Preventing a Gillespie SSA Stochastic Model From Running Negative

modeling

r

stochastic-process

stochastic

防止 Gillespie SSA 随机模型出现运行负值