ギブスサンプリング入門

前回２次元のガウス分布を条件付き確率により得られた１次元ガウス分布から推定する記事を書きましたが，今回はもっと単純な例で説明します． www.jonki.net

というのもGraham先生のチュートリアル資料（トピックモデル）にかなり分かりやすい例があったので，それを利用させていただいて説明します．チュートリアル資料の例題と対応しているので，それと合わせて本記事を見てみてください．

単純なサンプリング

まずサンプルする，ということはどういうことなのか，を説明します．これは何らかの確率分布から，その確率に応じてサンプル（標本）を抽出することを意味します．例えばある確率分布にはA, B, Cの３つの値がそれぞれ0.5, 0.3, 0.2で生起するような分布を考えます．そのような確率分布から１つののサンプルを得るには下記のようなコードになります．

import random

def sampleOne(probs):
    z = sum(probs.values())
    remaining = random.uniform(0, z)
    for k, v in probs.items():
        remaining -= v
        if remaining <= 0:
            return k

def main():
    probs = {
        'A': 0.5,
        'B': 0.3,
        'C': 0.2
    }

    N = 10000
    samples = {k: 0 for k in probs.keys()}
    for _ in range(N):
        sample = sampleOne(probs)
        # print('Try sample:', sample)
        samples[sample] += 1

    print('Result:')
    for k, v in samples.items():
        print('{}: {:.3f} ({}/{})'.format(k, v/N, v, N))

if __name__ == '__main__':
    main()

プログラムでは１万回サンプルを行い，A, B, Cがどれぐらいの確率で生起したか見てみましょう．だいたい元の分布の値と等しいですね．

$ python sample.py
Result:
A: 0.500 (5005/10000)
B: 0.302 (3016/10000)
C: 0.198 (1979/10000)

ちなみに上記コードのNを１００に変えて，１００回サンプルのときの例を見てみます．そうすると，先程よりも少しずれた値になってしまいました．このようにサンプリングでは基本的にサンプルする回数を多くすればするほど，元の分布に近しい値を獲得することができます．

Result:
A: 0.460 (46/100)
B: 0.300 (30/100)
C: 0.240 (24/100)

ギブスサンプリング

では，次にギブスサンプリングの説明をします．ギブスサンプリングでは，解析的に解けない同時確率を，条件付き確率に分解し，サンプリングを行うことで，元の確率（同時確率）を近似的に得ようというものです．その際に，各々の確率変量は同時には扱えないので，１つ以外を残して固定することで（観測値とする），サンプリングを行うものです．

例題として，Graham先生の資料を利用します．この問題では，親A（父or母）と子B（息子or娘）が買い物をしている状況を考え，それぞれの同時確率（例えばP(母，娘)）などをギブスサンプリングで求めます．Givenの情報として，AとBの条件付き確率は与えられているところからスタートしています．この例では同時確率は，解析的に計算できてしまう（後述）のですが，まずはギブスサンプリングで得てみましょう．

今回はAとBの２変量なので，Bを固定した状態でAをサンプル（A'〜P(A|B)），Aを固定した状態でBをサンプル（B'〜P(B|A）を交互に行います（〜は確率分布からのサンプルを表します）．そしてサンプルしたA'とB'が同時に買い物をしていた，としてカウントをしていき，各ケースの頻度を計算することで，擬似的に同時確率を求めます．このとき，AとBの初期値が必要になりますが，適当で良いです．

import random

# given probabilities
probs = {
    'mom|dag': 5/6,
    'fat|dag': 1/6,

    'mom|son': 5/8,
    'fat|son': 3/8,

    'dag|mom': 2/3,
    'son|mom': 1/3,

    'dag|fat': 2/5,
    'son|fat': 3/5
}

def sampleOne(probs):
    z = sum(probs.values())
    remaining = random.uniform(0, z)
    for k, v in probs.items():
        remaining -= v
        if remaining <= 0:
            return k

def gibbsSampling():
    samples = {
        'mom,dag': 0,
        'mom,son': 0,
        'fat,dag': 0,
        'fat,son': 0,
    }

    # initial values
    A = 'mom'
    B = 'dag'

    # sampling count
    N = 10000

    for _ in range(N):
        # Bを固定でAをサンプル
        _A = sampleOne({k:v for k, v in probs.items() if k.endswith(B)}).split('|')[0]
        # Aを固定でBをサンプル
        _B = sampleOne({k:v for k, v in probs.items() if k.endswith(A)}).split('|')[0]
        samples['{},{}'.format(_A, _B)] += 1

    print('Result:')
    for k, v in samples.items():
        print('P({}) = {:.3f} ({}/{})'.format(k, v/N, v, N))

if __name__ == '__main__':
    gibbsSampling()

これを実行すると，下記のようになりました．どうやらA=母，B=娘である確率が大きそうですね．

$ python gibbs_sampling.py
Result:
P(mom,dag) = 0.550 (5503/10000)
P(mom,son) = 0.278 (2776/10000)
P(fat,dag) = 0.111 (1107/10000)
P(fat,son) = 0.061 (614/10000)

ところでこの同時確率はギブスサンプリングを使わなくても解析的に解けるので，先程の事例の答え合わせができます．ベイズの定理より，下記が得られます．

$P(母, 娘) = P(母|娘)P(娘)$

$P(娘)$ は明に与えられていませんが，娘と息子が選ばれる内，娘が選ばれる確率であるので， $\frac{P(娘|母)+P(娘|父)}{P(娘|母)+P(娘|父)+P(息子|母)+P(息子|父)}=\frac{2/3+2/5}{2}=\frac{8}{15} \fallingdotseq 0.533$ で先程求めた値(0.55)と近しい値であり，正しくギブスサンプリングが出来ていそうです．他の同時確率についても同様に計算ができます．