我试图使用deSolve拟合一个合理的基本ODE模型,包括在模型中一个随时间变化的参数(感染力; FOI)。在没有此参数的情况下运行模型时,包含时间相关参数会产生错误(参见下文)。
我对R和数学建模相对较新,并且一直试图解决这个问题。
我已经将FOI参数创建为值矩阵,然后使用approxfun函数进行插值(我已经看到它可以使用强制函数,例如https://rdrr.io/rforge/deSolve/man/forcings.html)。
没有此时间相关参数的模型运行时没有任何错误,但尝试包含它会产生错误:
Error in checkFunc(Func2, times, y, rho) :
The number of derivatives returned by func() (200) must equal the
length of the initial conditions vector (2)
我无法弄清楚如何解决这个错误,因为我只有2个初始条件,似乎包含这个时间相关的FOI参数会产生更多的衍生物。
我知道其他人已经提出了类似的问题,但我没有发现这个强迫功能的问题。
非常感谢任何建议。
# Forcing function data
foi <- matrix(ncol=2,byrow=TRUE,data=c(
0, 0.003, 2, 0.03, 3, 0.08, 4,0.1, 5, 0.12, 6, 0.15,
8, 0.16, 10, 0.14,12, 0.12,14,0.08,15, 0.06,16, 0.03,
17, 0.01,18,0.003,19,0.003,20,0.003,30,0.003,40,0.003,
50,0.003,60,0.003,65,0.01, 70,0.08, 72,0.095,74,0.10,
76,0.1, 78,0.08, 80,0.06))
age <- seq(0, 80, by = 1)
input <- approxfun(x = foi[,1], y = foi[,2], method = "constant", rule = 2)
# Function
ab <- function(time, state, pars) {
with(as.list(c(state, pars)), {
import<-c(input(t))
diggP<- (import *iggN) - iggR*iggP
diggN<- (-import*iggN) + iggR*iggP
return(list(c(diggP, diggN)))
})
}
# Initial values
yini <- c(iggP=0, iggN=1)
# Parameters
pars <- c(iggR = 0, import)
# ODE solver
results<- ode(y=yini, times=age, func=foi_model, pars)
我希望制作一个模型,其中在每个时间点(或在这种情况下年龄),FOI根据我在FOI矩阵中输入的值而变化。因此,我想看看FOI随年龄的变化如何影响微分方程的输出。
你的主要问题是你将一个参数t
传递给input
,但是你的代码中不存在该变量。时间作为一个名为time
的参数传递给您的模型。 (另外,你的模型被称为ab
而不是foi_model
,正如在调用ode
中所述,加上pars
不需要import
,应该传递给ode
。)
# Load library
library(deSolve)
# Create FOI matrix
foi <- matrix(ncol=2,byrow=TRUE,data=c(
0, 0.003, 2, 0.03, 3, 0.08, 4,0.1, 5, 0.12, 6, 0.15,
8, 0.16, 10, 0.14,12, 0.12,14,0.08,15, 0.06,16, 0.03,
17, 0.01,18,0.003,19,0.003,20,0.003,30,0.003,40,0.003,
50,0.003,60,0.003,65,0.01, 70,0.08, 72,0.095,74,0.10,
76,0.1, 78,0.08, 80,0.06))
# Times for model solution
age <- seq(0, 80, by = 1)
# Linear interpolation function from FOI data
input <- approxfun(x = foi[,1], y = foi[,2], method = "constant", rule = 2)
# Model to be integrated
ab <- function(time, state, parms) {
with(as.list(c(state, parms)), {
##### IMPORTANT #####
import<-input(time) #<- 'time' was previously 't'
#####################
# Derivatives
diggP<- (import *iggN) - iggR*iggP
diggN<- (-import*iggN) + iggR*iggP
# Return results
return(list(c(diggP, diggN)))
})
}
# Initial values
yini <- c(iggP=0, iggN=1)
# Parameters
pars <- c(iggR = 0)
# Solve model
results<- ode(y=yini, times=age, func=ab, parms = pars)
# Plot results
plot(results)
由reprex package创建于2019-03-27(v0.2.1)