cocolog:91890615
「志村けん問題」を JAGS で解く…という記事を書いた。超有名人が感染したことがわかったら、感染者数はわかる以前より多いと見積もっても正しいか…という問題。「正しい」というのが私の結論で、JAGS によるシミュレーションをその証拠とする。 (JRF 6476)
JRF 2020年5月 9日 (土)
……。
[cocolog:91824909] で書いたように新型コロナウィルス(COVID-19 )関連の記事(↓)を読んで、JAGS に興味を持った。
《COVID-19 日本国内の潜在的な陽性者数を推定する試み - StatModeling Memorandum》
http://statmodeling.hatenablog.com/entry/covid19-estimate-total-number-of-positives-in-japan
JRF2020/5/94294
JAGS で「志村けん問題」を解いてみたらカッコイイとか思った。「志村けん問題」にアタリを付けながら、その前に、まず JAGS でモンティ・ホール問題を解いたりした([cocolog:91836101])。
[cocolog:91888253] で書いたように Kruschke 『ベイズ統計モデリング - R, JAGS, Stan によるチュートリアル 原著第2版』を読んで JAGS の勉強をしながら、この「志村けん問題」に取り組んだ。この本を読むのにとにかく時間がかかった。
JRF2020/5/97733
……。
「その2」の結論で書いたように、JAGS は確率に関する Prolog みたいなものという認識で、それがおもしろいと思っていた。しかし、そのような使い方だと ESS が小さくなり、実質的に評価ができないことがわかった。
JRF2020/5/93397
不等式を使う場合だけでなく、sk_problem24.jag のような不等式のトリックを使わないようなものでも、ESSがとても小さくなるのは意外だった。その「苦労」の結果、Prolog みたいに使うのは間違いだと結論せざるを得ない。
(そもそも Prolog 自体、万能推論機でないので、逆にそこが似てるのかもしれないが。)
JRF2020/5/93026
……。
……。
R の coda パッケージにバグを発見。
coda.samples から得られる mcmc.list …ここでは sk_problem01.jag を実行して post.list という名前を付けたとする。
JRF2020/5/98036
それについて plot(post.list[,"n.pos.sk"], trace=FALSE) とすると、内部では densplot というものが呼ばれる。重ねて plot するため、densplot に ylim を渡そうとするとエラーが出る。
JRF2020/5/90438
<pre>
> densplot(post.list[,"n.pos.sk"], ylim=c(0,0.02))
plot(yhist, xlab = xlab, ylab = ylab, ylim = ylim.par, main = main.par, でエラー:
オブジェクト 'ylim.par' がありません
</pre>
JRF2020/5/91036
英語だと…、
<pre>
> densplot(post.list[,"n.pos.sk"], ylim=c(0,0.02))
Error in plot(yhist, xlab = xlab, ylab = ylab, ylim = ylim.par, main = main.par, :
object 'ylim.par' not found
</pre>
JRF2020/5/97852
《coda/output.R at master - cran/coda》
https://github.com/cran/coda/blob/master/R/output.R
JRF2020/5/98644
で、解決法だが、グラフを重ねたければ ggplot2 を使えばいいということだと思う。そのほうがキレイだし。私はそうした。
JRF2020/5/99205
記事は↓。
《「志村けん問題」を JAGS で解く その1 結論編》
http://jrf.cocolog-nifty.com/software/2020/05/post-918b0f.html
《「志村けん問題」を JAGS で解く その2 試行錯誤編
http://jrf.cocolog-nifty.com/software/2020/05/post-e84089.html
JRF2020/5/96159