Mathtodon#julialang

mathtod.online/@tetu/1041798

おお!素晴らしいです。仕事速すぎ。

メルセンヌツイスターを使っている部分を も採用している dSFMT というメルセンヌツイスターの後継に置き換えることができれば、確実に速くなると思います。

私の計測では、メルセンヌツイスターよりdSFMTは3倍くらい速いです。

2Dイジングでは実用的には $10^{11}\approx 2^{36.5}$ 回以上乱数を発生させるような計算をすることになるので、周期がそれより長い擬似乱数を使わないと、科学的に意味のある計算にならない危険性があります。

dSFMTならその点も大丈夫だし、速いです。

2Dイジング模型のモンテカルロシミュレーションは「擬似乱数の質と生成の速さ」が科学的に非常に重要なことを納得させてくれます。

以前は擬似乱数発生にかかる時間は無視できるほど短いと誤解していた。

mathtod.online/@antimon2/10224

深夜に

[(p,q,r,s,t,u) for p ∈ 1:5 for q ∈ p+1:6 for r ∈ q+1:7 for s ∈ r+1:8 for t ∈ s+1:9 for u ∈ t+1:10 if mod(p+q+r+s+t+u, 5)==0]

してしまった.

不完全ベータ函数
\begin{align*}
& B_x(a,b)=\int_0^x t^{a-1}(1-t)^{b-1}\,dt, \\
& I_x(a,b)=\frac{B_x(a,b)}{B_1(a,b)}
\end{align*}の数値計算のコードが2千行以上あった件

twitter.com/genkuroki/status/9

twitter.com/genkuroki/status/9

と R で使われている不完全ベータ函数の数値計算は所謂 TOMS 708 のアルゴリズムに基づいています。

github.com/JuliaLang/Rmath-jul

これが2千行以上ある。先人の苦労がよくわかる内容で面白いです。

@satie どんな感じっすか? 面白そうなら、ぼくも試してみようかな。 への翻訳を始める危険があるので避けている(時間が取られる)。

nbviewer.jupyter.org/gist/genk
に新たに付け加わったプロットの例。

ベイズ自由エネルギーBFEとWBICによるモデル選択はほぼ一致している。

WBICのプロットとBFE(Beyesian Free Energy)のプロットの形もほぼ同じになっている。

mathtod.online/media/FBRTK7uHv

nbviewer.jupyter.org/gist/genk
に新たに付け加わったプロットの例。

サンプルを分散を1に固定した正規分布で生成したときの、分散を1に固定したprior_0と固定しない場合のprior_1の比較。

KL情報量KL(真の予測誤差)が小さい方の事前分布の方が予測精度が高い。

分散を1に決め打ちしたprior_0を採用している法が予測精度がおおむね高い。

真の予測誤差KL(これは観測可能量だけから計算できない)とWAIC(これは観測可能量だけから計算できる)によるモデル選択は、それらが逆相関していることが原因で、違いが出る場合がある。違いが出ている場合はWAICが間違っているということ。

同じ色の部分が対応。黄色とcyanの部分でWAICが正しいモデル選択に失敗している。

mathtod.online/media/5xOap4t1X

nbviewer.jupyter.org/gist/genk
正規分布の共役事前分布(正規ガンマ分布)

に「事前分布選択」の節を付け加えました。事前分布をWAIC, LOOCV, WBIC, iWBICで選択する話。

事前分布で一部のパラメーターをある特定の値に固定してしまうこともできるので、これはパラメーターの一部を固定した部分モデルともとのモデルを比較して選択することの特別な場合とみなせる。

この場合には、真の予測誤差(KL情報量)やベイズ自由エネルギーのexplicit formulaeがあります。

たから、真の予測誤差による選択という正解をWAIC、LOOCVによる選択がどれだけ再現できるかもわかる。

あと、ベイズ自由エネルギーの正確な値によるモデル選択とWBICによるモデル選択の違いがどの程度かもわかる。

