ポアソン分布とガンマ分布から事後分布を導出して,事後期待値や予測分布などを求める

この前までとある所で一般化線形モデルから始まるベイズモデリングの勉強に参加していました.一般化線形モデルを通して勉強したのは初めてかもしれません.ということでベイズモデリングについて体系的に勉強するために記事を書いてみようと思います.

この記事ではまず手始めとしてある例題をもとに,ベイズの定理を用いて事後分布を導出する過程について述べます.

こんな例題を考えてみる

この記事では割とよくある例題として,次のような観測値がポアソン分布,事前分布をガンマ分布という場合を考えてみます.

例題
ポアソン分布に従う5つの観測値[2,3,0,4,1]が得られたとする.またこのポアソン分布のパラメータ,すなわち平均は{\alpha=15,\beta=2 }のガンマ分布とする.このとき,

  • 事後期待値
  • 事後分散
  • 予測分布

を求める

得られる標本がポアソン分布に従うものとし,共役な事前分布としてガンマ分布を選んだときに事後分布もガンマ分布になる例を示します.

まず事前分布を決める

問題にあるように事前分布としてガンマ分布を考えます.ガンマ分布の密度関数は次のように与えられます.


\begin{aligned}
f(x|\alpha, \beta )=\frac{\beta ^ \alpha }{ \Gamma (\alpha)   }x^{\alpha-1}e^{-\beta x }
\end{aligned}
\alpha=15,β=2を与えた場合は次のようになります.
f:id:shu10038:20191222191203p:plain

観測値から尤度関数を求める

観測値[2,3,0,4,1]が与えられた場合は尤度関数は次のように与えられます.


f(x|\theta )= \prod_{i=1}^n  \frac{\theta ^{x_i} \exp(- \theta ) }{x_i! }
事前分布と尤度関数を合わせて描画すると次のようになります.
f:id:shu10038:20191222191155p:plain

事前分布と尤度関数から事後分布を求める

事後分布は次のようになります.


\begin{eqnarray*}f(\theta|x )&=& \frac{ f(x \mid  \theta ) f( \theta )}{ f(x) }  \\
&\propto& \prod_{i=1}^n  \frac{\theta ^{x_i} \exp(- \theta ) }{x_i! }\frac{\beta ^ \alpha }{ \Gamma (\alpha)   }\theta^{\alpha-1}e^{-\beta \theta } \\
&\propto &  \theta ^{ \alpha + \sum_{i=1}^{n} x_{i} -1 } e^{- n \theta -\theta  \beta}  \end{eqnarray*}
このように定数以外はガンマ分布同じ形状になりました.パラメータも次のように更新されます.


\begin{eqnarray*}\alpha' &=&\alpha + \sum_{i=1}^{n} x_{i}=25\\
\beta'&=&  n +  \beta=7
\end{eqnarray*}

事前分布と尤度関数と合わせて描画すると次のようになります.観測値をもとに,事前分布が更新されています.
f:id:shu10038:20191222191158p:plain

事後分布から事後期待値,事後分散,予測分布をもとめる

事後分布がガンマ分布になることが分かったので,ガンマ分布の期待値と分散の計算より,



\begin{eqnarray*}
E( \theta  \mid x) &=&  \alpha'  \beta' = 175\\
V( \theta  \mid x) &=&  \alpha'  \beta'^{2} = 1225
\end{eqnarray*}

予測分布は次のとおりです.


\begin{eqnarray*}
f(x^*|x ) &=&  \int_0^ \infty    { f(x^* \mid  \theta ) f( \theta\mid x )}  d \theta \\
&=&  \int_0^ \infty  \frac{\theta ^{x_i} \exp(- \theta ) }{x_i! }\frac{\beta ^ \alpha }{ \Gamma (\alpha)   }\theta^{\alpha-1}e^{-\beta \theta }d \theta   \\
&=& \frac{\beta^\alpha}{x^*!~  \Gamma (\alpha)}  \int_0^ \infty \theta^{x^* +\alpha -1} e^{-\theta(1+\beta) }d \theta  ~~~~~~  \ldots (1) \\
\end{eqnarray*}
ここで{\theta(1+\beta)=k }とおくと,{\frac{d \theta}{dk}=\frac{1}{1+\beta} }なので


