cocolog:91824909
新型コロナウィルス(COVID-19)関連のレポートの「追実験」を R や Python の練習としてやってみた。JAGS と Stan の使い方・変換の仕方がいまいちわからない。 (JRF 3031)
JRF 2020年4月13日 (月)
これは R で JAGS を用いているのだが、これを Python でできるようにするのを一つの目標にした。
とにかく、↑が動くように R のインストールからはじめる。JAGS、rjags は簡単にインストールできた。↑はすぐに動くようになった。
JRF2020/4/138151
……。
次に Python の側、pyjags をインストールするところでつまづく。pyjags は POSIX の環境でのみコンパイルできるとかなんとかで、Windows では現在インストールできないようだ。(たぶん)
諦めて、JAGS がダメなら Stan があるさ…と、pystan を入れる。しかし、conda istall pystan のみで↓でインストールしたのとほぼ同じ状況にできているようだ…。しかし example をやろうとするとエラーが出る。ImportError ではないが、↓を見て解決。
JRF2020/4/136083
《ImportError: DLL load failed: The specified module could not be found. · Issue #520 · stan-dev/pystan》
https://github.com/stan-dev/pystan/issues/520
JRF2020/4/132041
なんと、
素で、
<pre>
sm = pystan.StanModel(model_code=schools_code)
fit = sm.sampling(data=schools_dat, iter=1000, chains=4)
</pre>
…とかするとエラーが出るのに。
JRF2020/4/136812
<pre>
if __name__ == "__main__":
sm = pystan.StanModel(model_code=schools_code)
fit = sm.sampling(data=schools_dat, iter=1000, chains=4)
</pre>
…とかにするとエラーが出なくなった！
なんじゃそりゃ。
JRF2020/4/136958
……。
pystan なら使える。ならば、上のレポートの JAGS のコードを Stan に書き換え、それを PyStan に書き換えればいいだろうとあたりを付ける。ただし、上のレポートには「Stanは整数値を推定するのが少し手間」…と書いている。不吉な…しかし、やってみる。
JRF2020/4/136723
…まず、rstan を R に入れるところでつまづく。インストールはできるのだが、example コードを実行すると、Stan から g++ を呼ぶところでつまづいている。もちろん、Rtools を入れてもダメ。どうも "Program Files" というパスがダメそうな感じだが…。↓を見つける。
《rstan2.10.1が動かない！ | Sunny side up!》
https://norimune.net/2832
JRF2020/4/132569
…ということで、R を C:\Program Files 下ではなくルートの C:\ 以下にインストールしなおす。
JRF2020/4/136801
……。
上のレポートを stan に書き換えてみる。…しかしうまくいかない。どうも、Stan は parameters のなかに int を使えないらしい。generated quantities の中で、binomial_rng を使ってみると良いみたいな記事もあるが、上のレポートのモデルは、推論レベルで int を返す binomial を使うからうまくいかない。
「Stanは整数値を推定するのが少し手間」…どころか「できない」んじゃないか？
少なくとも今の私のレベルでは無理と判断し、これは今後の課題とする。…
JRF2020/4/137735
ちなみに…。PyStan も使えない。pyjags が入らない。…ので、「Python でできるようにする」目標については、PyperR を使って、R を Python から使う方針で行くことにする。PypeR 経由で rjags は使えたことから目標は一応達成ということに…。
JRF2020/4/130601
……。
……。
次に↓の記事について、
《「このままでは8割減できない」 「8割おじさん」こと西浦博教授が、コロナ拡大阻止でこの数字にこだわる理由》
https://www.buzzfeed.com/jp/naokoiwanaga/covid-19-nishiura
JRF2020/4/134585
こういうので使われるモデルって何なんだろう…、パーコレーションのコンタクト・プロセス？…とか考えていると、どこかで「SIR モデルだ」という記述を見る。Wikipedia には常微分方程式が書かれている。ほとんどしたことないが、常微分方程式のシミュレーションって興味ある。差分に直せばいいはずだが…とググっていると、↓を見つける。
JRF2020/4/137711
《感染病の数学予測モデルの紹介 (SIRモデル） - Qiita》
https://qiita.com/kotai2003/items/3078f4095c3e94e5325c
ここの Python コードは、そのまま動かせた。
JRF2020/4/137535
ここでは scipy.integrate の odeint を使っている。先に差分に直せばいいと書いたが、odeint を使わず簡単に差分にしたバージョンも作ってみると、それは odeint を使ったものと結果がほぼ一致した。
↓の記事によると SEIR モデルも使われているようだ。
《牧野淳一郎 新型コロナウィルスに関するメモ》
http://jun-makino.sakura.ne.jp/articles/corona/face.html
SEIR モデルも odeint で試してみたが、問題なくできるようだ。
JRF2020/4/136951
さらに牧野メモによると、↓のモデルが紹介されている。
《都道府県ごとのシミュレーションによる検討 – COVID-19情報共有 — COVID19-Information sharing》
https://www.fttsus.jp/covinfo/pref-simulation/
これは「佐藤レポート」として参照しよう。
JRF2020/4/138274
佐藤レポートでは、確率微分方程式が出てくる。確率微分方程式も Python で書いてみる…みたいな記事はあるが、dW(t - τ) のような遅れ要素を含むウィーナー過程をどう処理すればいいかは、ググっても見つからない。
しかたなく、dW[t] の列を先に発生させ、dW[t - tau] の乱数を使うということでやってみた。
JRF2020/4/134210
北海道の α と β の数値を使って確率微分方程式のシミュレーションのみ試してみる。いちおう、そうすると、グラフは佐藤レポートと同じような形になった。
ただし、I_0 から α と β の推定を行うことなどはやっていない。
JRF2020/4/135081
ただ、やってみて気になったのは、N の値。SIR モデルでは N = S + I + R は一定なのだけど、佐藤モデルでは遅れ要素があるため、dS に含まれる部分と dI に含まれる部分に若干の差がある。そのため N をプロットすると、確かに変化していることがわかる。「この程度」のことは別に気にしなくていいということだろうか？
JRF2020/4/133017
……。
あと、気になっているのが↓。
《志村けんのパラドックス - アスペ日記》
https://takeda25.hatenablog.jp/entry/2020/03/28/083257
これを JAGS か Stan で示せたらカッコイなとか夢想している。そこで、まずは JAGS と Stan の本を買って勉強したいな…と思っている。
JRF2020/4/133914
……。
追記。
おっと修正。
上で「conda istall pystan のみで↓でインストールしたのとほぼ同じ状況にできているようだ」と書いたが、「↓」の先を書いてなかった。それは↓。
《Pystanのインストール（anaconda環境） - Qiita》
https://qiita.com/Gyutan/items/c7d3c341efe09454a5e1
JRF2020/4/135729
これを見て追加でやったのは、conda install libpython m2w64-toolchain -c msys2 ぐらいだが、これは、conda install pystan をすればほぼすでにされているはず。distutils.cfg に付け加えろというコマンドもすでに付け加わっていた。
JRF2020/4/137092
……。
……。
追記。
「志村けん問題」については記事を書いた(↓+ [cocolog:91890615])。
《「志村けん問題」を JAGS で解く その１ 結論編》
http://jrf.cocolog-nifty.com/software/2020/05/post-918b0f.html
JRF2020/5/235767
で、↓を読んで離散確率分布を Stan で扱う方法のヒントを得た。
《西浦先生らによる実効再生産数の統計モデルを解説＆拡張する試み - StatModeling Memorandum》
http://statmodeling.hatenablog.com/entry/covid19-estimate-effective-reproduction-number
JRF2020/5/232529
…ので、この「ひとこと」の一番上の StatModeling Memorandum の記事の Stan 化に再挑戦してみる。
が、うまくいかない。
JRF2020/5/233984
<pre>
functions {
real continuous_binomial_lpdf(real k, real n, real p) {
real lp = if_else(n >= k, k * log(p) + (n - k) * log(1 - p)
- log(n + 1) - lbeta(n - k + 1, k + 1), 0);
return lp;
}
}
</pre>
JRF2020/5/235972
…みたいにして、binomial を使う替わりに continuous_binomial を使っていく。upper lower もできる限り定義して…。
JRF2020/5/237683
でも、うまくいかない。一番上の記事では n.pos.sum の平均が 7500 ぐらいなんだけど、私の Stan のモデルでは、5000 ぐらいになったり 270 ぐらいになったりで安定しない。4 つの chain で、50000 回ずつ回しても n.pos.sum の effectiveSize は 90 ぐらいで全々足りない感じ。
もっと回さないとダメなのかな…でも、これでも一時間以上かかるんだよね…。
JRF2020/5/238519
それとも何か間違ってるのかな？
JRF2020/5/235199
新型コロナウィルス関連でいくつかのレポートが出ている。
最初に興味を持ったのが、↓。
《COVID-19 日本国内の潜在的な陽性者数を推定する試み - StatModeling Memorandum》
http://statmodeling.hatenablog.com/entry/covid19-estimate-total-number-of-positives-in-japan
JRF2020/4/136652