本当はベイズ自由エネルギーのサンプルを動かしたときの平均でモデル選択できればいいのですが、通常の実験観察ではサンプルは1つだけなので無理。

iWBICは一見 $c$ による補正にあたることをやっているように見えるのですが、私の解釈での実装では、うまく補正できませんでした。

統計学で使われる道具はへたに補正を採用すると、複雑な処理になって何が問題なのかわかり難くなる上に、補正が原因で信頼性が下がるというようなことがありがちなように思えます。

数学的にシンプルな道具をその欠点を承知の上で使う方が安全なような気がします。

具体的には2×2の分割表では何の工夫もしていないカイ二乗検定が結構頑健で良かったという話。

私はG検定は使う気になれないし、フィッシャーの正確検定は正確じゃないので使いたくないし、わざわざその正確じゃない方にカイ二乗統計量を補正するというのも意味がわからん。

添付画像は

nbviewer.jupyter.org/gist/genk

より。WBICとベイズ自由エネルギーの比較の例です。

同一の正規分布で生成した千個のサンプルの各々についてWBICと自由エネルギーを計算し、まとめてプロットしています。大雑把に
$$
\mathrm{WBIC} \approx (\text{自由エネルギー}) - c
$$のような関係になっています。ここで、$c$ は事前分布とサンプルを生成している真の分布だけで決まる定数です。$c$ はサンプルの取り方にはよらない。

この $c$ の正体は何だろうか?

答えを知っている人がいれば教えて下さい。

もしも $c$ で WBIC を補正することが可能ならば素晴らしいことだと思います。 mathtod.online/media/5eyKd_7Ch

nbviewer.jupyter.org/gist/genk
正規分布の共役事前分布(正規ガンマ分布)

は込み入った手計算の結果をJulia言語で書き直したものとモンテカルロ法による計算を比較して誤りが残る可能性を減らすことによって作成されました。

まだバグが残っている可能性があります。見付けたら、教えて下さい。

私が公開しているJupyter notebookをコピー改変したものを公開しても全然問題ありません。そういう自由がないと誤りを見付けてもらい難くなると思う。

nbviewer.jupyter.org/gist/genk
正規分布の共役事前分布(正規ガンマ分布)

これは正規分布の場合に共役事前分布を用いて、すべて手計算で公式を作ってしまうという企画の Julia 言語で書かれた Jupyter notebook です。

WAICもWBICもexact formulaを使って計算できます。

論文が見付からないので自己流でiWBICらしきものを実装してみたのですが、「ベイズ自由エネルギーの真の値と差の平均値と標準偏差の大きさ」を見る限りにおいて、WBICよりiWBICが優れているようには見えませんでした。私の実装が間違っている?

どうして論文が見付からないんだろうか?賞をもらうような研究なのだから、論文が公開されていないとおかしいと思うのだが。

nbviewer.jupyter.org/gist/genk
正規分布の共役事前分布(正規ガンマ分布)

内容は実質的に、正規分布に関する共役事前分布とベイズ自由エネルギーとWBICとWAICなどの公式集。

大学1年レベルの微積分の計算は結構大事。みんな計算練習をおろそかにして後で後悔する。

「論より証拠」
「論よりコード」

で渡辺澄夫さんの反論はとても気持ちが良いのですが、MATLABのコードしか公開していない点は、不利な点だと思いました。

RやPythonのようなみんなが使っている言語や のような無料で気楽に使える高速計算言語で検証用コードを公開した方がよかったと思う。

渡辺さん以外の誰かがコードを公開して代わりに拡散するといいと思いました。

それを実行することには以下のように沢山のメリットがあり、とてもおいしい。

* 共役事前分布について勉強になる。

* WAICや交差検証についても勉強になる。

* プログラミングの勉強にもなる。

* WAICや交差検証などの普及に貢献できる。

