作ったプログラムは以下においたよ。

https://github.com/PhysicsEngine/cpi-stats

ここで作成した線形回帰モデルのエビデンス評価をしてみた。 モデルエビデンスはざっくり言うと、訓練データを与える尤度を周辺化した値。つまり、そのデータを与える可能性が最も高いものを選ぶことができればパラメタの複雑さなどを自分で判断しなくてもよいというメリットがある。

今回、実装したのは基底関数を多項式とした場合で各目標値には加法性のガウスノイズがのる前提で計算を行った。

超パラメタalpha(モデルパラメタのばらつき)とbeta(ガウスノイズの精度)はそれぞれ0.0001、15.0にした。推定の方法としてはベイズ推定の予測分布を計算して行った。

// Mはモデルパラメタの数
val bayesRegression = new BayesRegression(M, 0.0001, 15.0)

モデルエビデンスの評価

基本的にはPRML(3.86)式に従った。

// Alphaはパラメタの分布の事前分布のばらつきを表す超パラメタ
// Betaはガウスノイズの精度を表す超パラメタ
// eMnの最大尤度での二乗和誤差
(M/2.0) * Math.log(Alpha) + (xlist.length/2.0) * Math.log(Beta) - eMn - 1/2.0 * Math.log(det(inv(_sN))) - xlist.length / 2.0 * Math.log(2 * Math.PI)

モデルエビデンスをモデルパラメタの複雑さ2~100の間でプロットした。赤線はM=4の位置。M=4で最大となっていることがわかる。というわけでそれぞれの主要な位置で実際のモデルのグラフをプロットしてみた

M=3の場合

なんか直線っぽく見えるけれど、二次関数です。ほぼあってないのでエビデンスもすごい低い。二乗和誤差関数の値も大きい。 このモデルはエビデンスを求めなくても選ばれなさそう。

M=4の場合

エビデンスが最大のモデル。見た目かなりよくfittingしている。

M=9の場合

むむ、M=4の場合と対して変わらない。

M=15の場合

これもあんまり変わらない。対数尺度で見た時に10程度の差しかないからか。この差はエビデンスでみると大して大きな差ではないのかもしれない。

M=30の場合

over fittingし始めた。M=15とM=30で大きく異るところはエビデンスの値というよりは訓練データ集合が25個だからそれより多いか少ないかの違いな気がする。

直感的には基底関数を多項式にしているので最尤推定で行うと訓練データ集合よりもモデルパラメータの数が多いとすべての点を通るようなグラフになる。今回はベイズ推定で行っているので事前分布のお陰で全点を通るなんてことはないけれど、それでも過学習をし始める分岐点としては訓練データの数になるかもしれない。

x=1付近でしか過学習してないのが気になるけれど。

M=100の場合

もう右側(x=1付近)のグラフはめちゃめちゃ。でもx=0付近はそんなにおかしくないのはなぜだろう。事前分布における精度は常に一定なのでそんなに場所によってばらついたりしなさそうなんだけど。

まとめ

評価値xによって過学習したりしてなかったりするのが気になるけれど、全体としてはエビデンス関数を最大化するようなモデルパラメタMを選ぶことによって、最適なモデルを作ることができそう。PRMLにも書いてあったけれど、「目で見りゃいいじゃん」というのはもっと高次元になった場合に通用しなくなるので、解析的にモデルを評価できる手法は有用。しかも割と時間もそんなにかからない。(Scalaでモデル100個の評価を行ったけれど1~2秒ほど)

超パラメタ、alphaとbetaに関しても評価してエビデンスを最大化してみるとよりよいモデルが得られるのかも。