ベイズ推論の事後分布(離散)の更新について

パラメータの事前分布を自分で仮定して、観測データを元により適したパラメータを推定するのがベイズ学習といった話だった。
事前分布として固定値ではなく正規分布といった確率分布を与えてその時にちゃんとパラメータが学習されるのか確認する。

例としてある1枚のコインが存在して、そのコインで100回コイントスをする。
表が10回ほど、裏が90回ほど出たとする。このとき結果から1枚のコインの表が出る確率が求められるか。

100回中10回程度しか表が出てないんだから確率としては $\frac{1}{10}$ でいいじゃん＾＾と自分は思うがベイズ学習で似たような結果が得られるか確認する。

まずは上記の設定に似たコイントスの結果を生成はこんな感じ

true_mu = 0.16
X = np.random.binomial(n=1,p=true_mu,size=100)

# 結果の例
# (array([0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
#        0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
#        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
#        0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0,
#        0, 0, 0, 0, 0, 0, 0, 0]) )

コイントスは二値の値 $x \in$ { 0 , 1 } を取るので、{0,1}を生成してくれるベルヌーイ分布を使う。(1がコインの表とする)
ここではベルヌーイ分布に必要なパラメータ $\mu$ の真の値を 0.16 とした。この値に近いものを観測データから求めるのが目標。

前回の通りにやるとパラメータの推定は i.i.d な観測データ $X =$ { $x_1,x_2,...,x_N$ }としたとき、ベイズの定理から

$p(\mu |X) = \frac{p(X|\mu )p(\mu)}{p(X)} = p(\mu) \frac{\prod_{n=1}^N p(x_n | \mu)}{p(X)}$

とできる。

タイトル通り事前分布として固定値ではなく確率分布を与える。前提としてベルヌーイ分布のパラメータ $\mu$ の取りうる範囲は (0,1) の間でなくてはいけない。
つまり (0, 1) の間で実数を生成してくれるような確率分布を与えればいいんだな。

ベータ分布なるものは $\mu \in (0,1)$ となる実数を生成してくれる連続確率分布らしい。これを使ってみよう。定義は以下のようになっている

$Beta(\mu | a,b) = \frac{\Gamma (a+b)}{\Gamma (a)\Gamma (b)} \mu^{a-1} (1-\mu)^{b-1}$

a,bは非負の実数でなくてはならず、自分で設定するハイパーパラメータ。殺傷力ありそうな記号 $\Gamma()$ がついた関数はガンマ関数。階乗を表現してるそうだ。

a,bを適当に設定したベータ分布をプロットしてみた結果はこんなんになった

f:id:Owatank:20180412131347p:plain

ふわ＾〜

aとbの値が小さいと下に凹んだ形状で、大きいと正規分布みたいな山のような形になるのかな。いや数式見れば形状に説明つくんだろうけど・・・

観測データの見えない生成規則がベルヌーイ分布

$p(x| \mu) \sim Bern(x | \mu) = {\mu}^x (1-\mu)^{1-x}$

に従うとしているから、この式と上のベータ分布の式を事後分布の右辺に代入してみる

$p(\mu |X) = \frac{\Gamma (a+b)}{\Gamma (a)\Gamma (b)} \frac{(\prod_{n=1}^N {\mu}^{x_n} (1-\mu)^{1-x_n} ) \mu^{a-1} (1-\mu)^{b-1} }{p(X)}$

ヴワッ・・・果てしなく汚い式になってしまった・・・
ガンマ関数と分母はパラメータ $\mu$ によらないので、定数Aとおいて見ればこの式は

$p(\mu |X) = A { \mu^{(\sum_{n=1}^N x_n + a) -1} (1-\mu)^{(N - \sum_{n=1}^N x_n + b) -1} }$

こんな感じで書けるはずなのだ。ようは $\mu$ と $(1-\mu )$ でまとめられる。
この式の形はベータ分布と似ている。共役性とやらで実は $p(x| \mu)$ をベルヌーイ分布に設定して、事前分布 $p(\mu)$ をベータ分布にしたとき事後分布はベータ分布と同じ形の分布になるのだ・・・

何が嬉しいって分母、もとい周辺尤度 $p(X)$ を計算しなくて済む。
多分この周辺尤度ってパラメータ $\mu$ が取りうる範囲で積分、つまりは周辺化して得られるはずだから

$p(X) = \int p(X,\mu) d\mu = \int p(X|\mu)p(\mu) d\mu$

こうなるはずなんだよな・・・違うかな・・・
積分するものによっては難しいものあるんじゃなかろうか(´･ω･)

気を取り直して観測データの集合Xからパラメータ $\mu$ を推定する $p(\mu|X)$ を計算するには $p(x| \mu)$ をベルヌーイ分布にして、事前分布をベータ分布に設定すれば、ベータ分布のハイパーパラメータに観測データの情報を加えたものから得られることがわかった。

ベータ分布のハイパーパラメータを a = 0.1 , b = 0.1 に設定して、 $\mu$ が本当に推定されていくのか見ていく

for i in range(0,len(X),10):
    est_mu = 0
    est_mu = scipy.stats.beta.rvs(np.sum(X[0:i+10])+a[0], len(X[0:i+10]) - np.sum(X[0:i+10]) + b[0])
    print('観測データの数 : {0}件 、　パラメータμの推定値 : {1} '.format(i,est_mu))

