SIRモデルからはじめる微分方程式と離散時間確率過程（後編）

前の記事の続きです。

今回はSIRモデルを人数のまま扱い、確率過程で扱います。このことでモデルはより正確になって定量的になりますが、「時間がたったらどうなるのか？」などの定性的な理解は難しくなります。時間に関しては一日ごとに感染者数が発表されることを考えて離散時間にします。離散時間の場合はStanの得意領域です。連続時間の場合は確率微分方程式となります。Rでは{yuima}パッケージがよく使われるでしょうか。

それでは元の微分方程式を離散時間の確率モデルに変換します。 まずは適当な時間刻み幅Δtを持ってきて差分化します。

nは時点を表します。は時刻t-Δtから時刻tまでの間にSからIになった人数を表します。S人の人がΔtの間に、ある確率pでIになることを考えると、これは平均pSの二項分布に従うと考えることができます。平均pSは微分方程式モデルと対応していないといけないので、以下になります。

次に観測誤差についてです。感染者数であるIに観測誤差が入るとすれば、どのように入ると思いますか？僕はそれは病院に行く確率だと思いました。その病気の症状が重くてヤバい場合、1に近い確率で病院に行きます。まれに寝てれば治るとか、つらすぎて行けないとかあるかもしれないので、確率aで病院に行くとしました。

さてRでいくつか乱数を振ってシミュレーションしてみました。

library(reshape2)
library(ggplot2)

SIRsim <- function(beta, gamma, N, T, dt, a){
   S <- rep(0, T)
   I <- rep(0, T)
   R <- rep(0, T)
   S[1] <- N-1
   I[1] <- 1
   R[1] <- 0
   for (i in 2:T){
      d_StoI <- rbinom(1, S[i-1], beta/N*I[i-1]*dt)
      d_ItoR <- rbinom(1, I[i-1], gamma*dt)
      S[i] <- S[i-1] - d_StoI
      I[i] <- I[i-1] + d_StoI - d_ItoR
      R[i] <- R[i-1] + d_ItoR
   }
   I_obs <- rbinom(T, I, a)
   return(data.frame(T=1:T, S=S, I=I, R=R, I_obs=I_obs))
}

N <- 1000
T <- 50
dt <- 1
beta <- 1
gamma <- 0.3
a <- 0.9

set.seed(1234)
res <- SIRsim(beta, gamma, N, T, dt, a)
write.table(res, file='SIRmodel_stochastic-process.txt', sep='\t', quote=F, row.names=F)


N_seed <- 100
res_l <- lapply(1:N_seed, function(seed){
   set.seed(seed)
   res <- SIRsim(beta, gamma, N, T, dt, a)
   data.frame(seed=paste0('seed', seed), t=1:T, I=res$I, I_obs=res$I_obs)
})
res_df <- do.call("rbind", res_l)
res_melt <- melt(res_df, id=c("seed", "t"))

# p <- ggplot(data=subset(res_melt, variable=='I'), aes(x=t, y=value, group=seed, color=seed))
# p <- p + geom_line(alpha=I(1/5), size=0.5) + geom_point(alpha=I(1/5), size=2)
# p <- p + theme(legend.position="none")
# ggsave(plot=p, file="output/SP_I_seed.png", dpi=300, width=6, height=5)

p <- ggplot(data=subset(res_melt, variable=='I_obs'), aes(x=t, y=value, group=seed, color=seed))
p <- p + geom_line(alpha=I(1/5), size=0.5) + geom_point(alpha=I(1/5), size=2)
p <- p + theme(legend.position="none")
ggsave(plot=p, file="output/SP_I_obs_seed.png", dpi=300, width=6, height=5)

zero_count <- nrow(subset(res_df, t==20 & I==0))

• 8～10行目: 初期値はS=N-1人,I=1人,R=0人としました。
• 12～13行目: 二項分布で変化分を生成します。
• 18,27行目: 病院行く確率を0.9としてノイズを入れています。結局観測される人数はI_obsになります。
• 53行目: 後のグラフを見ると分かるのですが、今回の場合、乱数100個中22個の場合においてすぐにIが絶滅して何も広がらないことが分かります。

乱数の種ごとに色を変えてplotした結果（SP_I_obs_seed.png）は以下の通りです。

さて、上の29-31行目で乱数の種が1234の時に出力された結果をデータとして、パラメータをStanを推定してみます。とはいえ、Stanでは現状では離散的なパラメータが扱えませんのでsumming outしないといけないのですが、N以下の全ての自然数についてそれを実行するのは現実的に厳しいです。

ちなみにBUGSではできるかなーと思って以下のBUGSコードを試したのですが、「value of binomial I.obs[4] must be between zero and order of I.obs[7]」とかエラーが出て実行できませんでした。久保さんのブログ記事にもあるように、どうもそうは簡単にいかないような気がしてきました。

model {
   S[1] <- N-1
   I[1] <- 1
   R[1] <- 0
   p.ItoR <- gamma*dt
   for (t in 1:(T-1)){
      p.StoI[t] <- beta/N*I[t]*dt
      d.StoI[t] ~ dbin(p.StoI[t], S[t])
      d.ItoR[t] ~ dbin(p.ItoR, I[t])
   }
   for (t in 2:T){
      S[t] <- S[t-1] - d.StoI[t-1]
      I[t] <- I[t-1] + d.StoI[t-1] - d.ItoR[t-1]
      R[t] <- R[t-1] + d.ItoR[t-1]
   }
   for (t in 1:T){
      I.obs[t] ~ dbin(0.9, I[t])
   }
   beta ~ dunif(0, 5);
   gamma ~ dunif(0, 1);
}

今回のように感染した人数のように大きな人数が見込まれる時は二項分布を正規分布で近似してモデリングするのがよさそうです。

---

今回の記事に相当する研究を熱心にやられている方は統計数理研究所の斎藤正也先生です。プレスリリースはこちら。書籍もあります（以下の第5章のところ）。興味のある方は読んでみてはいかがでしょうか。またこんな催しものもありました。ご参考までに。

現象数理学入門