随机微分方程以及常微分方程MCMC参数估计在流行病学中的应用

参考文献：

  1. 《Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study》
  2. 《Early dynamics of transmission and control of COVID-19: a mathematical modelling study》
  3. 《Capturing the time-varying drivers of an epidemic using stochastic dynamical systems》
  4. 《Reconstruction of the full transmission dynamics of COVID-19 in Wuhan》
  5. MH采样

记一个一维随机微分方程的例子

常规的微分方程的构成：


随机微分方程：



随机微分方程即在常规微分方程后面添加一个布朗运动随机项即可

随机微分方程运用在传染病模型中


上图式子 dXt 这项代表上面常规微分方程的每一项改变为随机微分方程的形式，随机微分方程相较于常规微分方程来说体现了人群接触的随机性

利用R包Sim.DiffProc可计算随机微分方程，其中一维的Euler法解微分方程采取的是定义 dt = Δt = 某一微小值 然后采取马氏链的随机游走解法

扩散过程

核心代码：

# 重复扩散过程
 for (i in 1:N) {
    X[i + 1,] <- X[i,] + A(t[i],X[i,]) * Dt + S(t[i],X[i,]) * W[i,]  
  }

其中，W为随机矩阵，模拟布朗运动的随机性，A(t[i],X[i,])为参数，S(t[i],X[i,])为随机过程参数，Dt 为微小时间，为一确定值；X[i,] 为第 i 个点；X[i + 1,]为第 i + 1 个点

常微分方程利用MCMC估计参数的例子

1. 理论部分

其中：

  1. Di：Di = 2.9，I/A 状态到 R（康复）状态的周期（I 和 A 状态的患者直接康复的时间，没有专门去医院治疗）
  2. Dp：Dp = 2.3，P 状态到 I/A 状态的周期
  3. De：De = 2.9，E 状态到 P 状态的周期（潜伏期）
  4. Dq：Dq = c(21, 15, 10, 6, 3)，Dq 值分为五个阶段，I 状态到 H 状态的周期（确诊COVID19的患者到医院治疗所花的时间）
  5. Dh：Dh = 30，H 状态到 R（康复）状态的周期（确诊COVID19的患者在医院治疗到康复的时间）
  6. alpha：代表传播率，0.55
  7. N = 10000000，代表某城市的常住人口数
  8. n，分为5各阶段，代表每一阶段的流动人口数
  9. S 代表易感人群
  10. E 代表被COVID19病毒入侵后的无症状人群
  11. P 代表由类似症状，但尚未确诊的人；分为两部分，一部分发展为确诊病例，另外一部分发展为有类似症状但不是COVID19的病例，大概有为 r 比例的 P 人群会确诊为COVID19患者
  12. A 为类似症状但不是COVID19的病例人群
  13. I 为确诊COVID19病例的人群，这一部分分为两步，一部分是自己没去医院痊愈了，治愈周期为 Di ；另一部分为去医院治疗的患者，所花费时间为 Dq（确诊COVID19的患者到医院治疗所花的时间）+ Dh（确诊COVID19的患者在医院治疗到康复的时间）
  14. H 为在医院接受治疗的确诊患者
  15. R 为康复人群

动力学系统为：


其中：

  1. α：代表未确诊患者转为确诊病例的比例
  2. b：代表传播率；bS × αP 代表 P 群体确诊部分的人群使得 S（易感人群）exposed 变为 E 状态；bS × αA 代表 A 群体确诊部分的人群使得 S（易感人群）exposed 变为 E 状态；bS × I 代表 I 群体（确诊人群）使得 S（易感人群）exposed 变为 E 状态
  3. nS/N 代表流出本市的人群
  4. n 代表流入本市的人群，均为易感者 S
  5. 形如 E/De 这样的代表在 De 时期内平均每天由 E 状态转移到 P 状态的人群数量，其他类似理解
  6. r 代表有 r 比例的 P 人群会确诊为COVID19患者

