前回２次元のガウス分布を条件付き確率により得られた１次元ガウス分布から推定する記事を書きましたが，今回はもっと単純な例で説明します．
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)
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
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,
}
A = 'mom'
B = 'dag'
N = 10000
for _ in range(N):
_A = sampleOne({k:v for k, v in probs.items() if k.endswith(B)}).split('|')[0]
_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)
ところでこの同時確率はギブスサンプリングを使わなくても解析的に解けるので，先程の事例の答え合わせができます．ベイズの定理より，下記が得られます．
は明に与えられていませんが，娘と息子が選ばれる内，娘が選ばれる確率であるので，
で先程求めた値(0.55)と近しい値であり，正しくギブスサンプリングが出来ていそうです．他の同時確率についても同様に計算ができます．
まとめ
サンプリング及びギブスサンプリングについて，Graham先生の資料を元に説明してみました．条件付き確率にうまく分解し，ギブスサンプリングを行うことで，元の確率に近しいものを獲得できることを，実例で見てみました．名前は仰々しいですが，意外と単純なしくみですよね．
今回のプログラムは下記にまとめています．
github.com
