R中ODE求解器的问题

Question

假设我们在 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

    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）或两列的 Matrix，但它不起作用。我必须做什么才能完成这项工作？

Answer 1

您的代码有几个问题，一个是时间从零开始，而上限需要从一开始，并且参数名称也有些混乱。在下文中，我展示了使用approxfuns 而不是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)