3

我使用 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

任何帮助，将不胜感激

4

1 回答 1

2

我相信问题在于您在 ssa 函数中设置的方法（“ETL”）。ETL 方法最终会产生负数。您可以尝试“OTL”方法，该方法基于tau-leaping 模拟方法的有效步长选择- 您可以在其中调整更多参数，但基本命令是：

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")
于 2016-09-13T15:53:31.697 回答