由于政府采取不同的措施，因此参数传播率 b，以及参数 r ， 即 r 比例的 P 人群会确诊为COVID19患者 (b12，b3，b4，b5，r12，r3，r4，r5；角标12，3，4，5分别代表第一，二阶段，第三阶段，第四阶段和第五阶段) 将分为5个不同的阶段；并且流动人口数 n 也分为五个阶段，即正常情况的人口流动，春运时的人口流动，封锁时的人口流动；I 状态到 H 状态的周期（确诊COVID19的患者到医院治疗所花的时间）Dq 也分为5个阶段，前几个阶段将确诊为COVID19的患者送达医院的时间较慢长。采取措施后，后几个阶段将确诊COVID19的患者送达医院的时间逐步地加缩短了，其送医院周期分别为 21, 15, 10, 6, 3天

2. 核心代码解析

数据源为：


  1. index 表示不同的日期
  2. 第一列表示当天确诊COVID19的人数
  3. 第二列表示累计确诊COVID19的人数

脚本：https://github.com/hzaurzli/SAPHIRE/blob/master/R/

2.1 初始化参数：

# 定义初始值
R0 <- 0
H0 <- 27
r0 = 0.23
Di = 2.9
Dp = 2.3
De = 2.9
Dq = c(21, 15, 10, 6, 3)
alpha = 0.55
Dh = 30
N = 10000000

realData_all <- read.csv("../data/Covid19CasesWH.csv", row.names = 1)  
realData <- realData_all[-c(1:24), ] # the 25th row correspond to 1 Jan
jan1_idx = 25

#cases from Jan 1 to Feb 29，选择分析的时间段
daily_new_case <- realData[1:60, 1]
daily_new_case_all <- realData[, 1]

# 初始化的时间窗口是 Jan 1 及以后的时间段
## Jan 3-5 for De=2.9 and Dp=2.3，定义周期为 Jan 3-5
## 定义初始 E 群体， 在 Jan 3-5 确诊的患者在 Jan 1 中都称为 r0 × E 群体（时间周期为 De + Dp），r0 代表初始的确诊率
## E0 <- (40 + 23 + 47) / r0，E 群体中有 r0 比例左右的人确诊为 COVID19，E为被病毒入侵的群体
E0 <- sum(realData_all[(jan1_idx+round(Dp)):(jan1_idx+round(Dp)+round(De)-1),1]) / r0 

## Jan 1-2 for Dp=2.3，定义初始 P 群体， 从Jan 1-2 确诊的患者在 Jan 1 中都称为 r0 × P 群体（时间周期为 Dp）
## P0 <- (41 + 34) / r0，P 群体中有 r0 比例左右的人确诊为 COVID19，P为有症状但没有确诊                       
P0 <- sum(realData_all[jan1_idx:(jan1_idx+round(Dp)-1),1]) / r0                     

# 定义 I 群体，Di 周期前（Dec 29-31）确诊的患者为 Jan 1 的 I 群体，因为在 Di 周期内尚未康复，因此到 Jan 1 仍然属于确诊患者
## Dec 29-31 for Di=2.9，I0 <- 11 + 13 + 10 
I0 <- sum(realData_all[(jan1_idx-round(Di)):(jan1_idx-1),1])

## 有类似症状但不是COVID19的群体 A，r0 × P = I0，(1 - r0) × P = A0，联立的得 A0 = I0 × ( 1 - r0) / r0                            
A0 <- I0 * (1 - r0) / r0

# 定义 易感人群 E
S0 <- N - E0 - P0 - I0 - A0 - H0 - R0
init_states <- round(c(S = S0, E = E0, P = P0, I = I0, A = A0, H = H0, R = R0), 0)

其中：
  1. N 为总人口，初始的康复人群 R0 = 0，研究的起始时间点为 Jan 1
  2. 初始化的时间窗口是 Jan 1 及以后的时间段，Jan 3-5 for De=2.9 and Dp=2.3，定义周期为 Jan 3-5，定义初始 E 群体， 在 Jan 3-5 确诊的患者在 Jan 1 中都称为 r0 × E 群体（时间周期为 De + Dp），r0 代表初始的确诊率
  3. Jan 1-2 for Dp=2.3，定义初始 P 群体， 从Jan 1-2 确诊的患者在 Jan 1 中都称为 r0 × P 群体（时间周期为 Dp），P0 <- (41 + 34) / r0，P 群体中有 r0 比例左右的人确诊为 COVID19，P为有症状但没有确诊
  4. 定义 I 群体，Di 周期前（Dec 29-31）确诊的患者为 Jan 1 的 I 群体，因为在 Di 周期内尚未康复，因此到 Jan 1 仍然属于确诊患者，Dec 29-31 for Di=2.9，I0 <- 11 + 13 + 10
  5. 有类似症状但不是COVID19的群体 A，r0 × P = I0，(1 - r0) × P = A0，联立的得 A0 = I0 × ( 1 - r0) / r0
  6. 定义 易感人群 E 为 N - E0 - P0 - I0 - A0 - H0 - R0

  7. 结果：

