2015-09-22

PRML 第3章のEvidence近似

PRML Python

PRML第3章では、線形回帰モデルを扱っています。
その後半で、これまで固定だとしていた、重みwの分布 ${N(w|0,\alpha^{-1} I)}$ の精度パラメータ ${\alpha}$ と、実際に観測される値tのモデルである ${N(t|w^{T}\phi(x),\beta^{-1})}$ に現れるノイズの精度を表す ${\beta}$ の値もいい感じの数値にしちゃいましょうという話が出てきます。

これを真面目にやろうとすると、α,βに関する事前分布を定義してやって、事後分布である
${p(\alpha,\beta|t) \propto p(t|\alpha,\beta)p(\alpha,\beta)}$
を最大化する ${\alpha,\beta}$ を探すことになります。

じゃあαとβの分布はどうするの？ということになります。
ここで、仮定を追加して、われわれは、αやβについて何も情報を持っていない、とします。

これはどういうことかというと、どのαやβに対しても、だいたい同じぐらいの確率があるよね、という風に考えるということです。
なので、(α,β)=(0.2,0.1)だと思っているモデルM1も、(3.0,5.0)だと思っているモデルM2も、同じぐらいの確率で選ばれるよね、とするわけです。

こうなると、p(α,β)はどのαやβに対しても同じ値となるので、先ほどの最大化すべき関数が
${p(t|\alpha,\beta)}$
となり、これは単にwを積分消去すれば計算できて、
${p(t|\alpha,\beta)=\int{p(t|w,\beta)p(w|\alpha)dw}}$
となり、特に ${p(t|w,\beta),p(w|\alpha)}$ がガウス分布に従うと仮定しているので、これは解析的に解くことが出来て、とても嬉しいねという流れになります。これをEvidence近似とか経験ベイズとか第二種最尤推定とか言ったりします。第二種最尤推定と呼ばれるのは、普通の最尤推定がwの分布が存在していることまで仮定しているのに対して、今度はαやβも固定されているのではなくて、（平らという条件はついていますが）確率的な分布があると考えているからです。（これをもっと発展させて、αやβも分布をもたせることも可能で、それはPRML下巻の変分推論のあたりに書いてあります。）

さてでは今度は実際に解析的に解けたものを最大化していこうと思いますが、実はこれにも二つのアプローチがあります。一つは、今回とる解析的に解いて微分をする方法。もう一つがEMアルゴリズムを用いる方法です。EMアルゴリズムは、最大化したい関数が隠れ変数の積分で表されている時に使えるので、今回だとwは積分によって消えているので、EMアルゴリズムを適用することが出来ます。
でも今回は、せっかく解析的に解くことができるので、前者のアプローチを取ることにします。（計算すると、EMアルゴリズムと解析的に解く方法は一致することも確認できます。PRML下巻を参照）

どうやって解くのかを式変形を書いても良いのですが、それは何処にでも載っているので、ここではその気持について書きたいと思います。（自分が後で見て復習（あの時の考え方ってどうだったんだろうとか）に使える為）
まず大事なのが有効パラメータ数（well-determined parameter)の事です。有効パラメータとはカジュアルに言えば特徴量ベクトル ${phi(x)}$ によって計算できる行列 ${\beta*\sum_i{\phi(x)\phi(x)^T}}$ の $n$ 個の固有値 ${\lambda_i}$ のうち事前分布 $\alpha$ と比べて大きい値を持つものの数といえます。

この行列の固有値が小さい固有ベクトルというのは特徴量の空間において、あまり動きのない方向であると言えます。この固有値はデータ数が増えると単調に増加していくということを考慮すれば、まだデータが十分に無いためにその方向の情報量がない(別の固有ベクトルで十分に説明がつく) という状況を表しています。すなわち各データ ${phi(x)}$ がそちら方向の変動に乏しい、ということです。

