在R中寻找SIR模型参数的问题

我正在尝试使用 R 中的 SIR 模型来模拟冠状病毒感染率。但是，如您所见，我的 beta（控制易感者和感染者之间的过渡）和 gamma（控制感染者和康复者之间的过渡）值不正确。

      beta     gamma 
 1.0000000 0.8407238

这是我的代码：

library(deSolve)


Infected <- c(994, 307, 329, 553, 587, 843, 983, 1750, 2950, 4569, 5632, 4848, 9400, 10311, 11166, 13451, 17388, 18743, 19452, 20065, 20732,24914)
Day <- 1:(length(Infected))
N <- 331002651 # population of the us

这是SIR功能

SIR <- function(time, state, parameters) {
  par <- as.list(c(state, parameters))
  with(par, {
    dS <- -beta/N * I * S
    dI <- beta/N * I * S - gamma * I
    dR <- gamma * I
    list(c(dS, dI, dR))
  })
}

然后我找到RSS

init <- c(S = N-Infected[1], I = Infected[1], R = 0)
RSS <- function(parameters) {
  names(parameters) <- c("beta", "gamma")
  out <- ode(y = init, times = Day, func = SIR, parms = parameters)
  fit <- out[ , 3]
  sum((Infected - fit)^2)
}

Opt <- optim(c(0.5, 0.5), RSS, method = "L-BFGS-B", lower = c(0, 0), upper = c(1, 1)) # optimize with some sensible conditions
Opt$message

然后我找到 beta 和 gamma

Opt_par <- setNames(Opt$par, c("beta", "gamma"))
Opt_par
##      beta     gamma 
## 1.0000000 0.8407238

我的代码出了什么问题以及如何修复它？先感谢您！

优化通常很棘手，有几件事可以改进：

• 模型非常简单，您可以考虑使用SEIR模型代替（仍然很简单），
• 扩大盒子限制，这对你的情况有帮助，
• 使用另一个优化器，例如来自像这样的元包FME，还有额外的诊断，
• 将变量重新调整到更方便的范围，例如将人口除以 1e3，
• 或者，作为替代方案，调整公差（atol）适合原始尺度。

library(deSolve)

Infected <- c(994, 307, 329, 553, 587, 843, 983, 1750, 2950, 4569, 5632, 4848, 
              9400, 10311, 11166, 13451, 17388, 18743, 19452, 20065, 20732,24914)
Day <- 1:(length(Infected))
N <- 331002651 # population of the us

SIR <- function(time, state, parameters) {
  par <- as.list(c(state, parameters))
  with(par, {
    dS <- -beta/N * I * S
    dI <- beta/N * I * S - gamma * I
    dR <- gamma * I
    list(c(dS, dI, dR))
  })
}

init <- c(S = N-Infected[1], I = Infected[1], R = 0)
RSS <- function(parameters) {
  names(parameters) <- c("beta", "gamma")
  out <- ode(y = init, times = Day, func = SIR, parms = parameters)
  fit <- out[ , 3]
  sum((Infected - fit)^2)
}

# optimize with some sensible conditions
Opt <- optim(c(0.5, 0.5), RSS, method = "L-BFGS-B", 
             lower = c(0, 0), upper = c(10, 10))
Opt$message

Opt_par <- setNames(Opt$par, c("beta", "gamma"))
Opt_par
# beta     gamma 
# 0.6512503 0.4920399 

out <- ode(y = init, times = Day, func = SIR, parms = Opt_par)

plot(out)
plot(out, obs=data.frame(time=Day, I=Infected))

可以找到基本的 SEIR 模型和其他链接here https://github.com/tpetzoldt/covid。但请注意，不仅模型太简单，数据也太简单。确诊病例数为not感染人数，取决于检测方式，可能比感染人数少得多。