对初始化部分的返回结果进行讨论：

return(list(Di=Di,
              Dp=Dp,
              De=De,
              Dq=Dq,
              alpha=alpha,
              Dh=Dh,
              N=N,
              flowN=flowN,
              daily_new_case = daily_new_case, 
              daily_new_case_all = daily_new_case_all, 
              init_states = init_states,
              days_to_fit=1:60,
              stage_intervals=list(
                c(start=1, end=9),
                c(start=10, end=22),
                c(start=23, end=32),
                c(start=33, end=47),
                c(start=48, end=60)
              ),
              var_trans_fun=transform_var_main_5stage,
         par_lower = c(b12 = 0, b3 = 0, b4 = 0, b5 = 0, r12 = 0, delta3 = -10, delta4 = -10, delta5 = -10),
         par_upper = c(b12 = 2, b3 = 2, b4 = 2, b5 = 2, r12 = 1, delta3 = 10, delta4 = 10, delta5 = 10)))


## 待估计的参数转换公式
transform_var_main_5stage=function(pars) {
    b_vec <- pars[1:4]
    b_vec <- c(b_vec[1], b_vec[1], b_vec[2:4])
    ##
    r12 <- pars[5]
    r3 <- 1 / (1 + (1 - r12) / (r12 * exp(pars[6])))
    r4 <- 1 / (1 + (1 - r3) / (r3 * exp(pars[7])))
    r5 <- 1 / (1 + (1 - r4) / (r4 * exp(pars[8])))
    r_vec <- c(r12,r12,r3,r4,r5)
    
    return(list(b_vec, r_vec))
  }

其中：作者研究的时间段为 Jan 1 - Feb 29 为 60 d

  1. Di：Di = 2.9，I/A 状态到 R（康复）状态的周期（I 和 A 状态的患者直接康复的时间，没有专门去医院治疗）
  2. Dp：Dp = 2.3，P 状态到 I/A 状态的周期
  3. De：De = 2.9，E 状态到 P 状态的周期（潜伏期）
  4. Dq：Dq = c(21, 15, 10, 6, 3)，Dq 值分为五个阶段，I 状态到 H 状态的周期（确诊COVID19的患者到医院治疗所花的时间）
  5. Dh：Dh = 30，H 状态到 R（康复）状态的周期（确诊COVID19的患者在医院治疗到康复的时间）
  6. alpha：代表传播率，0.55
  7. N = 10000000，代表某城市的常住人口数
  8. flowN，flowN = c(500000, 800000, 0, 0, 0)，分为5各阶段，代表每一阶段的流动人口数
  9. daily_new_case 和 daily_new_case_all 分别代表每日新增的确诊COVID19人数和每日累计的COVID19人数
  10. init_states ：代表初始各个状态人群的初始值，init_states <- round(c(S = S0, E = E0, P = P0, I = I0, A = A0, H = H0, R = R0), 0)
  11. days_to_fit=1:60 代表研究的时间一个60 d
  12. stage_intervals，代表 stage_intervals=list(c(start=1, end=9), c(start=10, end=22),c(start=23, end=32),c(start=33, end=47),c(start=48, end=60))，代表 60 d 内五个不同阶段对应的时间
  13. par_lower 和 par_upper 代表参数的上下限
  14. transform_var_main_5stage()，这里有5个阶段的参数，（r12，b12） 代表第一，第二阶段的；（r3，b3）代表第三阶段，（r4，b4）代表第四阶段；（r5，b5）代表第五阶段；（r12，r3，r4，r5）代表5个不同阶段 P 群体（有类似症状但尚未确诊的人群）中确诊COVID19患者（I 群体）的比例；（b12，b3，b4，b5）代表5个不同阶段病毒的传播率
    待估计的参数转换公式

