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の条件付き確率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と独立であるということを利用していきます。

Marsaglia法とBox-Muller法

一様分布からガウス分布を作るアルゴリズムとして有名なBox-Muller法というのがあります．
式としては、{z_1,z_2\in(0,1)}とした時に
{\displaystyle y_1 = \sqrt{\mathstrut -2\ln z_1}\cos 2 \pi z_2 \\ y_2 = \sqrt{\mathstrut -2\ln z_1}\sin 2 \pi z_2 }
という変換を行う．これがBox-Muller
ボックス＝ミュラー法 - Wikipedia


それと似たもので、一様分布する変数の範囲がちょっと変わって{z_1,z_2\in(-1,1)}としてやって、更に
{ z_1^2+z_2^2 \leq 1 }
を満たすものだけを取り出して（要するにそうじゃないやつを棄却して）
{ y_1 = z_1(-2 \ln r^2 / r^2)^{1/2}\\ y_2 = z_2(-2 \ln r^2 / r^2)^{1/2} }
という変換をしてやるのが、Marsaglia
Marsaglia polar method - Wikipedia, the free encyclopedia

なんか読むとBox-Mullerのほうがcosとsinの計算を並列でできる？から早いとかなんとか．まあ棄却する部分があるMarsagliaのほうが遅そうやなぁという雰囲気はあるよね．

というわけで実装

import numpy as np
import math as math
import time
import matplotlib.pyplot as plt

pears = 1000
##compare marsaglia with box-muller
def boxmuller():
    """
    using box-muller method
    return:points random sampling following N(0,1)
    """
    points = []
    for i in range(pears):
        z1 = np.random.uniform(0,1)
        z2 = np.random.uniform(0,1)
        y1 = np.power((-2 * math.log(z1)),0.5) * math.cos(2 * np.pi * z2)
        y2 = np.power((-2 * math.log(z1)),0.5) * math.sin(2 * np.pi * z2)
        points.append([y1,y2])
    return points

def marsaglia():
    """
    marglia method
    return:points random sampling following N(0,1)
    """
    points = []
    for i in range(pears):
        z1,z2 = 0.,0.
        r = 0.
        while 1:
            z1 = np.random.uniform(-1,1)
            z2 = np.random.uniform(-1,1)
            r = z1 * z1 + z2 * z2
            if r <= 1:
                break
        y1 = z1 * np.power(-2. * math.log(r) / r,0.5)
        y2 = z2 * np.power(-2. * math.log(r) / r,0.5)
        points.append([y1,y2])
    return points

methods = [boxmuller,marsaglia]
for i in methods:
    start_time = time.time()
    for k in range(10):
        i()
    end_time = time.time()
    print (end_time-start_time)/100,i

points = np.array(marsaglia())
plt.plot(points[:,0],points[:,1],"ro",alpha=0.5)
points = np.array(boxmuller())
plt.plot(points[:,0],points[:,1],"go",alpha=0.5)
plt.show()

結果

0.00240999937057 <function boxmuller at 0x04357D70>
0.00396000146866 <function marsaglia at 0x04357DB0>

やっぱboxmullerのほうが早いみたい．

f:id:dette:20150802093017p:plain

どっちも標準正規分布にはなってるみたい．

PRML 第5章

ニューラルネットワークのところを読んでいますが詰まりました．

ニューラルネットワークでは、最終的な出力が

{\displaystyle y_k(x,w)=\sigma(\sum^M_jw^{2}_{ik}h(\sum^D_iw^{1}_{ji}x_i)) }

という風になっていて、

そこで疑問なのが、
例えばKクラスの分類問題をニューラルネットワークに説かせようと思った時、
この最後の{\sigma}の部分をソフトマックス関数

{\displaystyle \sigma=\exp(a_k)/\sum_j\exp(a_j) }
とすれば良い
（ただし{a_k}は出力ユニットへの活性である）


という説明があるのだけれど、疑問がひとつ．

最終の出力ユニットの活性{a_k}は{a_k=\sum^M_kw^{2}_{kj}a_j}となっているのだから、
出力ユニットjへの入力は既に隠れユニットの出力{a_j}と重み{w_{kj}}の積の和になっているはずで、

そこからどうやって各々の{a_j}を計算しなおしてるんだろうか…という

なんか多分どうでもいい勘違いなんだろうけど凄くイライラ．

prml 第4章

PRML

先週から読み始めてようやく第四章までたどり着きました。噂には聞いていたのですが、ここまで重たいものだとは正直予想外で、なんとかかんとか頑張ってついていっているのか、はたまたわかったような気になっているのかは謎ですが、とりあえずちょっとづつ前進している感じです。

読んでいて思うのは、式変形はわかるんだけど、その意味合いを考えるとちょっとこんがらがるというか、冷静に考えると意味わかってないなって気づくところが多いということ。これはベイズの定理を使うのに慣れ始めてから少しづつ改善はしているのですが、やっぱりなぜこの式をここで使うのか？とかをちゃんと考えながら読まないと、ビショップさんの意図を汲み取れなくなっていくので、そこがしんどいです。

あと自分が理解できていないからかも知れませんが、微妙に天下り的に結果を先に与えていたり、すでにわかっているパラメータを条件付き確率の既知の部分から省いていたり、そういう微妙な罠がそこかしこにあって、記述は綺麗だけど意味を汲み取るのに時間が掛かる時があって、そこも苦労している点です。

まあわかってる人からすれば大したことはない省略だとは思うんですけども。モデルの選択の部分でモデル選択のMiが突然なくなった時とか、結構悩まされました。みんなどうやってクリアしてるのか気になるところです。

やっぱ誰か一緒に読む人がいたほうがいいなぁ…

あと、今まで研究室にあったやつを借りて読んでいたのですが、やっぱ高いけどこの本は買わないとダメだと決心がついたので、そろそろ買いに行こうかと思います。買ったら投げ出さんやろうという、退路を断つ作戦。そうでもしないと最後まで絶対やらんしな。