#julialang #jupyter #pyplot #matplotlib #baysian #mcmc #mamba
Julia言語による混合正規分布モデルの場合の最尤推定とベイズ推定にフィッティング (AIC も WAIC も計算してある)
https://gist.github.com/genkuroki/e99426b5968fb2ef9dc95573d33c6391
添付画像の説明
・ヒストグラムはテストサンプルのヒストグラム。
サンプルは混合正規分布 0.5N(0,1)+0.5N(1,1) で生成されている。サイズはN=100.
・オレンジの線はベイズ推定による予測分布
・見難い緑の点線は最尤推定による予測分布
モデルのパラメーターは a,b でモデルは
(1-a)N(0,1)+aN(b,1).
#julialang #jupyter #pyplot #matplotlib #baysian #mcmc #mamba
混合正規分布のモデルで最尤推定とMCMCによるベイズ推定で色々遊んでみたのですが、そのケースでは大体以下のようになっている感じでした。(注意：事前分布は適当に一様分布にしているので、尤度函数＝事後分布と思って構いません。)
(1) 事後分布が1つの狭い山になっているケースでは、最尤推定とベイズ推定の結果の違いは微小で、MCMCによる推定の収束も速い。(iterationsの回数を大きくする必要が全然ない。)
(2) 事後分布の形が奇妙な形になるケース(モデルの特異点が事後分布の形に見えているケース)では、MCMCでベイズ推定を収束させることがとても大変になる。
こういうのはやってみないとわからないことが結構たくさんある感じ。
百聞は一見に如かず。
#julialang #jupyter #pyplot #matplotlib #baysian #mcmc #mamba
あと、事後分布の形は surf() で 3D plot するよりも、pcolormesh() で色だけでをプロットした方がわかりやすいということにも気付きました。
https://mathtod.online/@genkuroki/382542
に投稿した事後分布の pcolormesh() によるプロットを添付画像に投稿しておきます。
それぞれサンプルサイズはN=1000とN=10000なのですが、正則モデルで使える中心極限定理の効果が全然見えて来ていません。
$N\to\infty$ で消える効果が有限の $N$ で見えている。有限の $N$ で見えているのは「特異点」の効果です。
#julialang #jupyter #pyplot #matplotlib #baysian #mcmc #mamba
ちなみにサンプルサイズ N=10^5 でもこんな感じで「1つの山」になるには程遠い感じ。
サンプルを生成した確率分布は
0.5N(0,1)+0.5N(0.1,1)
で推定のための確率モデルは
(1-a)N(0.1)+aN(b,1)
です。 https://mathtod.online/media/ceepuRfQb1IOUICS30Y https://mathtod.online/media/qiobfEZm36ZeOJt-pdA
#julialang #jupyter #pyplot #matplotlib #baysian #mcmc #mamba
まだ試していないこと。
事後分布がプロットできているくらいだから、それを使って、MCMCに頼らずに直接積分を計算できるケースもあるように思えました。
誰かやってくれないですかね。
#julialang #jupyter #pyplot #matplotlib #baysian #mcmc #mamba
今朝、サンプルサイズ N=10^6 の事後分布＝尤度函数 (事前分布は一様分布)のグラフを描かせる函数を起動してから出勤した。
帰宅すると計算が終わっているはず。終わって無かったり、トラぶっていたらがっかりなのだ。
ところで、コンピューターで計算する話はコードを公開してくれないと、内容が曖昧になってとてもわかりにくいと思う。
ただし、昔から「困ったな」と思うケースも結構ある。matlabのコードだと公開してくれていても、matlabを使えないので困る。
個人的にはmatlabが駆逐されて、juliaが流行ればよいと思っています。
たぶん、juliaを宣伝している人達の多くがmatlabを嫌い。
#julialang #jupyter #pyplot #matplotlib #baysian #mcmc #mamba
Dropbox の共有フォルダに今朝走らせて来た函数のJupyter notebookが自動保存されていた！
無事計算が終わっていました。
添付画像がサンプルサイズN=10^6=百万の事後分布＝尤度函数のグラフです。
N=10^6でも真のパラメーターに集中せずに島ができている。