一方で固有値が大きいというのは上記の逆で、その方向への変動が大きかったということを表します。これらを加味すると、行列の固有値のうち比較的数値が大きい物の数、すなわちデータを表すのに必要な次元の数を計算して、その数の次元程度でデータを予測するのが良さそうです。

これを表すのが有効パラメータ数 ${\gamma = \sum_i{\frac{\lambda_i}{\alpha+\lambda_i}}}$ になります。
これはｗの事前分布の精度αと、先ほどのデータの固有値（重要度）λとの比になっていて、αに比べて大きければ各∑の中身は1になり、小さければ0になります。なので、すべての和を取ると、今あるデータ点は、αとくらべて重要度が高いものが何個ぐらいあるか、という指標になります。

今回は、ガウス基底と多項式基底を用いた時に、

・普通にMAP推定をやった時の予測平均値
・Evidence近似によって尤度を最大化した時に得られる予測平均値
・その時の有効パラメータγの変化

をプロットするやーつを作ってみました。

まずはガウス基底から。基底関数には、ガウス基底を50次元用意しています。
訓練データは、y=sin(x/2)+ノイズ（精度パラメータ5.)で生成していて、グラフ上には青い点として表示しました。青の実線はノイズの無い目的の値y=sin(x/2)です。
f:id:dette:20150922144805p:plain

・データ点が多くなると、有効パラメータγも大きくなること
・MAP推定だけだと、選んだαとβに依存するので、今回初期の値としてα=0.1としているために、データ点にめちゃくちゃ引っ張られていること
・データ点が多くなると、両者は一致してくること
とかがわかると思います。

データ点が多くなればML推定に全部近づいていくので問題ないのですが、次元に比べてデータ点が少ない時（今回は基底に50個のガウス基底を用いていますから、10,20,30の時とかが顕著だと思います）に、Evidence近似では、データにひっぱられることなくうまく平滑化ができています。

つぎは多項式基底。10次式まで用意しました。こちらはy=sin(x)がターゲットです。
f:id:dette:20150922150110p:plain

データ点が少ない時は繰り返し回数が足りてませんね。雰囲気サンプル数が20以下の時は、γは最終的に0になりそうです。
（データが足りないからこれからは何も言えないので全部零にしちゃえということでしょうか）

2015-08-27

EMアルゴリズムについての殴り書き

PRML

EMアルゴリズムって何？

隠れ変数が存在するモデルに対して、モデル変数を変化させた時に尤度関数を最大化させる方法の事．

隠れ変数って何？

実際に観測されないけれど、観測される値がどういう分布に従うのかを決定する変数のこと．
具体例で言うと、

回帰問題の重み変数だったり、

クラス分け問題の時にその変数がどのクラスに入るのかを表したりする変数だったり．

ざっくりいうと別になくてもいいんだけどあったほうが観測値を説明しやすくするための変数

結局何が起こるの