\begin{eqnarray*}
式(1)&=&\frac{\beta^\alpha}{x^* !~ \Gamma (\alpha)}  \int_0^ \infty \left(\frac{1}{1+\beta}\right)^{x^* +\alpha -1} e^{-\theta(1+\beta) } \frac{1}{1+\beta} dk  \\
&=&\frac{\beta^\alpha}{x^*!~  \Gamma (\alpha)} \left(\frac{1}{1+\beta}\right)^{x^*+ \alpha }\int_0^ \infty k^{x^* +\alpha -1} e^{- k}dk ~~~~~~  \ldots (2) \\
\end{eqnarray*}

上式について,ガンマ関数の定義より,


\begin{eqnarray*}
\int_0^ \infty k^{x^* +\alpha -1} e^{- k}dk =  \Gamma (x^* +\alpha)
\end{eqnarray*}
であり,またガンマ関数の性質より

\begin{eqnarray*}
x^*!= \Gamma (x^* +1)
\end{eqnarray*}
なので,


\begin{eqnarray*}
式(2)&=&\frac{ \Gamma (x^* +\alpha)}{ \Gamma (x^*+1)  \Gamma (\alpha)}   \beta ^ \alpha (1+ \beta )^{-(x^*+ \alpha )} ~~~~~~  \ldots (3) 
\end{eqnarray*}
ここで,

\begin{eqnarray*}
 \beta = \frac{1-p}{p} 
\end{eqnarray*}
とおくと,

\begin{eqnarray*}
式(3)&=&\frac{ \Gamma (x^* +\alpha)}{ \Gamma (x^*+1)  \Gamma (\alpha)}   \left(\frac{1-p}{p}\right) ^ \alpha \left( 1+ \frac{1-p}{p} \right)^{-(x^*+ \alpha )}\\
&=&\frac{ \Gamma (x^* +\alpha)}{ \Gamma (x^*+1)  \Gamma (\alpha)}   \left(\frac{1-p}{p}\right) ^ \alpha \left(  \frac{1}{p} \right)^{-(x^*+ \alpha )}\\
&=&\frac{ \Gamma (x^* +\alpha)}{ \Gamma (x^*+1)  \Gamma (\alpha)} p^{x^*}  \left(1-p\right) ^ \alpha
\end{eqnarray*}
という具合に,負の二項分布になります.

参考

PyMC3をインストールする際に起きたトラブル色々(’19年9月秋の陣)

普段はPyMCを使っているんですが,とある勉強会でPyStanを推奨しているのでインストールしようとしたらバージョンの違いか何かでC++コンパイラが動かず.もう諦めてPyMC3で頑張ろうと思ったら今度はPyMCも謎のエラー.Anacondaを再インストールしてPyMC3を再インストールしようとしたら色々とトラブルが発生しました.これのトラブルシューティングを記録しておきます.
結論を先にいうと,Anacondaと全モジュールのバージョンを上げました.

インストールしたバージョンなど

まず最初に,こいつらを再インストールしました.

  • Anaconda3-5.2.0-Windows-x86_64
  • PyMC3 Version: 3.6

その後に次のような症状が出ました

主な症状

PCを複数台持っているのですが,それぞれ違う症状が出ました.

PC1

PyMC3をインポートしようとしたら次のメッセージが.

ImportError: ArviZ is not installed. In order to use plot_trace:
pip install arviz

一応こういうのを参考にしましたが,解決せず..
discourse.pymc.io
arviz-devs.github.io

PC2

pm.traceplotを実行しようとすると次のようなエラーが出てくる.

AttributeError: module 'theano' has no attribute 'gof'

