R 中 ODE 求解器的问题
Problems with ODE solver in R
假设我们在 R 中有一个任意的 ODE 系统,我们想要求解它,例如 SIR 模型
dS <- -beta * I * S
dI <- beta * I * S - gamma * I
dR <- gamma * I
我希望 beta 和 gamma 具有随时间变化的参数,例如
beta_vector <- seq(0.05, 1, by=0.05)
gamma_vector <- seq(0.05, 1, by=0.05)
用户@Ben Bolker 给了我在梯度函数中使用 beta <- beta_vector[ceiling(time)] 的建议
sir_1 <- function(beta, gamma, S0, I0, R0, times) {
require(deSolve) # for the "ode" function
# the differential equations:
sir_equations <- function(time, variables, parameters) {
beta <- beta_vector[ceiling(time)]
gamma <- gamma_vector[ceiling(time)]
with(as.list(c(variables, parameters)), {
dS <- -beta * I * S
dI <- beta * I * S - gamma * I
dR <- gamma * I
return(list(c(dS, dI, dR)))
})
}
# the parameters values:
parameters_values <- c(beta=beta, gamma = gamma)
# the initial values of variables:
initial_values <- c(S = S0, I = I0, R = R0)
# solving
out <- ode(initial_values, times, sir_equations, parameters_values)
# returning the output:
as.data.frame(out)
}
sir_1(beta = beta, gamma = gamma, S0 = 99999, I0 = 1, R0 = 0, times = seq(0, 19))
当我执行它时,出现以下错误
Error in checkFunc(Func2, times, y, rho) :
The number of derivatives returned by func() (1) must equal the length of the initial
conditions vector (3)
问题一定出在这里:
parameters_values <- c(beta=beta, gamma = gamma)
我曾尝试将 paramters_values 更改为具有两行(第一行是 beta,第二行是 gamma)或两列的矩阵,但没有成功。我需要做什么才能完成这项工作?
你的代码有几个问题,一个是时间从零开始,而天花板需要从一开始,还有一些参数名称混淆。在下文中,我展示了使用 approxfun
s 而不是 ceiling
的(几种)可能方法之一。这更稳健,即使 ceiling
也有一些优势。然后参数是作为列表传递给 ode
的函数。一种更简单的方法是使用全局变量。
另一个需要考虑的问题是时间相关的 gamma
和 beta
应该线性插值还是逐步插值。 approxfun
函数允许两者,下面我使用线性插值。
require(deSolve) # for the "ode" function
beta_vector <- seq(0.05, 1, by=0.05)
gamma_vector <- seq(0.05, 1, by=0.05)
sir_1 <- function(f_beta, f_gamma, S0, I0, R0, times) {
# the differential equations
sir_equations <- function(time, variables, parameters) {
beta <- f_beta(time)
gamma <- f_gamma(time)
with(as.list(variables), {
dS <- -beta * I * S
dI <- beta * I * S - gamma * I
dR <- gamma * I
# include beta and gamma as auxiliary variables for debugging
return(list(c(dS, dI, dR), beta=beta, gamma=gamma))
})
}
# time dependent parameter functions
parameters_values <- list(
f_beta = f_beta,
f_gamma = f_gamma
)
# the initial values of variables
initial_values <- c(S = S0, I = I0, R = R0)
# solving
# return the deSolve object as is, not a data.frame to ake plotting easier
out <- ode(initial_values, times, sir_equations, parameters)
}
times <- seq(0, 19)
# approxfun is a function that returns a function
f_gamma <- approxfun(x=times, y=seq(0.05, 1, by=0.05), rule=2)
f_beta <- approxfun(x=times, y=seq(0.05, 1, by=0.05), rule=2)
# check how the approxfun functions work
f_beta(5)
out <- sir_1(f_beta=f_beta, f_gamma=f_gamma, S0 = 99999, I0 = 1, R0 = 0, times = times)
# plot method of class "deSolve", plots states and auxilliary variables
plot(out)
假设我们在 R 中有一个任意的 ODE 系统,我们想要求解它,例如 SIR 模型
dS <- -beta * I * S
dI <- beta * I * S - gamma * I
dR <- gamma * I
我希望 beta 和 gamma 具有随时间变化的参数,例如
beta_vector <- seq(0.05, 1, by=0.05)
gamma_vector <- seq(0.05, 1, by=0.05)
用户@Ben Bolker 给了我在梯度函数中使用 beta <- beta_vector[ceiling(time)] 的建议
sir_1 <- function(beta, gamma, S0, I0, R0, times) {
require(deSolve) # for the "ode" function
# the differential equations:
sir_equations <- function(time, variables, parameters) {
beta <- beta_vector[ceiling(time)]
gamma <- gamma_vector[ceiling(time)]
with(as.list(c(variables, parameters)), {
dS <- -beta * I * S
dI <- beta * I * S - gamma * I
dR <- gamma * I
return(list(c(dS, dI, dR)))
})
}
# the parameters values:
parameters_values <- c(beta=beta, gamma = gamma)
# the initial values of variables:
initial_values <- c(S = S0, I = I0, R = R0)
# solving
out <- ode(initial_values, times, sir_equations, parameters_values)
# returning the output:
as.data.frame(out)
}
sir_1(beta = beta, gamma = gamma, S0 = 99999, I0 = 1, R0 = 0, times = seq(0, 19))
当我执行它时,出现以下错误
Error in checkFunc(Func2, times, y, rho) :
The number of derivatives returned by func() (1) must equal the length of the initial
conditions vector (3)
问题一定出在这里:
parameters_values <- c(beta=beta, gamma = gamma)
我曾尝试将 paramters_values 更改为具有两行(第一行是 beta,第二行是 gamma)或两列的矩阵,但没有成功。我需要做什么才能完成这项工作?
你的代码有几个问题,一个是时间从零开始,而天花板需要从一开始,还有一些参数名称混淆。在下文中,我展示了使用 approxfun
s 而不是 ceiling
的(几种)可能方法之一。这更稳健,即使 ceiling
也有一些优势。然后参数是作为列表传递给 ode
的函数。一种更简单的方法是使用全局变量。
另一个需要考虑的问题是时间相关的 gamma
和 beta
应该线性插值还是逐步插值。 approxfun
函数允许两者,下面我使用线性插值。
require(deSolve) # for the "ode" function
beta_vector <- seq(0.05, 1, by=0.05)
gamma_vector <- seq(0.05, 1, by=0.05)
sir_1 <- function(f_beta, f_gamma, S0, I0, R0, times) {
# the differential equations
sir_equations <- function(time, variables, parameters) {
beta <- f_beta(time)
gamma <- f_gamma(time)
with(as.list(variables), {
dS <- -beta * I * S
dI <- beta * I * S - gamma * I
dR <- gamma * I
# include beta and gamma as auxiliary variables for debugging
return(list(c(dS, dI, dR), beta=beta, gamma=gamma))
})
}
# time dependent parameter functions
parameters_values <- list(
f_beta = f_beta,
f_gamma = f_gamma
)
# the initial values of variables
initial_values <- c(S = S0, I = I0, R = R0)
# solving
# return the deSolve object as is, not a data.frame to ake plotting easier
out <- ode(initial_values, times, sir_equations, parameters)
}
times <- seq(0, 19)
# approxfun is a function that returns a function
f_gamma <- approxfun(x=times, y=seq(0.05, 1, by=0.05), rule=2)
f_beta <- approxfun(x=times, y=seq(0.05, 1, by=0.05), rule=2)
# check how the approxfun functions work
f_beta(5)
out <- sir_1(f_beta=f_beta, f_gamma=f_gamma, S0 = 99999, I0 = 1, R0 = 0, times = times)
# plot method of class "deSolve", plots states and auxilliary variables
plot(out)