nykergoto’s blog

機械学習とpythonをメインに

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

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

 

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

前提条件

まず初めに目的変数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

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],

        """

ギブスサンプリング:やってることのイメージ

統計の計算とかをやろうとすると、サンプリングという方法をとって計算をさせるという場面がよく起こります。

どういうときに使うのかというと、例えば確率密度関数{p(z)}に従う確率変数zを引数に取るある関数{f(z)}の期待値が計算したい場合などです。この場合計算するべきは
{ \displaystyle
E[f(z)]=\int f(z)p(z)dz
}
の計算になるわけですが、実際に計算するとなると、この関数がきれいな形をしていないと解析的に解くことが難しくなります。要するに不定積分をして数値を代入するだけでは解けなくなります。

もしどうしてもこれを計算したかったら、近似的にこの積分に一致するように計算を行っていく必要が出てきます。
一番シンプルに考えた時に思いつくのは「積分区間を均一に区切ってやって、その時の関数の値で近似を行う」という方法です。

しかしこの方法にも問題があります。それは確率変数がすごく大きな次元を取っていたときに、計算に必要な区切りの数が爆発してしまうという問題です。(次元の呪いと言います)
ですから全部の領域を区切っていく方法には限界があります。

そこで新たに考えられるのは、確率p(z)によってzをたくさん選んできて(これをz_i(iは何番目に選ばれたか)とあらわすことにします)、そのz_iを使って各々でf(z)を計算してやって、それの平均値をとるという方法です。数式にすると
\displaystyle
E[f(z)] \simeq \frac{1}{L}\sum_i^L f(z_i)
という風になります。こうするともはや計算するときにすべての領域を見る必要がなくなっています。また感覚的にもわかるように、選ぶ数をどんどん増やしていくとだんだんと求めたい値に収束していく、という性質も持ち合わせています。しかしこの方法も問題がないわけではなく、今度は全空間をなめる必要がなくなった代わりに、確率p(z)からzを選ばなくてはならなくなります。

もし確率がすべての空間で同じである(一様分布)ならば、適当な乱数をその範囲内の値から0の間をとるように発生させて、それをzとして採用すればよいですがふつう世の中にあるような問題を解こうとすると、もっと複雑な確率分布に従った現象は多く存在していますから、どうにかしてそのような難しい確率分布にそってzを取り出してくることが必要になります。

このzを取り出してくる操作のことをサンプリングと言い、いろいろなサンプリングの方法が提案されています。
今回はそのうちの一つである、ギブスサンプリングというものを作ってみました。

ギブスサンプリングはサンプリング法の中でもマルコフ連鎖モンテカルロ法(MCMC)という部分に分類されます。
マルコフ連鎖モンテカルロは、サンプル点を直前に取り出したサンプルに依存して新しいサンプルを取得します。これがマルコフ連鎖という名前のもとになっています。その中でもギブスサンプリングでは、その新しい点を取り出す操作をある次元にだけにそって行います。

別の言い方をすると、例えば3次元であればx方向だけを変える、と決めれば直前のy,zはそのままにして、xだけを新しく選びます。そしてその選ばれたx'を使って(x',y,z)を新しいサンプルとして採用します。この操作をいろいろな次元にそって行っていくことで、たくさんのサンプル点を次々に生成していきます。なので一回のサンプリングでは、変数は一つしか変わりません。

ではそのxはどう選ぶのかというと、次のxの条件付き確率p(x|y,z)によって選びます。

ギブスサンプリングが有効になるのは、この条件付き確率が簡単な場合です。なぜなら、結局は新しい点は、その条件付き確率(提案分布)から選んでこないといけませんから、この条件付き確率が難しい形をしていると意味がないからです。具体的な例で言うと、多次元のガウス分布などがあります。多次元ガウス分布はその条件付き確率も必ずガウス分布になることが知られていて、ガウス分布に従うサンプリングは、一様分布を用いて簡単に作成することができる(Marsaglia,Box-Muller法など)ので、一回の繰り返しは一様分布に従うサンプリングができれば簡単に行えます。ですからギブスサンプリングを用いれば、多次元ガウス分布によるサンプリングを行うことができるのです。*1*2

これを使って、2次元の場合にギブスサンプリングを行うプログラムをpythonで作成しました。

# coding: utf-8
__author__ = "nyk510"
"""
ギブスサンプリングを用いた, 二次元ガウス分布からのサンプリング
"""

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats


def gibbs_sampling(nu, cov, sample_size):
    """
    ギブスサンプリングを用いて与えられた共分散, 平均値を持つ
    多次元ガウス分布からのサンプリングを行う関数

    :param np.ndarray nu: 平均値
    :param np.ndarray cov: 共分散
    :param int sample_size: サンプリングする数
    :return:
    :rtype: np.ndarray
    """
    samples = []
    n_dim = nu.shape[0]
    # start point of sampling
    start = [0, 0]
    samples.append(start)
    search_dim = 0

    for i in range(sample_size):
        if search_dim == n_dim - 1:
            """
            search dimension selection is cyclic.
            it can be replaced random choice.
            """
            search_dim = 0
        else:
            search_dim = search_dim + 1

        prev_sample = samples[-1][:]
        A = cov[search_dim][search_dim - 1] / float(cov[search_dim - 1][search_dim - 1])  # A*Σ_yy = Σ_xy
        _y = prev_sample[search_dim - 1]  # previous values of other dimension

        # p(x|y) ~ N(x|nu[x]+A(_y-nu[y]),Σ_zz)
        # Σ_zz = Σ_xx - A0*Σ_yx

        mean = nu[search_dim] + A * (_y - nu[search_dim - 1])
        sigma_zz = cov[search_dim][search_dim] - A * cov[search_dim - 1][search_dim]

        sample_x = np.random.normal(loc=mean, scale=np.power(sigma_zz, .5), size=1)
        prev_sample[search_dim] = sample_x[0]
        samples.append(prev_sample)

    return np.array(samples)


if __name__ == '__main__':
    # 2 dimension normal distribution
    nu = np.ones(2)
    covariance = np.array([[0.5, 0.5], [0.5, 3]])

    # eig_values: 固有値
    # eig_vectors: 固有ベクトル
    eig_values, eig_vectors = np.linalg.eig(covariance)
    average_eigs = np.average(eig_values)
    sample = gibbs_sampling(nu, covariance, 1000)

    fig, ax1 = plt.subplots(figsize=(6, 6))
    ax1.scatter(sample[:, 0], sample[:, 1],
                marker="o", facecolor="none", alpha=1., s=30., edgecolor="C0", label="Samples"
                )

    # 答え合わせ
    # scipy.stats を用いて多次元ガウス分布の確率密度関数を計算
    multi_norm = stats.multivariate_normal(mean=nu, cov=covariance)
    X, Y = np.meshgrid(np.linspace(nu[0] - average_eigs * 2, nu[0] + average_eigs * 2, 1000),
                       np.linspace(nu[1] - average_eigs * 2, nu[1] + average_eigs * 2, 1000))
    Pos = np.empty(X.shape + (2,))
    Pos[:, :, 0] = X
    Pos[:, :, 1] = Y
    Z = multi_norm.pdf(Pos)
    ax1.contour(X, Y, Z, colors="C0", label="True Probability Density Function")
    ax1.legend()
    ax1.set_title("Gibbs Sampling")
    fig.tight_layout()
    fig.savefig("figures/gibbs_sampling.png", dpi=150)

f:id:dette:20170715071245p:plain
ちゃんと目標の分布に従ってサンプルが選ばれていることがわかります。

*1:
p(x,y) = N((x,y)|(\tilde{x},\tilde{y}),\Sigma)\\
ただし \Sigma = \begin{pmatrix}
\Sigma_{xx} & \Sigma_{xy}\\ 
\Sigma_{yx} & \Sigma_{yy}
\end{pmatrix}
\\
このとき条件付き確率\\
 p(x|y)=N(x|\tilde{x}+A_0(y-\tilde{y}),\Sigma')\\
ただし\\
\left\{\begin{matrix}
A_0\Sigma_{yy} = \Sigma_{xy}\\
\Sigma' = \Sigma_{xx} - A_0\Sigma_{yx}
\end{matrix}\right.
である

*2:これの証明はいろんな方法があるんですが、個人的に一番きれいだと思うのは特性関数を使ったものじゃないかと思います。スケッチとしては、z=x-\tilde{x}-A_0(y-\tilde{y})っておいてあげて、このzがyと独立であるということを利用していきます。