一応こういうのを参考にしましたが解決せず..
github.com
stackoverflow.com

解決方法

Pythonのバージョンを上げずに,Anacondaをアップデート

conda update -n base conda

次にモジュールをアップデート

conda update --all

結局どれが効いていたのか分からないけど,とにかくこれで直りました.

エラーメッセージ;visual studio code '$ file ' cannot be resolved. please open an editor

トラブル

Visual studio codeでjupyterを動かしたあとDebugを走らせようとすると,

visual studio code '$ file ' cannot be resolved. please open an editor

というダイアログが立ち上がった.
f:id:shu10038:20190812174859j:plain

対処法

ダイアログにある「Opne launch.json」を開けたてから再びDebugしたら,何故か正常に動いた

エラーメッセージ「OverflowError: Exceeded cell block limit (set 'agg.path.chunksize' rcparam)」の対処

久々の更新

matplotlibを使ってみたらこんなメッセージが出た

OverflowError: Exceeded cell block limit (set 'agg.path.chunksize' rcparam)

該当するところは描画のところ

plt.plot([3万行くらいのリスト],[3万行くらいのリスト])

サイズのところが悪さをしていたらしい

plt.figure(figsize=(10, 6))
plt.subplot(2, 1, 1)

コイツらをコメント化したら治った

標準化と中心化

Pythonで標準化とか中心化はどうやってやるんだろうという話です.

結論だけ書きます.(そのうち合間を見て加筆します)

このページではこういうデータを使います
f:id:shu10038:20181006173233j:plain

標準化

from scipy import stats
x = df.iloc[:, 2:6].apply(stats.zscore, axis=0)
x.head()

f:id:shu10038:20181006173231j:plain

標準化されているかを調べます
f:id:shu10038:20181006173229j:plain
平均ほぼ0,分散(標準偏差)ほぼ1

中心化

y = df.iloc[:, 1].apply(lambda v:v-df.iloc[:, 0].mean() ).astype(float)
y.head()

バイアス-バリアンス分解の導出を真面目にやる

どうもです.

只今,下記の書籍にを用いて勉強中です.

応用をめざす 数理統計学 (統計解析スタンダード)

応用をめざす 数理統計学 (統計解析スタンダード)

P118のあたりに推定の偏りの話や,不偏推定量の事が述べられています.そして推定量の真の母数と推定量の二乗誤差の期待値を計算する所では次の式が書いてありました.

f:id:shu10038:20180920215318j:plain
さぁ,なんということでしょう~,式の行間が開きすぎてよくわかりません.他の書籍で調べてもこんな感じで端折ってありますし,ネットで調べても「期待値の計算すれば上手いこと消えてほしい項が0になってくれます」みたいな感じで,きちんとと書いている記事が直ぐには見当たりませんでした.
なので,折角なので勉強がてらに式の展開をしてみます

展開

二乗誤差を計算すると

f:id:shu10038:20180920215319j:plain
ここで,第三項について, \theta  E\left(  \hat {  \theta } _{n} \right)  は定数であることを考慮すると,
f:id:shu10038:20180920215320j:plain
よって,元の式は
f:id:shu10038:20180920220322j:plain

Windows環境下でPython版のglmnetを使うとどうなるか?

夏のある日のこと,
私はLasso回帰を用いてあるデータを分析しようとしていた.

glmnetを呼び出して実行したところ,

x input must be a scipy float64 ndarray

とエラーメッセージが帰ってきた.

あれ?っと思ったけど,その日に扱っていたのは整数値だったことを思い出した.
しゃぁ無いから浮動小数点数にキャストしてみたんだ.
これで間違いは無いはず.
スクリプトを実行してみると・・・

raise ValueError('loadGlmlib does not currently work for windows')
ValueError: loadGlmlib does not currently work for windows

なん・・・だと・・・・
これwindows環境では動かんのかい!!
諦めて,教えてもらったAzureNotebookを使うかな.