種別[statuses] cocolog:91890615
セクションJRF のひとこと
日時2020年05月09日
元URLhttp://jrf.cocolog-nifty.com/statuses/2020/05/post-2d4490.html

「志村けん問題」を JAGS…

「志村けん問題」を JAGS で解く…という記事を書いた。超有名人が感染したことがわかったら、感染者数はわかる以前より多いと見積もっても正しいか…という問題。「正しい」というのが私の結論で、JAGS によるシミュレーションをその証拠とする。
JRF 2020年5月9日

記事は↓。

《「志村けん問題」を 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

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

JRF 2020年5月9日

JAGS で「志村けん問題」を解いてみたらカッコイイとか思った。「志村けん問題」にアタリを付けながら、その前に、まず JAGS でモンティ・ホール問題を解いたりした([cocolog:91836101])。

[cocolog:91888253] で書いたように Kruschke 『ベイズ統計モデリング - R, JAGS, Stan によるチュートリアル 原著第2版』を読んで JAGS の勉強をしながら、この「志村けん問題」に取り組んだ。この本を読むのにとにかく時間がかかった。

JRF 2020年5月9日

……。

「その2」の結論で書いたように、JAGS は確率に関する Prolog みたいなものという認識で、それがおもしろいと思っていた。しかし、そのような使い方だと ESS が小さくなり、実質的に評価ができないことがわかった。

JRF 2020年5月9日

不等式を使う場合だけでなく、sk_problem24.jag のような不等式のトリックを使わないようなものでも、ESSがとても小さくなるのは意外だった。その「苦労」の結果、Prolog みたいに使うのは間違いだと結論せざるを得ない。

(そもそも Prolog 自体、万能推論機でないので、逆にそこが似てるのかもしれないが。)

JRF 2020年5月9日

……。

……。

R の coda パッケージにバグを発見。

coda.samples から得られる mcmc.list …ここでは sk_problem01.jag を実行して post.list という名前を付けたとする。

JRF 2020年5月9日

それについて plot(post.list[,"n.pos.sk"], trace=FALSE) とすると、内部では densplot というものが呼ばれる。重ねて plot するため、densplot に ylim を渡そうとするとエラーが出る。

JRF 2020年5月9日

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

JRF 2020年5月9日

英語だと…、

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

JRF 2020年5月9日

《coda/output.R at master - cran/coda》  
https://github.com/cran/coda/blob/master/R/output.R

JRF 2020年5月9日

で、解決法だが、グラフを重ねたければ ggplot2 を使えばいいということだと思う。そのほうがキレイだし。私はそうした。

JRF 2020年5月9日

後方参照 (2 件)