先日はSIRモデルについて書いたが、本県のコロナ動態を予測しようという大それたことを思った。素人考えであり、web上には引用しきれないほどの良質なサイトがあるので、きちんとした推定はそちらをごらんください。
例えば「COVID-19情報共有 — COVID19-Information sharing」。ただしwebではそれぞれの手法の妥当性について喧々轟々。私にはほとんど理解できません。
SIRモデルは三元の連立微分方程式で、
健常者
感染者（無症状含む）
回復者＋死亡者
Web情報を見ると、パラメータの推定方法はいろいろあって目が回りそうである。後退差分などで求めるらしいが：
本県70万人、発症者数がとても少ないので (54例、4 Apr 2020）、bの式の分子はほぼ I(t)とみて、
、
r についてはこの段階で上の方法は難しいので、推奨隔離期間の２週間をもとに、
とした。
（感染者数はわからないので、まだ十分初期とみなして発症者数で推定した。）
S(0) = 700,000 人、I(0) = 0.1 人, R(0) = 0 としてRのdeSolveパッケージのodeでシミュレートすると、
左が全体、右が最初のほうを拡大したものである。横軸は時間、縦軸は人口(指数表記)。青が健常者、赤が感染者、緑が回復者、水色は感染者数の実測値である。
うん。全然合っていないし、本県が崩壊するような気がするので少し修正を試みる。
SIRモデルは感染者と接触した瞬間に健常者が一定の割合 b で感染するということだが、本県の山がちな地形を考えると、それぞれの地域が十分隔離されており、ある場所で増加した感染者のウイルスは他の地域になかなか広がらないかもしれない。その広がりにくさの効果を考えてパラメータaを導入し、
健常者
感染者
回復者
としてみる（本来はパラメータ推定するべきだが、実測値の傾向に合うように自分でゴリゴリと合わせる。）
当てはまりはなかなかよさそうである(当てはめたのだから当たり前かもしれないが)。ピーク時には1e+05 (= 100,000人)の感染者がいるようである。ただし、無症状者を含んだ数字であることに注意。
基本再生産数は一人の感染者が瞬間的（この想定では一日）に感染させる人数
で、コロナは一般に2から3と報告されているから、まず妥当な推定になっているのではないか。
ただし、ほんの少し初期値やパラメータをずらしただけで相当違う挙動をするので、本当は感度や確率を考える必要がありそうです。
結論
今回はかなりラフな見積もりをしましたが、当分は感染のリスクがあるので、自覚症状がなくても拡散を防ぐことは大事そうです。
コード
# SIR-kochi# SIR-kochi
# k.i 2020-04-12
library (deSolve)
# parameters
days <- 400 # simulation periodpop <- 700000 # population size
(b <- 1/pop)
(r <- 1/14)
(R0 <- b/r * pop)
# reported 2.28 (2.06–2.52) by Zhang et al. (2020) Int. J. Infectious Deseases
parameters <- c(b, r)
initial <- c (S = pop, I = 0.1, R = 0)
times <- seq (0, days, by = 1)
a <- 0.14 # degree of isolation
SIR <- function (t, state, parameters){
with (as.list (c(state, parameters)),{
dS <- - a * b * S * I
dI <- a * b * S * I - r * I
dR <- r * I
list ( c(dS, dI, dR) )
}
)
}
out <- ode (y = initial, times = times, func = SIR, parms = parameters)
out1 <- as.data.frame (out)
# plot
xmax <- max(times)
ymax <- pop
matplot (x = times, y = out1[,2:4],type = "l",
xlim = c(0, xmax), ylim = c(0, ymax),
xlab = "Time", ylab = "Susceptible and Recovered",
main = "SIR Model",
lwd = 1, lty = 1, bty = "l",
col = c("blue","red","green")
)
# 以下略
