nykergoto’s blog

機械学習とpythonをメインに、雑多な内容をとりとめなく扱うブログです。

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

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

どういうときに使うのかというと、例えば確率密度関数{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で作成しました。

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

## 2dimention normal distribution
nu = np.ones((2))
covariance = np.array([[1,0.5],[0.5,3]])
# la:   固有値array
# v:    固有ベクトル array
la,v = np.linalg.eig(covariance)
avr_sigma = np.average(la) #あとでプロット範囲を決めるときの指標に

def gibssampling(nu,cov,sample_size):
    """
    Gibbs sampling 
    @nu     :average vector
    @cov    :covariance matrix
    @sample_size    :size of sample
    return type     :numpy.array
    length          :sample_size
    """

    samples = []
    dim = len(nu)
    # start point of sampling
    start = [0,0]
    samples.append(start)
    search_dim = 0
    for i in range(sample_size):
        if search_dim == dim-1:
            """
            search dimension select is cyclic
            it can replace randomly
            """
            search_dim = 0
        else:
            search_dim = search_dim +1
        #new-sampling

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

        # 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)

samples = gibssampling(nu,covariance,1000)
plt.plot(samples[:,0],samples[:,1],alpha=.1)

#答え合わせ
malti_normal = stats.multivariate_normal(mean=nu,cov=covariance)
X,Y = np.meshgrid(np.linspace(nu[0]-avr_sigma*2,nu[0]+avr_sigma*2,100),
                  np.linspace(nu[1]-avr_sigma*2,nu[1]+avr_sigma*2,100))
Pos = np.empty(X.shape + (2,))
Pos[:,:,0]=X
Pos[:,:,1]=Y
Z=malti_normal.pdf(Pos)
plt.contour(X,Y,Z,colors="k")
plt.show()

f:id:dette:20150806032152p: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と独立であるということを利用していきます。