利用rjags估计微分方程组的参数

2023-08-31  本文已影响0人  小潤澤

项目地址：https://github.com/lilywang1988/eSIR

本项目利用前30d已知的数据先估计微分方程组的各项参数，然后利用估计好的参数来估计后170d各个状态病人的人数走势（估计时，以前30d最后1d的各个状态病人的人数作为起点，估计后170d各个状态病人的人数走势）

安装必要的环境：

这里需要安装 R 包 rjags 和对应的可执行文件JAGS，JAGS的下载地址为：https://sourceforge.net/projects/mcmc-jags/files/JAGS/4.x/Windows/，根据自己的R版本选择合适的JAGS版本，本案例下载的是JAGS-4.3.0.exe

安装时取消32位即可：


安装路径可以自定义，本例子为：


安装rjags：

Sys.setenv(JAGS_HOME='D:/tools/JAGS-4.3.0')
install.packages('rjags')
library(rjags)

微分方程模型以及 RK4 估计数值

eSIR微分方程组:


其中:

  1. θtS, 代表易感患者 S 的人数
  2. θtI, 代表感染者 I 的人数
  3. θtR, 代表恢复的人群 R
  4. β 为 transmission rate
  5. γ 为 治愈速率

作者利用Runge-Kutta(RK4)来估计微分方程的解:
Step 1:写出三个状态的迭代关系式子

公式1
变量分别如下：
公式2
公式3
公式4

1.对于一阶微分方程来说RK4的迭代估计如下：

RK4
  1. m 为迭代次数
  2. h 为步长
  3. 其中 x 代表自变量的估计区间，x0 < x < xn，x(m) 代表第 m 次迭代的时自变量 x 的值，例如步长 h = 0.1，x(0) = 0，x(1) = 0.1，x(2) = 0.2

2.对于一阶微分方程组来说RK4的迭代估计如下：

  1. k 为迭代次数
  2. 其中 x1，x2，...，xm 代表微分方程组不同的函数变量，t 代表 关于函数 x 的自变量，如果 f() 中含有 t 变量，则在迭代中需要加入 （tk，tk+h/2，tk+h）这几项变量。同理如果 f() 中含有 x1，x2，...，xm 这些变量，则在迭代中需要加入对应的变量，如果没有则不需要加入；例如没有 x2 那么后续迭代的 f() 中无需加入（x2k，x2k+h/2•f，x2k+h•f）这几项变量
  3. h 为步长

微分方程模型以及 RK4 估计数值的代码部分