2.2 定义先验分布，采样函数和似然函数

首先是先验分布和采样函数

# default_pars_density() 为先验分布的似然函数，pars 为每次采样的待估计参数值
## 例如 d_vec[6] <- dnorm(delta_3, delta_mean, delta_sd, log = T)，代表在正态分布中(均值为delta_mean，方差为 delta_sd) ，自变量为 delta_3 时所对应的似然值

# 定义先验分布的均值和标准差
delta_mean <- 0
delta_sd <- 1
beta_shape1 <- 7.3
beta_shape2 <- 24.6

default_pars_density <- function(pars) {
  d_vec <- rep(NA, 8)
  ##b12, b3, b4, b5
  for(i in c(1:4)) {
    d_vec[i] <- dunif(pars[i], 0, 2, log = T)
  }
  ## r12
  r12 <- pars[5]
  d_vec[5] <- dbeta(r12, beta_shape1, beta_shape2, log = T)
  ## r3
  delta_3 <- pars[6]
  d_vec[6] <- dnorm(delta_3, delta_mean, delta_sd, log = T)
  ## r4
  delta_4 <- pars[7]
  d_vec[7] <- dnorm(delta_4, delta_mean, delta_sd, log = T)
  ## r5
  delta_5 <- pars[8]
  d_vec[8] <- dnorm(delta_5, delta_mean, delta_sd, log = T)
  ##
  return(sum(d_vec))
}

# 定义采样函数，MCMC 每迭代一次，分别在待估计参数对应的分布中采样一次
## 采样的分布与上面先验分布函数是一样的
default_pars_sampler <- function(n = 1) {
  s_vec <- matrix(NA, n, 8)
  ## b12, b3, b4, b5
  for(i in c(1:4)) {
    s_vec[, i] <- runif(n, 0, 2) 
  }
  ## r12 
  r12 <- rbeta(n, beta_shape1, beta_shape2)
  s_vec[, 5] <- r12
  ## r3
  s_vec[, 6] <- rnorm(n, delta_mean, delta_sd)
  ## r4
  s_vec[, 7] <- rnorm(n, delta_mean, delta_sd)
  ## r5
  s_vec[, 8] <- rnorm(n, delta_mean, delta_sd)
  return(s_vec)
}

其中：

  1. default_pars_density() 为先验分布的似然函数，pars 为每次采样的待估计参数值。例如 d_vec[6] <- dnorm(delta_3, delta_mean, delta_sd, log = T)，代表在正态分布中(均值为delta_mean，方差为 delta_sd) ，自变量为 delta_3 时所对应的似然值（dnorm代表正态分布的似然函数）
  2. 定义采样函数，MCMC 每迭代一次，分别在待估计参数对应的分布中采样一次。采样的分布与上面先验分布函数是一样的（ rnorm(n=1) 代表在正态分布中随机采样 n=1 次）
  3. 这里的迭代的次数为 180000 次，因此要采 180000 次样
  4. 定义先验分布函数的返回值时，返回值需要写成各个参数似然值之和的形式
  5. 定义采样函数的返回值时，可以用一个向量来储存各个参数的一次采样的值

其次是似然函数
作者原文对似然函数的定义是：


即待估计参数 b12，b3，b4，b5，r12，r3，r4，r5 的似然值服从上图所示的泊松分布

# 定义似然函数
loglh_func <- function(pars){
    ## 估计值为 ypred，真实值为 onset_obs (60 d，Jan 1 - Feb 29)
    ## 函数 SEIRpred() 将在 2.3-2.4 章节进行叙述
    ypred <- SEIRpred(pars, init_settings = init_sets_list)
    ypred <- ypred[, "Onset_expect"]
    
    # meant to suppress warnings when ypred is negative
    suppressWarnings(p <- dpois(onset_obs, ypred))
    
    if(any(p == 0) || any(is.nan(p))){
      logL <- -Inf
    }else{
      logL <- sum(log10(p))
    }
    return(logL)
  }