* 単なる個人的勉強の枠を超えて、Gelman et al. vs. Watanabe の論争に参加できる。いきなりメジャーデビューする感じ(笑)。

で、渡辺澄夫さん曰く

【CV と WAIC を参考にしてみてください(プログラム含む)。

PSISCV と WAIC を参考にしてみてください(プログラム含む)。】

watanabe-www.math.dis.titech.a

watanabe-www.math.dis.titech.a

プログラムの中身に目を通しました。matlabなので嫌だなと思って目を通さずにすませていたのですが、実際に中身を見たら面白かったです。

共役事前分布で求めた公式を使って計算しています。重いMCMC法を使っていません。軽い計算ですみます。

PythonやR(や )での同様のプログラムを公開拡散する人が出て来れば、WAICや交差検証の普及に貢献することになると思います。

社会貢献になる。ちょっと数学とプログラミングができるだけで、ちょっとした社会貢献をできるネタは大量に転がっていると思う。

MambaでのMCMCシミュレーションのために使う確率分布の型には一部の函数が実装されているだけで十分です。ただし、並列計算の都合で {at}everywhere で定義しておかなければいけません。詳しい説明が

mambajl.readthedocs.io/en/late

にあります。

Distrubutions パッケージで定義されている確率分布を使うのではなく、自前で定義した確率分布の型を使った方がMambaでのMCMCが速くなることがあります。

Distributions パッケージの MixtureModel type の実装は非常に遅いです。

とあるMCMCで MixtureModel type 使わずにすませたら、10倍近くMCMCが速くなりました。

nbviewer.jupyter.org/gist/genk

擬似確率分布を直接定義することによって高速化

の節を見て下さい。

高速計算が必要な場合に、Distributions パッケージの MixtureModel type の使用は要注意だと思います。

誰か原因を突き止めてプルリクエストを送るとよいと思います。

gist.github.com/genkuroki/eea1
Julia言語におけるやってはいけないタイプの定義の仕方

特別に型を定義せずに直接計算すると

2.832 ms, 8.85 MiB

で終わる計算が、ある方法で定義した型を経由して計算すると、

1.986 s, 130.77 MiB

もかかってしまう!

計算時間が670倍になり、メモリー効率も大幅に悪化。

実はこれをやらかしていることに気付いて、自分の馬鹿さ加減に嫌気がさした。

docs.julialang.org/en/stable/m

に書いてある注意を守れば、型経由で計算しても計算効率は悪化しません。

ここを修正したら、激遅で困っていたMambaでのとあるMCMCの速さが60倍になった。

今まで数分かかっていたのが、数秒に縮まった。

gist.github.com/genkuroki/8e97
逆温度βでのベイズ推定

分散1の2つの正規分布の混合モデルをDistributionsパッケージのMixtureModelで定義していたのですが、どうもその部分がMambaでのMCMCの律速段階になっているらしいことに気付いたので、確率モデルを直接書き下してみました。

MCMCのスピードが一桁上がりました。

Julia言語のMambaパッケージは計算速度面でかなりピーキーな感じ。注意しないと簡単に計算速度が一桁遅くなる。Mambaの問題ではなく、主にDistributionsパッケージの方の問題なのかもしれないと思いました。

DistributionsパッケージのMixtureModelの実装のどこに問題があるかを明らかにして改良すると喜ぶ人がいると思う。

今日作ったJuputer notebookのリスト訂正。以下の2つ。

nbviewer.jupyter.org/gist/genk
逆温度βでのベイズ推定
【逆温度βでのMCMCをJulia言語のMambaで簡単に行う方法。β→∞でベイズ推定が最尤法に近付くこと。】

nbviewer.jupyter.org/gist/genk
Bayes自由エネルギーの精密な評価
【事前分布のサンプルを使ったモンテカルロ積分で自由エネルギーを計算。分散を1に固定した場合の自由エネルギーの正確な表式に基いた計算。】