Mastodawn

Jupyter notebooksがどこかで公開されているとよいと思いました。

渡辺さんの本を読んでMCMCの計算結果からWAICを計算するコードを自分で書くと勉強になります。

読んでなくても(事後)予測分布を計算してプロットするところまでなら誰でもできるはず。

階層ベイズでない単純なケースならWAICの計算は簡単です。

しかし、ベイズ統計は気軽に階層化できる点が便利なので階層ベイズのケースでWAICが計算できないとありがたみが薄い。

階層ベイズのケースではWAICの計算には外に出て来ないパラメータについての積分$$
p(x|z)=\int dy\, p_1(x|y)p_2(y|z)
$$を繰り返し計算する必要があり、これが非常に重い計算になります。続く

MCMC chainの長さに比例した回数の数値積分が必要になります。

精度を高めるとめちゃくちゃ遅くなり、精度を下げると計算結果の信頼度が下がる。結構面倒です。

私の実力では、そういう計算を自動化することは困難で、階層ベイズモデルごとに高速数値積分の工夫が必要なのではないかと思いました。

渡辺澄夫さんのフォロワーは、階層ベイズのケースで実際に計算するための方法が普及していないことにもっと注意を払うべきだと思いました。

統計学では、理論的だけでは不十分で、実際に計算する方法を普及させて全世界で共有するところまで行くことが必要だと思います。

続く

私はJulia言語で遊んでいるので、階層ベイズの数値積分の大量の繰り返しが必要なケースでのWAICの計算をJulia+Mamba.jl+Cubature.jlでやってみました。コードは

で公開されています。試行錯誤をしまくったので私のパソコン上には少しずつ違うJupyter notebooksが大量にある。そのうちの一つです。

そこでやっていることはRStanを使っている

のモデル1の「(2)のStanコード」と本質的に同じです。そこでは数値積分をStanの側でシンプソンの公式を使ってやっています。R側で数値積分を実行すると遅いので、Stan側で数値積分のプログラムを書いている。

PyMCだとどうするのでしょうか？

私のケースでは、Julia言語なので計算速度のことを心配せずに全部Juliaでやれたので楽でした。しかも、全部、既存のパッケージですますことができた。

しかし、RStanのケースでは、高速化のために、Stanによる数値積分のプログラミングが必要になってしまいました。おそらく、そういうことは普通の人には気楽にはできないと思う。

階層ベイズのケースでのWAICを気軽にかつ高速に計算できる環境は私が知る限りまだ存在しないです。

あるなら教えて欲しいところです。作って発表すると感謝する人がたくさんいると思う。

私は最も簡単な階層ベイズや統計物理でのIsing模型やPotts模型のMCMCをJulia言語でやってみたのですが、MCMCの周辺では重い計算が必要であり、実用面で「高速化」の必要性を非常に納得できました。「速いは正義」です。

Potts模型のJulia言語によるシミュレーションのコードはたとえば

で公開されています。

わざと試行錯誤の様子を残してあります。クリーンなものだけを公開すると、気軽に汚いnotebookを公開するべきだという私の主義に反する。

個人的に、MCMCでベイズ推定できるソフトを用意し、渡辺澄夫さんの本またはウェブサイトの情報を見て、WAICを計算するコードを自分で書くという勉強の仕方は非常におすすめ。(私もやってみた。)

WAICの理論的導出部分は理解できなくても、正しいWAICの値を出すプログラムが書けるところまで各種定義と意味を理解していれば実用的に役に立つ知識になります。

WAICを計算するために必要な定義は渡辺さんの本の中には分散して書かれているので、それらを適切にピックアップできるようになっていれば、渡辺さんの本を実際の仕事で役に立てることができます。

さらに、プログラムを書いて、正しいWAICの値を計算する手続きをすでに知っていれば、各種定義を正しく認識できているということになり、理論的な議論の部分も理解し易くなると思います。