其中，作者在这里仅仅预测 60 d 的数据（Jan 1 - Feb 29）

  1. 泊松分布的参数与各个状态的人群有关，而各个状态的人群数量与待估计的参数 b12，b3，b4，b5，r12，r3，r4，r5 有关，因此各个状态的人群服从泊松分布即为 b12，b3，b4，b5，r12，r3，r4，r5 服从泊松分布
  2. 每进行一次采样，更新一次参数则重新计算一次 SEIRpred() 函数

2.3 先验分布，采样函数和似然函数之间的关系

1. 先验分布： 所谓先验分布就是事先定义的待估计参数的分布

每一次对参数的采样都是基于先验分布进行随机采样

2. 采样函数： 顾名思义采样函数即为MCMC每一次采样所对应的函数，采样函数的分布就是先验分布，MCMC每一次采样都是从先验分布中随机采样

3. 似然函数： MCMC的采样是接受-拒绝采样，因此每次采样获得的参数值带入似然函数计算出来的似然值，需要与均匀分布的似然值（均匀分布的似然值在每次采样中也是随机获得的）作比较，再决定是否接受该参数的采样。该例子中的似然值定义为以 Jan 1 - Feb 29 期间每日真实的新增的确诊人数为泊松分布的期望。每采样一次所获得的参数，带入动力学方程计算出来的估计值（pred），再带入泊松分布中可以计算 Jan 1 - Feb 29 每一天的对数似然值，求和后作为总的对数似然值再与均匀分布的似然值（均匀分布的似然值在每次采样中也是随机获得的）作比较，判断是否接受本次采样

 ## 估计值为 ypred，真实值为 onset_obs (60 d，Jan 1 - Feb 29)
 ## 函数 SEIRpred() 将在 2.3-2.4 章节进行叙述
ypred <- SEIRpred(pars, init_settings = init_sets_list)
ypred <- ypred[, "Onset_expect"]
 
p <- dpois(onset_obs, ypred)

## 似然值求对数和
if(any(p == 0) || any(is.nan(p))){
     logL <- -Inf
}else{
     logL <- sum(log10(p))
}

2.4 解析SEIRpred()函数

# 其中一组参数作为例子
pars = c(1.4, 0.4, 0.1, 0.1, 0.5, -1, 0, 0)
init_settings = init_sets_list

SEIRpred <- function(pars, 
                     init_settings) {

  tmp_ret=init_settings$var_trans_fun(pars)
  b_vec=tmp_ret[[1]]
  r_vec=tmp_ret[[2]]
  stage_intervals=init_settings$stage_intervals
  n_stage=length(stage_intervals)
  
  ## 读取初始值
  Di <- init_settings$Di
  Dp <- init_settings$Dp
  De <- init_settings$De
  Dq_vec <- init_settings$Dq
  alpha <- init_settings$alpha
  Dh <- init_settings$Dh
  N <- init_settings$N
  flowN_vec <- init_settings$flowN
  init_states <- init_settings$init_states
  days_to_fit <- init_settings$days_to_fit
 
  ## matrix for results
  states_mat <- matrix(0, length(days_to_fit), length(init_states) + 2)
  states_mat[, 1] <- days_to_fit
  colnames(states_mat) <- c("time", "S", "E", "P", "I", "A", "H", "R", "Onset_expect")
  
  myold_states <- init_states
  
## 分 5 个阶段进行计算，因为 5 个阶段的参数是不一样的
  for (i_stage in 1:n_stage) {
    stage_pars_setings <- c(b = b_vec[i_stage], r = r_vec[i_stage], Dq = Dq_vec[i_stage], n = flowN_vec[i_stage])
    ## 对于每一个阶段计算每一天各个状态的人数
    for (d in stage_intervals[[i_stage]][["start"]]:stage_intervals[[i_stage]][["end"]]) {
      states_mat[d, -1] <- update_func(stage_pars = stage_pars_setings, states_old = myold_states)
      myold_states <- states_mat[d, -1]
    }
  }
  return(states_mat)
}


