我使用 Gillespie SSA 生成了一个随机感染模型(寄生虫)。该模型使用了“GillespieSSA”包 ( https://cran.r-project.org/web/packages/GillespieSSA/index.html )。简而言之,代码模拟了一群离散的隔间。隔间之间的移动取决于用户定义的速率方程。 SSA 算法用于计算给定时间步长 (tau) 内每个速率方程产生的事件数量,并相应地更新总体,过程重复到给定时间点。问题是,假设事件数是泊松分布的(泊松(rate[i]*tau)),因此当比率为负时会产生错误,包括当人口数变为负数时。# Parameter Valuessir.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 Valuessir.x0 <- c(W=20,M=10,L=0.02)# Rate Equationssir.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 evensir.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 <- 10set.seed(1)# Data Frame of outputsir.out <- data.frame(time=numeric(),W=numeric(),M=numeric(),L=numeric())# Multiple runs and combining data and SSA methodsfor(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 任何帮助,将不胜感激 最佳答案 我相信问题在于您在 ssa 函数中设置的方法(“ETL”)。 ETL 方法最终会产生负数。您可以尝试基于 Efficient step size selection for the tau-leaping 模拟方法的“OTL”方法 - 其中还有一些参数可以调整,但基本命令是: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")关于r - 防止 Gillespie SSA 随机模型运行为负值,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/37748255/
10-12 17:19