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

作者: 小潤澤 | 来源:发表于2023-01-20 22:24 被阅读0次

随机微分方程以及常微分方程MCMC参数估计在流行病学中的应用
【常微分方程】$1.1 微分方程及其解的定义
高等数学(七)常微分方程
10-Python 科学计算_scipy 篇
常/偏微分方程求解
二阶变系数线性常微分方程的求解
你学高数了吗？我是焦头烂额！Python居然能解决高数问题？
前向欧拉方程
常微分方程
常微分方程

参考文献：

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

常规的微分方程的构成：

随机微分方程：

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

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

上图式子 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. 理论部分

其中：

Di：Di = 2.9，I/A 状态到 R（康复）状态的周期（I 和 A 状态的患者直接康复的时间，没有专门去医院治疗）

Dp：Dp = 2.3，P 状态到 I/A 状态的周期

De：De = 2.9，E 状态到 P 状态的周期（潜伏期）

Dq：Dq = c(21, 15, 10, 6, 3)，Dq 值分为五个阶段，I 状态到 H 状态的周期（确诊COVID19的患者到医院治疗所花的时间）

Dh：Dh = 30，H 状态到 R（康复）状态的周期（确诊COVID19的患者在医院治疗到康复的时间）

alpha：代表传播率，0.55

N = 10000000，代表某城市的常住人口数

n，分为5各阶段，代表每一阶段的流动人口数

S 代表易感人群

E 代表被COVID19病毒入侵后的无症状人群

P 代表由类似症状，但尚未确诊的人；分为两部分，一部分发展为确诊病例，另外一部分发展为有类似症状但不是COVID19的病例，大概有为 r 比例的 P 人群会确诊为COVID19患者

A 为类似症状但不是COVID19的病例人群

I 为确诊COVID19病例的人群，这一部分分为两步，一部分是自己没去医院痊愈了，治愈周期为 Di ；另一部分为去医院治疗的患者，所花费时间为 Dq（确诊COVID19的患者到医院治疗所花的时间）+ Dh（确诊COVID19的患者在医院治疗到康复的时间）

H 为在医院接受治疗的确诊患者

R 为康复人群

动力学系统为：

其中：

α：代表未确诊患者转为确诊病例的比例

b：代表传播率；bS × αP 代表 P 群体确诊部分的人群使得 S（易感人群）exposed 变为 E 状态；bS × αA 代表 A 群体确诊部分的人群使得 S（易感人群）exposed 变为 E 状态；bS × I 代表 I 群体（确诊人群）使得 S（易感人群）exposed 变为 E 状态

nS/N 代表流出本市的人群

n 代表流入本市的人群，均为易感者 S

形如 E/De 这样的代表在 De 时期内平均每天由 E 状态转移到 P 状态的人群数量，其他类似理解

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. 核心代码解析

数据源为：

index 表示不同的日期

第一列表示当天确诊COVID19的人数

第二列表示累计确诊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)

其中：

N 为总人口，初始的康复人群 R0 = 0，研究的起始时间点为 Jan 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 代表初始的确诊率

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

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

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

定义易感人群 E 为 N - E0 - P0 - I0 - A0 - H0 - R0

结果：

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

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

Di：Di = 2.9，I/A 状态到 R（康复）状态的周期（I 和 A 状态的患者直接康复的时间，没有专门去医院治疗）

Dp：Dp = 2.3，P 状态到 I/A 状态的周期

De：De = 2.9，E 状态到 P 状态的周期（潜伏期）

Dq：Dq = c(21, 15, 10, 6, 3)，Dq 值分为五个阶段，I 状态到 H 状态的周期（确诊COVID19的患者到医院治疗所花的时间）

Dh：Dh = 30，H 状态到 R（康复）状态的周期（确诊COVID19的患者在医院治疗到康复的时间）

alpha：代表传播率，0.55

N = 10000000，代表某城市的常住人口数

flowN，flowN = c(500000, 800000, 0, 0, 0)，分为5各阶段，代表每一阶段的流动人口数