## ODE function based on deterministic SEIR model
## 更迭函数
update_func <- function(stage_pars, states_old) {
  ## stage pars
  b <- stage_pars[1]
  r <- stage_pars[2]
  Dq <- stage_pars[3]
  n <- stage_pars[4]
  ## old states number: c(S, E, P, I, A, H, R)
  ## 读取上一个时间点各个状态的人数数据
  ## c(S, E, P, I, A, H, R)  为上一个时间点各个状态的人数数据
  S <- states_old[1]
  E <- states_old[2]
  P <- states_old[3]
  I <- states_old[4]
  A <- states_old[5]
  H <- states_old[6]
  R <- states_old[7]
  ## new values
  ## 计算新时间点 (本时间点) 各个状态的人数数据
  ## c(S, E, P, I, A, H, R)  为上一个时间点各个状态的人数数据
  ## c(S_new, E_new, P_new, I_new, A_new, H_new, R_new) 为新时间点 (本时间点) 各个状态的人数数据
  S_new <- S - b * S * (alpha * P + I + alpha * A) / N + n - n * S / N
  E_new <- E + b * S * (alpha * P + I + alpha * A) / N - E / De - n * E / N
  P_new <- P +  E / De  - P / Dp - n * P / N
  I_new <- I + r * P / Dp - I / Di - I / Dq
  A_new <- A + (1 - r) * P / Dp - A / Di - n * A / N
  H_new <- H + I / Dq - H / Dh
  R_new <- R + H / Dh + (A + I) / Di - n * R / N
  Onset_expect <- r * P / Dp
  ##
  return(c(S_new, E_new, P_new, I_new, A_new, H_new, R_new, Onset_expect))
}

其中：ODE 迭代的方程为
  1. 这里迭代方程新时间点新时间点 (本时间点) 各个状态的人数数据由上一个时间点各个状态的人数数据通过对应公式进行计算，是因为这里将微分方程的 Δt = 1 d；例如方程第一项 S_new<-S-b*S*(alpha*P+I+alpha*A)/N+n-n*S/N 等于 S_new<-(S-b*S*(alpha*P+I+alpha*A)/N+n-n*S/N)*Δt，此时微分方程的 Δt = 1 d
  2. 利用 MCMC 对参数 b12，b3，b4，b5，r12，r3，r4，r5 进行采样180000次；每采样一次得到的新一轮参数带入似然函数 loglh_func 中计算似然值，利用极大似然估计的思想，18000次采样中若其中某一次采样所得到的一组参数使得似然值最大，那么这一组参数即为最佳的期望参数组合，
loglh_func <- function(pars){
   ## 估计值为 ypred，真实值为 onset_obs (60 d，Jan 1 - Feb 29) I 状态的人数
   ## 函数 SEIRpred() 将在 2.3-2.4 章节进行叙述
   ## 180000 次参数采样，每采样一次得到一组参数，则带入SEIRpred进行 Jan 1 - Feb 29 期间每一个时间点 I 状态 的人数计算
   ypred <- SEIRpred(pars, init_settings = init_sets_list)
   ypred <- ypred[, "Onset_expect"]
   
   # meant to suppress warnings when ypred is negative
   ## 判断标准为带入每一组参数计算 Jan 1 - Feb 29 期间每一个时间点 I 状态的人数，然后和这 60 d 的真实值做比较，dpois(onset_obs, ypred)最大的参数组合为期望的参数组合
   suppressWarnings(p <- dpois(onset_obs, ypred))
   
   if(any(p == 0) || any(is.nan(p))){
     logL <- -Inf
   }else{
     logL <- sum(log10(p))
   }
   return(logL)
 }
  1. 180000 次参数采样，每采样一次得到一组参数，则带入 SEIRpred 进行 Jan 1 - Feb 29 进行每一个时间点 I 状态的人数计算；判断标准为带入每一组参数计算 Jan 1 - Feb 29 期间每一个时间点 I 状态 的人数，然后和这 60 d 的真实值（I 状态） 做比较，dpois(onset_obs, ypred) 似然值最大的参数组合为期望的参数组合；注意这里的采样并非每一个时间点对参数 b12，b3，b4，b5，r12，r3，r4，r5 采样180000次；而是先采样180000次得到180000个参数组合，然后分别带入函数 SEIRpred 进行计算，得到 I 状态 的dpois(onset_obs, ypred) 似然值最大时参数组合，则定义为期望的参数组合
    该图点坐标即为期望的参数组合所计算的每一个时间点 I 状态 的人群数量，阴影部分为其他参数组合所计算的置信区间