# 定义RK4估计式子
## T_prime 为步长
model1.string <- paste0("
  model{
     for(t in 2:(T_prime+1)){
      # 公式2, 各个参数的意义见下面, pi 为一列数值为 1 的向量
       Km[t-1,1] <- -beta*pi[t-1]*theta[t-1,1]*theta[t-1,2]
       Km[t-1,2] <- -beta*pi[t-1]*(theta[t-1,1]+0.5*Km[t-1,1])*(theta[t-1,2]+0.5*Km[t-1,5])
       Km[t-1,3] <- -beta*pi[t-1]*(theta[t-1,1]+0.5*Km[t-1,2])*(theta[t-1,2]+0.5*Km[t-1,6])
       Km[t-1,4] <- -beta*pi[t-1]*(theta[t-1,1]+Km[t-1,3])*(theta[t-1,2]+Km[t-1,7])

       # 公式3, 各个参数的意义见下面       
       Km[t-1,5] <- -Km[t-1,1]-Km[t-1,9]
       Km[t-1,6] <- -Km[t-1,2]-Km[t-1,10]
       Km[t-1,7] <- -Km[t-1,3]-Km[t-1,11]
       Km[t-1,8] <- -Km[t-1,4]-Km[t-1,12]
       
       # 公式4, 各个参数的意义见下面
       Km[t-1,9] <- gamma*theta[t-1,2]
       Km[t-1,10] <- gamma*(theta[t-1,2]+0.5*Km[t-1,5])
       Km[t-1,11] <- gamma*(theta[t-1,2]+0.5*Km[t-1,6])
       Km[t-1,12] <- gamma*(theta[t-1,2]+Km[t-1,7])
      
       # 公式1, 各个参数的意义见下面
       alpha[t-1,1] <- theta[t-1,1]+(Km[t-1,1]+2*Km[t-1,2]+2*Km[t-1,3]+Km[t-1,4])/6
       alpha[t-1,2] <- theta[t-1,2]+(Km[t-1,5]+2*Km[t-1,6]+2*Km[t-1,7]+Km[t-1,8])/6
       alpha[t-1,3] <- theta[t-1,3]+(Km[t-1,9]+2*Km[t-1,10]+2*Km[t-1,11]+Km[t-1,12])/6
      
        # theta 服从 Dirichlet 分布
       theta[t,1:3] ~ ddirch(k*alpha[t-1,1:3])
        # Y, R 分别表示 在时间 t 时, I 状态和 R 状态人数所占比例
       Y[t-1] ~ dbeta(lambdaY*theta[t,2],lambdaY*(1-theta[t,2]))
       R[t-1] ~ dbeta(lambdaR*theta[t,3],lambdaR*(1-theta[t,3]))
     }

    # 定义方程的初值条件
    theta0[1:3]<-c(", 1 - Y[1] - R[1], ",", Y[1], ",", R[1], ")
    theta[1,1:3] ~ ddirch(theta0[1:3])
    gamma ~  dlnorm(", lognorm_gamma_parm$mu, ",", 1 / lognorm_gamma_parm$var, ")
    R0 ~ dlnorm(", lognorm_R0_parm$mu, ",", 1 / lognorm_R0_parm$var, ")
    beta <- R0*gamma
    k ~  dgamma(2,0.0001)
    lambdaY ~ dgamma(2,0.0001)
    lambdaR ~ dgamma(2,0.0001)
  }
")
  
 # 将上述方程(字符串)转换为rjags可识别的语言
  model.spec <- textConnection(model1.string)
  
  Sys.setenv(JAGS_HOME='D:/tools/JAGS-4.3.0')
  library(rjags)
  
# 定义采样函数(模型)
  posterior <- suppressWarnings(
    jags.model(
      model.spec,
      data = list(
        "Y" = Y,
        "R" = R,
        "T_prime" = T_prime,
        "pi" = pi
      ),
      n.chains = nchain,
      n.adapt = nadapt
    )
  )
  
  suppressWarnings(
    update(posterior, nburnin)
  ) # burn-in
  
# 对模型进行MCMC采样
  jags_sample <- suppressWarnings(
    jags.samples(
     # 采样函数
      posterior,
     # 需要采样的参数
      c("theta",
        "gamma",
        "R0",
        "beta",
        "Y",
        "lambdaY",
        "lambdaR",
        "k"),
      # 迭代次数
      n.iter = M,
      thin = thn
    )
  )

其中：
1.Km[t-1,1]代表 kS1t， Km[t-1,1] <- -beta*pi[t-1]*theta[t-1,1]*theta[t-1,2] 代表


2.Km[t-1,2]代表 kS2t，Km[t-1,2] <- -beta*pi[t-1]*(theta[t-1,1]+0.5*Km[t-1,1])*(theta[t-1,2]+0.5*Km[t-1,5]) 代表
  1. Km[t-1,3]代表 kS3t，Km[t-1,3] <- -beta*pi[t-1]*(theta[t-1,1]+0.5*Km[t-1,2])*(theta[t-1,2]+0.5*Km[t-1,6]) 代表
  2. Km[t-1,4]代表 kS4t，Km[t-1,4] <- -beta*pi[t-1]*(theta[t-1,1]+Km[t-1,3])*(theta[t-1,2]+Km[t-1,7]) 代表
  3. Km[t-1,5]代表 kI1t，Km[t-1,5] <- -Km[t-1,1]-Km[t-1,9] 代表
  4. Km[t-1,6]代表 kI2t，Km[t-1,6] <- -Km[t-1,2]-Km[t-1,10] 代表
  5. Km[t-1,7]代表 kI3t，Km[t-1,7] <- -Km[t-1,3]-Km[t-1,11] 代表
  6. Km[t-1,8]代表 kI4t，Km[t-1,8] <- -Km[t-1,4]-Km[t-1,12] 代表
  7. Km[t-1,9]代表 kR1t，Km[t-1,9] <- gamma*theta[t-1,2] 代表
  8. Km[t-1,10]代表 kR2t，Km[t-1,10] <- gamma*(theta[t-1,2]+0.5*Km[t-1,5]) 代表
  9. Km[t-1,11]代表 kR3t，Km[t-1,11] <- gamma*(theta[t-1,2]+0.5*Km[t-1,6]) 代表
  10. Km[t-1,12]代表 kR4t，Km[t-1,12] <- gamma*(theta[t-1,2]+Km[t-1,7]) 代表
  11. alpha[t-1,1]代表 θSt，alpha[t-1,1] <- theta[t-1,1]+(Km[t-1,1]+2*Km[t-1,2]+2*Km[t-1,3]+Km[t-1,4])/6 代表
  12. alpha[t-1,2]代表 θIt，alpha[t-1,2] <- theta[t-1,2]+(Km[t-1,5]+2*Km[t-1,6]+2*Km[t-1,7]+Km[t-1,8])/6 代表
  13. alpha[t-1,3]代表 θRt，alpha[t-1,3] <- theta[t-1,3]+(Km[t-1,9]+2*Km[t-1,10]+2*Km[t-1,11]+Km[t-1,12])/6 代表
    参数的采样：
  14. theta[t,1:3] ~ ddirch(k*alpha[t-1,1:3]) 代表 参数 θ 的采样分布，alpha指代内容见 13，14，15
  15. Y 表示在时间 t 时, I 状态人数所占比例，Y[t-1] ~ dbeta(lambdaY*theta[t,2],lambdaY*(1-theta[t,2]))，服从beta分布
  16. R 表示在时间 t 时, R 状态人数所占比例，R[t-1] ~ dbeta(lambdaR*theta[t,3],lambdaR*(1-theta[t,3]))，服从beta分布
    参数的初值条件：
  17. theta0[1:3]<-c(1 - Y[1] - R[1], Y[1], R[1]); theta[1,1:3] ~ ddirch(theta0[1:3]) 代表 表征 θ 的初值条件 θ0
  18. gamma ~ dlnorm(lognorm_gamma_parm$mu, 1 / lognorm_gamma_parm$var) 代表
  19. R0 ~ dlnorm(lognorm_R0_parm$mu, 1 / lognorm_R0_parm$var) 代表
  20. k ~ dgamma(2,0.0001)，lambdaY ~ dgamma(2,0.0001) 和 lambdaR ~ dgamma(2,0.0001) 分别代表 其中 lambdaY 代表 λI，lambdaR 代表 λR

微分方程模型获得MCMC估计的参数

由于估计参数所使用的时间点为30d，并且

posterior <- suppressWarnings(
    jags.model(
      model.spec,
      data = list(
        "Y" = Y,
        "R" = R,
        "T_prime" = T_prime,
        "pi" = pi
      ),
      n.chains = nchain,
      n.adapt = nadapt
    )
  )

jags_sample <- suppressWarnings(
    jags.samples(
     # 采样函数
      posterior,
     # 需要采样的参数
      c("theta",
        "gamma",
        "R0",
        "beta",
        "Y",
        "lambdaY",
        "lambdaR",
        "k"),
      # 迭代次数
      n.iter = M,
      thin = thn
    )
  )

上述参数 n.chains=4, M=500, thin=10 这里每个参数的MCMC相当于做了500次迭代，每10次迭代产生的数值取一次均值，一共产生50个数值。由于 n.chains=4，所以每个参数有200个数值

# 整理MCMC估计的参数
## 提取MCMC估计的参数 R0
R0_p <- unlist(as.mcmc.list(jags_sample$R0))
## 提取MCMC估计的参数 gamma
gamma_p <- unlist(as.mcmc.list(jags_sample$gamma))
## 提取MCMC估计的参数 beta
beta_p <- unlist(as.mcmc.list(jags_sample$beta))
## 提取MCMC估计的参数 lambdaY (λI)
lambdaY_p <- unlist(as.mcmc.list(jags_sample$lambdaY))
## 提取MCMC估计的参数 lambdaR (λR)
lambdaR_p <- unlist(as.mcmc.list(jags_sample$lambdaR))
## 提取MCMC估计的参数 k
k_p <- unlist(as.mcmc.list(jags_sample$k))

## 提取MCMC估计的参数 theta
### 将其整理为3个数组，每个数组行为200，列为31 (对应每一个时间点(步长)) 被估计的参数
### 这三个数组分别代表 θS, θI, θR
theta_p <- array(
  Reduce(
    rbind, as.mcmc.list(jags_sample$theta)
  ),
  dim = c(len, T_prime + 1, 3)
)

## 取最后一个时间点三个状态 (θS, θI, θR) 均值 
theta_p_mean <- apply(theta_p[, T_prime + 1, ], 2, mean)

## 查看最后一个时间点三个状态 (θS, θI, θR) 的分位数
theta_p_ci <- as.vector(
  apply(
    theta_p[, T_prime + 1, ],
    2,
    quantile,
    c(0.025, 0.5, 0.975)
  )
)

# 取R0的均值
R0_p_mean <- mean(R0_p)

## 查看 R0 的分位数
R0_p_ci <- quantile(
  R0_p,
  c(0.025, 0.5, 0.975)
)

# 取gamma的均值
gamma_p_mean <- mean(gamma_p)

## 查看 gamma 的分位数
gamma_p_ci <- quantile(
  gamma_p,
  c(0.025, 0.5, 0.975)
)

## 取 beta 的均值
beta_p_mean <- mean(beta_p)

## 查看 beta 的分位数
beta_p_ci <- quantile(
  beta_p,
  c(0.025, 0.5, 0.975)
)

## 取 lambdaY (λI) 的均值
lambdaY_p_mean <- mean(lambdaY_p)

## 查看 lambdaY (λI) 的分位数
lambdaY_p_ci <- quantile(
  lambdaY_p,
  c(0.025, 0.5, 0.975)
)

## 取 lambdaR (λR) 的均值
lambdaR_p_mean <- mean(lambdaR_p)

## 查看 lambdaR (λR) 的分位数
lambdaR_p_ci <- quantile(
  lambdaR_p,
  c(0.025, 0.5, 0.975)
)

## 取 k 的均值
k_p_mean <- mean(k_p)

## 查看 k 的分位数
k_p_ci <- quantile(
  k_p,
  c(0.025, 0.5, 0.975)
)

其中：

  1. theta_p 代表被拆分成的三个数组分别代表三个状态 θS，θI，θR，分数组的原理是基于每个时间点MCMC估计的数组数量（200）作为行数，一共31个时间点（第一个时间点为初值条件）作为列数，theta 估计三个状态（S，θI，θR）因此分为三个数组。并且这三个状态 θS，θI，θR 估计的数组仅限于30d内的，如果像往后续估计参见下一章节
  2. "gamma"，"R0"，"beta"，"Y"，"lambdaY"，"lambdaR"，"k" 参数估计的结果仅有一列向量（长度为200）

利用已估计的参数预测未来 θS，θI，θR 的走势

#### Forecast ####
library(gtools)

theta_pp <- array(0, dim = c(len, T_fin - T_prime, 3))
Y_pp <- matrix(NA, nrow = len, ncol = T_fin - T_prime)
R_pp <- matrix(NA, nrow = len, ncol = T_fin - T_prime)
for (l in 1:len) {
  ## 选择30d最后一个时间点作为后续时间段估计的起点
  thetalt1 <- theta_p[l, T_prime + 1, 1]
  thetalt2 <- theta_p[l, T_prime + 1, 2]
  thetalt3 <- theta_p[l, T_prime + 1, 3]

  betal <- c(beta_p)[l]
  gammal <- c(gamma_p)[l]
  kt <- c(k_p)[l]
  lambdaYl <- c(lambdaY_p)[l]
  lambdaRl <- c(lambdaR_p)[l]
  if (betal < 0 | gammal < 0 | thetalt1 < 0 |
      thetalt2 < 0 | thetalt3 < 0) next
  for (t in 1:(T_fin - T_prime)) {
    # T_fin = 200，T_prime = 30
    Km <- NULL
    alpha_pp <- NULL

   ## RK4 估计法
    Km[1] <- -betal * pi[t + T_prime] * thetalt1 * thetalt2
    Km[9] <- gammal * thetalt2
    Km[5] <- -Km[1] - Km[9]
    
    Km[2] <- -betal *
      pi[t + T_prime] * (thetalt1 + 0.5 * Km[1]) * (thetalt2 + 0.5 * Km[5])
    Km[10] <- gammal * (thetalt2 + 0.5 * Km[5])
    Km[6] <- -Km[2] - Km[10]
    
    Km[3] <- -betal *
      pi[t + T_prime] * (thetalt1 + 0.5 * Km[2]) * (thetalt2 + 0.5 * Km[6])
    Km[11] <- gammal * (thetalt2 + 0.5 * Km[6])
    Km[7] <- -Km[3] - Km[11]
    
    Km[4] <- - betal *
      pi[t + T_prime] * (thetalt1 + Km[3]) * (thetalt2 + Km[7])
    Km[12] <- gammal * (thetalt2 + Km[7])
    Km[8] <- -Km[4] - Km[12]
    
    alpha_pp[1] <- thetalt1 + (Km[1] + 2 * Km[2] + 2 * Km[3] + Km[4]) / 6
    alpha_pp[2] <- thetalt2 + (Km[5] + 2 * Km[6] + 2 * Km[7] + Km[8]) / 6
    alpha_pp[3] <- thetalt3 + (Km[9] + 2 * Km[10] + 2 * Km[11] + Km[12]) / 6
    
    # 服从dirichlet 分布
    thetalt_tmp <- rdirichlet(1, kt * c(alpha_pp))

    ## thetalt1 (theta_pp[l, t, 1]) 代表后 170 d 的 θS
    thetalt1 <- theta_pp[l, t, 1] <- thetalt_tmp[1]
    ##  ## thetalt1 (theta_pp[l, t, 1]) 代表后 170 d 的 θS
    thetalt2 <- theta_pp[l, t, 2] <- thetalt_tmp[2]
    ## thetalt1 (theta_pp[l, t, 1]) 代表后 170 d 的 θR
    thetalt3 <- theta_pp[l, t, 3] <- thetalt_tmp[3]
    
    Y_pp[l, t] <- rbeta(1, lambdaYl * thetalt2, lambdaYl * (1 - thetalt2))
    R_pp[l, t] <- rbeta(1, lambdaRl * thetalt3, lambdaRl * (1 - thetalt3))
  }
}

其中：

  1. thetalt1 <- theta_p[l, T_prime + 1, 1]; thetalt2 <- theta_p[l, T_prime + 1, 2]; thetalt3 <- theta_p[l, T_prime + 1, 3] 代表以第30天为后续估计的起点，估计后 170d θS，θI，θR 的数值，T_prime + 1 代表前 30d 最后一天
  2. thetalt1 (theta_pp[l, t, 1]) 代表后 170 d 的 θS，thetalt1 (theta_pp[l, t, 1]) 代表后 170 d 的 θS，thetalt1 (theta_pp[l, t, 1]) 代表后 170 d 的 θR

最终估计任然是利用RK4估计后170 d的数据

上一篇下一篇

猜你喜欢

热点阅读