argmax p(X|θ）．

ここでXは観測データの値で、θはモデルが持っている変数
（例：線形回帰問題の観測データ各点に乗っかるノイズの大きさβとか、重みwの事前分布の分散の大きさαとかそんなの）

じゃあ具体的には何をするの？

EステップとMステップと呼ばれている計算を交互にします．

Eステップって何

現在持っているモデル変数θ(old)を用いて、隠れ変数wの事後分布を計算します．

どういうことかというと、モデル変数がθ(old)で、それがわかっている時にdata:Xが得られる可能性を計算します．数式にするとp(Z|X,θ(old))になります．

（注：ここでoldとしているのはあとのMステップでこのθを更新するので、それと違いを表すために添字としておいています）

Mステップって何

Eステップで計算した事後分布p(Z|X,θ(old))を用いて、ある関数Q(θ,θ(old))を最大化します．

Qって何

完全データ対数尤度と呼ばれる関数を、先ほどのEステップで計算した事後分布で期待値を取った関数です．いみわからないと思うので数式を書いたほうが分かり良いと思います．

$Q(\theta,\theta_{old})=E_z[p(X,Z|\theta)]=\int p(X,Z|\theta)\times p(Z|X,\theta_{old})dZ$

2015-08-09

ガウス過程による事前分布からのサンプル

PRML Python

前提条件

まず初めに目的変数 $y$ は、ある重みのベクトル $w$ と、dataから得られた特徴ベクトル $\phi(x)$ の線形結合であらわされるとします。
要するに
$\displaystyle y = w^T \phi(x)$
となっているわけです。このとき、 $w$ の事前分布が
$p(w)= N(w|0,\alpha^{-1}I)$
で与えられているとすると、ガウス分布の線形和もガウス分布なので、yの分布もガウス分布となってくれます。

ガウス分布ということは、期待値と、共分散行列さえ求めてしまえばよいので、yの平均と共分散行列を求めると
$\displaystyle E[y]=\Phi E[w]=0\\ E[y y^T]=E[\Phi w (\Phi w)^T]=\Phi E[w w^T]\Phi = \Phi \alpha^{-1} \Phi^T=\alpha^{-1}\Phi\Phi^T$
となります。

これのイメージは、いくつかのxでのyの値がどういう分布になるのか知りたければ、今から知りたいと思うxの特徴ベクトル $\phi(x)$ をすべて計算して、それを掛け算して計算することのできるグラム行列 $K$ を計算すれば、yの同時分布がどういう多次元のガウス分布に従っているのか計算できますよ、ということです。

例えば $x=1,2,3$ でのyの値 $y(1),y(2),y(3)$ がどういう分布に従うのか計算したければ
1.まず各々の特徴ベクトルを計算して
2.グラム行列を作り
3.同時分布を得る

この同時分布が得られるというのが、ちょっと直感的にわかりづらいところだ個人的には思っていて、すべての点がこの場所にある確率がいくつ、という形で得られるので例えば $x\in[0,1]$ でどういう風なyの値になるのかを知りたいとすると、自分でxをある程度の数に区切って離散化して、そのグラム行列を計算してやる必要があります。なので厳密には連続的なyの分布を得ることはできません。（という理解でいます）

なんか要するに10次元とか、100次元とかの多次元ガウス分布を作ってそこからyを適当によってくると、それがガウス過程の事前分布から得られるあるサンプルのyの値ってことになるわけです。だから個人的にはPRML6章の図6.5のグラフはちょっとなめらかすぎて勘違いしそうだなぁ（というか僕が最初グラフの意味が分からなかったので…）とも思いました。はい。

ということを踏まえて図6.5を作るコードを書いてみました。
numpfspanっていうのが、離散化する数になっています。

# coding: UTF-8
import numpy as np
import matplotlib.pyplot as plt
class kernelfunction():
    def __init__(self,theta1,theta2,theta3,theta4):
        """
        kernel = Θ1*exp(-Θ2/2*||x1-x2||^2) + Θ3 + Θ4*<x1|x2>
        defaultはガウシアンカーネルを使用
        別途追加すれば別のカーネルに変更することもできるので
        """
        self.theta1 = theta1
        self.theta2 = theta2
        self.theta3 = theta3
        self.theta4 = theta4
        self.norm = "gaussian" # default
    def _kernel(self,x1,x2):
        if self.norm == "gaussian":
            return np.exp(-self.theta2/2.*np.inner(x1-x2,x1-x2))
        if self.norm == "ornstein process": #オルンシュタインウーレンベック過程
            return np.exp(-self.theta2*np.power(np.inner(x1-x2,x1-x2),.5))
    def kernel(self,x1,x2):
        """
        calculate kernel
        x1,x2: numpy.array, has same dimention
        """
        val = self.theta1 * self._kernel(x1,x2) + self.theta3 + self.theta4 * (np.inner(x1,x2))
        return val

kernel1 = kernelfunction(1,4.,0.,0.)#thata1=1,theta2=4 よくあるガウスカーネルになってます。
#kernel1.norm = "ornstein process"
numofspan = 100 #離散化する数 増やせばなめらかさは増しますが、計算コストも増えます。
gram_matrix = np.identity(numofspan)
x_range = np.linspace(-1,1,numofspan)
for i in range(numofspan):
    for k in range(numofspan):
        x1 = x_range[i]
        x2 = x_range[k]
        gram_matrix[i][k] = kernel1.kernel(np.array([x1]),np.array([x2]))
#plt.plot(range(len(gram_matrix[0])),np.abs(np.linalg.eig(gram_matrix)[0]))
#plt.show()
"""
#グラム行列の固有値がどうなってるか調べたかったので。
#本来はグラム行列が正定値である条件を満たしていないと、カーネルの条件に反してしまうので
#分布を作れないんですが、np.random.multivariate_normalは勝手に疑似逆行列を生成してくれるみたいで
#（正定値じゃないぞコラ）というエラーを吐きながらちゃんとサンプリングしてくれます。えらい！
"""

color = ["g","r","b","y"]
for i in range(10):
    y = np.random.multivariate_normal(np.zeros(numofspan),gram_matrix,1)
    plt.plot(x_range,y[0],color[i % len(color)])
plt.show()

(1.,4.,0.,0.)
f:id:dette:20150809050548p:plain
(1.,64.,0.,0.)
f:id:dette:20150809050551p:plain
(1.,4.,10.,0)
f:id:dette:20150809051244p:plain
どれもいい感じですね。傾向としては、ガウスカーネルのノルムにかかる係数が大きくなると、二つのxの値が離れた時に急速に0に共分散が近づいてしまうために、近いところのyの値にも相関がどんどんなくなっていることがわかります。逆に、近くないxにも大きな共分散を与えるような条件でやると、近い値どころか遠い値にも影響を及ぼしあって、あんまりyの値が変化しなくなっていることもわかります。

おまけ

オルンシュタインウーレンベック過程
これってウィナー過程の派生なんですかね？よくわかりません。
f:id:dette:20150809051100p:plain

2015-08-06

pythonのメッシュ切りのやり方

pythonでグラフを書きたいときに、[-1,1]の範囲で等間隔にメッシュを切りたいなと思うことがあると思いますがそれのやり方メモ。

## numpy.meshgridを使うやり方
## (始まり、終わり、個数)

X,Y = numpy.meshgrid(numpy.linspace(-1,1),numpy.linspace(-1,1))

"""
>>> X
array([[-5.        , -4.79591837, -4.59183673, ...,  4.59183673,
         4.79591837,  5.        ],
       [-5.        , -4.79591837, -4.59183673, ...,  4.59183673,
         4.79591837,  5.        ],
       [-5.        , -4.79591837, -4.59183673, ...,  4.59183673,
         4.79591837,  5.        ],
       ..., 
       [-5.        , -4.79591837, -4.59183673, ...,  4.59183673,
         4.79591837,  5.        ],
       [-5.        , -4.79591837, -4.59183673, ...,  4.59183673,
         4.79591837,  5.        ],
       [-5.        , -4.79591837, -4.59183673, ...,  4.59183673,
         4.79591837,  5.        ]])
"""


## numpy.mgridを使うやり方
## 0.10区切り
X,Y = numpy.mgrid[-1:1:.1,-1:1:.1]

"""
>>> X
array([[-1. , -1. , -1. , -1. , -1. , -1. , -1. , -1. , -1. , -1. , -1. ,
        -1. , -1. , -1. , -1. , -1. , -1. , -1. , -1. , -1. ],
       [-0.9, -0.9, -0.9, -0.9, -0.9, -0.9, -0.9, -0.9, -0.9, -0.9, -0.9,
        -0.9, -0.9, -0.9, -0.9, -0.9, -0.9, -0.9, -0.9, -0.9],
>>> Y
array([[-1. , -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1,  0. ,
         0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  0.9],
       [-1. , -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1,  0. ,
         0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  0.9],

        """

nykergoto’s blog

機械学習とpythonをメインに