2.5 参数MCMC的过程

# 初始化
init_sets_list = generate_init_condi(r0 = 0.23)
delta_mean <- 0
delta_sd <- 1
beta_shape1 <- 7.3
beta_shape2 <- 24.6
init_sets_list = get_init_sets_list 
startValue=NA
panel_B_R_ylim=4
plot_combined_fig=T
pars_density=default_pars_density
pars_sampler=default_pars_sampler
pars_name=c("b12", "b3", "b4", "b5", "r12", "delta3", "delta4", "delta5")
calc_clearance=T
n_burn_in=4000
n_iterations=180000

randomize_startValue = T
run_id = "main_analysis" 
output_ret = T
skip_MCMC=F

# 得到初始的数据
## onset_obs 60 d 真实的新增患者数据 I 状态
onset_obs <- init_sets_list$daily_new_case
init_states <- init_sets_list$init_states
n_pars = length(pars_name)
n_stage = length(init_sets_list$stage_intervals)


# 定义先验分布和采样函数
pars_prior <- createPrior(density = pars_density, sampler = pars_sampler, 
                          lower = init_sets_list$par_lower, upper = init_sets_list$par_upper)

# 定义似然函数
bayesSEIR <- createBayesianSetup(loglh_func, prior = pars_prior)
startValue=pars_sampler()
while (is.infinite(loglh_func(startValue))) {
  startValue=pars_sampler()
}

# 初始化 MCMC
mh_settings = list(startValue = startValue,
                   adapt = T, DRlevels = 2, iterations = 180000, thin = 10)

# runMCMC
mh_out <- runMCMC(bayesianSetup = bayesSEIR, sampler = "Metropolis", settings = mh_settings)
# 得到每次采样的参数
mcmc_pars_estimate <- getSample(mh_out, start = n_burn_in+2, thin = 1)  ## set start = 2002 as the burn in period
mcmc_pars_estimate <- round(mcmc_pars_estimate, 3)
colnames(mcmc_pars_estimate) <- pars_name 

# 分阶段获取最佳（期望）参数组合值
## n_stage 代表 5 个不同的阶段
 if (n_stage>1) {
    for (i_stage in 1:n_stage) {
     ## 最佳参数值为各个阶段每个参数多次采样的均值
      r_str[[i_stage]]=paste0(round(mean(estRt_mat[i_stage,]),2), " (",
                              round(quantile(estRt_mat[i_stage,],0.025),2)," - " , 
                              round(quantile(estRt_mat[i_stage,],0.975),2), ")")
    }
  } else {
    ## 最佳参数值为各个阶段每个参数多次采样的均值，这里计算的是第一阶段每个参数多次采样的均值
    r_str[[1]]=paste0(round(mean(estRt_mat),2), " (",
                            round(quantile(estRt_mat,0.025),2)," - " , 
                            round(quantile(estRt_mat,0.975),2), ")")
  }

得到最佳（期望）参数组合值和每次采样的参数组合后即可轻易计算每一个时间点各个状态的人数了，一般情况下，180000次采样的最佳（期望组合值）组合值可以分阶段取各个参数的均值得到

因为作者采用的是MH接受-拒绝采样，因此在所接受的参数群体中，每一阶段对应的每个参数应该是服从正态分布的，因此期望即为每一阶段各个参数的均值

当然，将SEIRpred()函数换为deSolve包所计算的数值也是可以的，只需要确定步长后（假设步长是1）将每一个点对应的确诊COVID19人数获取到，然后和每天确诊COVID19真实值比较p <- dpois(onset_obs, ypred)，然后进行似然值比较，决定是否进行本次采样

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 204,732评论 6赞 478
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 87,496评论 2赞 381
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 151,264评论 0赞 338
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,807评论 1赞 277
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,806评论 5赞 368
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,675评论 1赞 281
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 38,029评论 3赞 399
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,683评论 0赞 258
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 41,704评论 1赞 299
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,666评论 2赞 321
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,773评论 1赞 332
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,413评论 4赞 321
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 39,016评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,978评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,204评论 1赞 260
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 45,083评论 2赞 350
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,503评论 2赞 343