daily_new_case 和 daily_new_case_all 分别代表每日新增的确诊COVID19人数和每日累计的COVID19人数

init_states ：代表初始各个状态人群的初始值，init_states <- round(c(S = S0, E = E0, P = P0, I = I0, A = A0, H = H0, R = R0), 0)

days_to_fit=1:60 代表研究的时间一个60 d

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 内五个不同阶段对应的时间

par_lower 和 par_upper 代表参数的上下限

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)
}

其中：

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

定义采样函数，MCMC 每迭代一次，分别在待估计参数对应的分布中采样一次。采样的分布与上面先验分布函数是一样的（ rnorm(n=1) 代表在正态分布中随机采样 n=1 次）

这里的迭代的次数为 180000 次，因此要采 180000 次样

定义先验分布函数的返回值时，返回值需要写成各个参数似然值之和的形式

定义采样函数的返回值时，可以用一个向量来储存各个参数的一次采样的值

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

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

# 定义似然函数
loglh_func <- function(pars){
    ## 估计值为 ypred，真实值为 onset_obs (60 d，Jan 1 - Feb 29)
    ## 函数 SEIRpred() 将在 2.3 章节进行叙述
    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）

泊松分布的参数与各个状态的人群有关，而各个状态的人群数量与待估计的参数 b12，b3，b4，b5，r12，r3，r4，r5 有关，因此各个状态的人群服从泊松分布即为 b12，b3，b4，b5，r12，r3，r4，r5 服从泊松分布

每进行一次采样，更新一次参数则重新计算一次 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 章节进行叙述
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 迭代的方程为

这里迭代方程新时间点新时间点 (本时间点) 各个状态的人数数据由上一个时间点各个状态的人数数据通过对应公式进行计算，是因为这里将微分方程的 Δ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

利用 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 章节进行叙述
   ## 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)
 }
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)，然后进行似然值比较，决定是否进行本次采样

随机微分方程以及常微分方程MCMC参数估计在流行病学中的应用
参考文献：《Nowcasting and forecasting the potential domestic ...
【常微分方程】$1.1 微分方程及其解的定义
1.常微分方程、偏微分方程概念常微分方程：未知函数是一元函数。偏微分方程：未知函数是多元函数，微分方程...
高等数学(七)常微分方程
（一）常微分方程的基本概念微分方程微分方程的阶/解微分方程的通解/特解初始条件积分曲线（二）一阶微分方...
10-Python 科学计算_scipy 篇
课程概要：1、非线性方程组2、数值积分3、常微分方程组 1、非线性方程组 2、数值积分（定积分） 3、常微分方程组
常/偏微分方程求解
常微分方程 1、四阶经典步长Runge-Kutta法求解高阶常微分方程 2、用ode15s函数、ode23s函数和...
二阶变系数线性常微分方程的求解
之前总结过一次常见的常微分方程的解法，不过里面并没有二阶变系数线性常微分方程的求解方法，现在补充如下。对于形式的...
你学高数了吗？我是焦头烂额！Python居然能解决高数问题？
本文将学习Python的Sympy库，通过该库来求极限，微分，定积分与不定积分，常微分方程、偏微分方程。求解方程...
前向欧拉方程
欧拉方法是一种一阶数值方法，用以对给定初值的常微分方程（即初值问题）求解。它是一种解决数值常微分方程的最基本的一类...
常微分方程
import numpy as np import matplotlib.pyplot as plt from s...
常微分方程
一、一阶线性微分方程一般形式：当时，则方程为齐次方程，反之为非齐次。通解：通解公式：特解：特解是将通解中的根...

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

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

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

1. 理论部分

2. 核心代码解析

2.1 初始化参数：

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

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

2.4 解析SEIRpred()函数

2.5 参数MCMC的过程

相关文章

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

【常微分方程】$1.1 微分方程及其解的定义

高等数学(七)常微分方程

10-Python 科学计算_scipy 篇

常/偏微分方程求解

二阶变系数线性常微分方程的求解

你学高数了吗？我是焦头烂额！Python居然能解决高数问题？

前向欧拉方程

常微分方程

常微分方程

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读