f:id:Owatank:20180412150324p:plain

観測するデータを10件ごと増やしていくと50件以降から真のパラメータの値 0.16 に近くなっていくのが若干わかる
学習できていると見ていいんじゃなかろうか

もういっちょ別の確率分布を使ってパラメータの推定をやってみる
ベルヌーイ分布は二値に対してのみだったが、サイコロの目のようなKパターンの確率も推定してみたい
このようなK次元の確率分布としてカテゴリ分布なるものがある

$Cat(\boldsymbol{s} | \boldsymbol{\pi}) = \prod_{k=1}^K {\pi_k}^{s_k}$

$\boldsymbol{s}$ は $\boldsymbol{s} = (0,0,1,0,..,0)$ といった、K次元のうちk番目の値が１で他の値が０といったonehot表記であること、 $\boldsymbol{\pi}$ は $\pi_k \in (0,1)$ で、かつ $\sum_{k=1}^K \pi_k = 1$ となることが条件

ベルヌーイ分布でのパラメータ $\mu$ のように、このカテゴリ分布では $\boldsymbol{\pi}$ がパラメータになる
次元数をK=6 とし、 $\pi_k = \frac{1}{6}$ とすれば等確率のサイコロが表せるんだな

同様に

$p(\boldsymbol{s} | \boldsymbol{\pi}) \sim Cat(\boldsymbol{s} | \boldsymbol{\pi})$

としたとき事前分布 $p(\boldsymbol{\pi})$ はどうすればいいだろう

ベータ分布は $\mu \in (0,1)$ 、つまり変数1個しか生成してくれない。同じ (0,1) 区間で実数を多次元で、かつそれらの和が 1 となってくれるようなものを生成してくれる分布はないだろうかと悩むところでディクリレ分布というものがその役割を果たしてくれる。すげえ！！！

$Dir( \boldsymbol{\pi}| \boldsymbol{\alpha} ) = \frac{\Gamma(\sum_{k=1}^K \alpha_k)}{\prod_{k=1}^K \Gamma(\alpha_k)} \prod_{k=1}^K {\pi_k}^{\alpha_k - 1}$

$\boldsymbol{\alpha}$ の要素 $\alpha_k$ は正の実数が条件のハイパーパラメータ

何はともあれ事後分布に代入だ。観測データを $S =$ { $\boldsymbol{s_1} , \boldsymbol{s_2}, ... , \boldsymbol{s_N}$ }とすると

$p(\boldsymbol{\pi} |S) = \frac{p(S|\boldsymbol{\pi} )p(\boldsymbol{\pi})}{p(S)} = \frac{(\prod_{n=1}^N Cat(\boldsymbol{s_n} | \boldsymbol{\pi}))Dir( \boldsymbol{\pi}| \boldsymbol{\alpha} )}{p(S)} = \frac{(\prod_{n=1}^N{(\prod_{k=1}^K {\pi_k}^{s_n,k}}))\frac{\Gamma(\sum_{k=1}^K \alpha_k)}{\prod_{k=1}^K \Gamma(\alpha_k)}\prod_{k=1}^K {\pi_k}^{\alpha_k - 1}}{p(S)}$

ヴォエッ・・。こっちも汚くなったが分母とガンマ関数の部分は定数とみなせるのでAとして、総乗の部分はどちらもkのものだからまとめることができる。整理すると

$p(\boldsymbol{\pi} |S) = A\prod_{k=1}^K {\pi_k}^{\sum_{n=1}^N s_{n,k} + \alpha_k - 1}$

こんな感じで事後分布をディクリレ分布と似た形の式にみることができるのだ。共役性ありがてえ・・・

等確率でないサイコロを作って、そのサイコロの確率、もといパラメータを推定してみる
まずサンプルとなるインチキなサイコロを投げた結果を作る

weights = [0.1,0.1,0.35,0.1,0.25,0.1]
dice = [1,2,3,4,5,6]

obs_dice = np.random.choice(dice,p=weights,size=300)

# カテゴリ分布に投げるためにonehot表記する
obs_dice_onehot = []
t = np.zeros(6)
for d in obs_dice:
    t = np.zeros(6)
    t[d-1] = 1
    obs_dice_onehot.append(t)

np.random.choice()で任意の重みで設定した結果をsizeの数だけくれるようだ。
ここでweights = [0.1,0.1,0.35,0.1,0.25,0.1]が左からサイコロで１がでる確率、２が出る確率、として、これらの重みを結果から推定する

ディクリレ分布に必要なハイパーパラメータ $\boldsymbol{\alpha}$ の値を全て 0.2 とする
データを10件ずつ増やしていってパラメータ $\boldsymbol{\pi}$ の推定値がどう変わるか見ていく

for i in range(0,300,10):
    s = [0] * 6
    alpha = [0.2,0.2,0.2,0.2,0.2,0.2]
    for k in range(0,6):
        s[k] = np.sum(list(map(lambda x: x[k], obs_dice_onehot[0:i+10])))
    alpha = np.add(alpha,s)
    p = scipy.stats.dirichlet.rvs(alpha)
    print(i,p)