使用ODE中的循环以图形方式比较不同的参数R.

Question

我正在使用该deSolve软件包绘制几个微分方程(如果感兴趣,请阅读http://www.maa.org/press/periodicals/loci/joma/the-sir-model-for-spread-of-disease-the-微分方程模型).

我最终的目标是创建一个迭代函数或过程(for循环)来绘制某些参数(beta和gamma)的变化将如何影响解决方案.首选输出是一个列表,其中包含ode循环中每个指定beta值的所有解决方案.我遇到了将循环集成到deSolve软件包所需的设置中的问题ode.

在下面的代码中,我试图绘制参数beta中值的范围(1到2的增量为0.1)将如何影响微分方程的图.

for(k in seq(1,2,by=0.1)){ #range of values for beta

    init <- c(S=1-1e-6, I=1e-6, R=0) #initial conditions for odes
    time <- seq(0,80,by=1) #time period
    parameters <- c(beta=k, gamma=0.15) #parameters in ode

SIR <- function(time,state,parameters){ #function containing equaations
    with(as.list(c(state,parameters)),{
        dS <- -beta*S*I
        dI <- beta*S*I-gamma*I
        dR <- gamma*I

        return(list(c(dS,dI,dR)))
    })
}

ode(y=init,times=time,func=SIR()[beta],parms=parameters[k])}

}

我得到的第一个错误表明缺少SIR函数中的参数参数

as.list中的错误(c(init,parameters)):缺少参数"parameters",没有默认值

我不明白为什么在我之前分配parameters的行中报告了这个错误.

Answer 1

您也可以在循环外定义渐变函数(以及其他非变化元素):

SIR <- function(time,state,parameters) {
  with(as.list(c(state,parameters)),{
    dS <- -beta*S*I
    dI <- beta*S*I-gamma*I
    dR <- gamma*I
    return(list(c(dS,dI,dR)))
  })
}
init <- c(S=1-1e-6, I=1e-6, R=0) #initial conditions for odes
time <- seq(0,80,by=1) #time period

现在定义要尝试的值向量(不必要但方便):

betavec <- seq(1,2,by=0.1)

并定义一个列表来保存结果:

res <- vector(length(betavec),mode="list")

library(deSolve)
for (k in seq_along(betavec)){ #range of values for beta
    res[[k]] <- ode(y=init,times=time,func=SIR,
           parms=c(beta=betavec[k], gamma=0.15))
}

现在您有一个列表,其中每个元素都包含一次运行的结果.您可以sapply或lapply在此列表上,例如,从每次运行中获取最后状态的矩阵:

t(sapply(res,tail,1))

或者,如果您希望将结果作为一个长数据框...

names(res) <- betavec  ## to get beta value incorporated in results
dd <- dplyr::bind_rows(lapply(res,as.data.frame),.id="beta")
dd$beta <- as.numeric(dd$beta)

do.call(rbind,...)会像以前一样工作bind_rows(),但是bind_rows.id参数可以方便地将beta值添加到每个数据框.您还可以将结果作为列表保留并循环显示它们,同时使用单独的lines()调用进行绘图,或者(例如)将感染列绑定在一起并使用matplot()它们同时绘制它们.这只是风格和成语的问题.

library(ggplot2); theme_set(theme_bw())
library(viridis)
ggplot(dd,aes(x=time,y=I,colour=beta))+
    geom_line(aes(group=beta))+
    scale_color_viridis()+
    scale_y_log10()