人間だったら考えて

考えて考えて人間だったら考えて

「ベイズ推論による機械学習入門」を読んだので実験してみた (その1)

この記事は何?

ベイズ推論による機械学習入門」を読んでベイズ推論を理解した気がするので、本の3章で紹介されているアルゴリズムを実験してみました。
www.kspub.co.jp

ベルヌーイ分布の学習

本の3.2.1節で紹介されているベルヌーイ分布の学習と予測UCI Machine Learning Repository: Breast Cancer Wisconsin (Diagnostic) Data Setに適用してみます。

このデータセット乳がん診断の結果に関するデータセットであり、診断結果が良性か悪性かという情報が与えられています。
データセットのサイズはN=569であり、良性のレコード数は357、悪性のレコード数は212でした。

確率変数 x \in \{0, 1\}を、診断結果が良性のときにx=0、悪性のときにx=1と定義します。
2値を取る離散の事象を表現するために、ベルヌーイ分布を考えます。

 p(x|\mu)=\mathrm{Bern}(x|\mu)
ここで、\muはベルヌーイ分布のパラメータであり、\muの分布を訓練データX=\{x_1, \dots, x_N\}から推論します。

ベルヌーイ分布の共役事前分布であるベータ分布を考えます。

 p(\mu)=\mathrm{Beta}(\mu|a, b)
ここで、a,bはベータ分布のパラメータであり、事前に適当に与えることにします。

すると、事後分布p(\mu|X)は以下で与えられます、詳しくは本のP78を参照してください。

p(\mu|X)=\mathrm{Beta}(\mu|, \hat{a}, \hat{b}) \\ \hat{a}=\sum_{n=1}^{N}x_n + a \\ \hat{b}=N - \sum_{n=1}^{N}x_n + b

それではベルヌーイ分布の学習を実装してみます。といっても、\hat{a},\hat{b}を計算するだけですが……

def fit(X, a, b):
    # X: numpy array ([0, 1, 0, 1, ..., ])
    # a, b: hyper parameter of the prior distribution
    N = X.shape[0]
    hat_a = np.sum(X) + a
    hat_b = N - np.sum(X) + b
    return hat_a, hat_b

ベルヌーイ分布の学習を上のデータセットに適用します。事前分布のパラメータはa=2, b=2としました。
学習の結果、p(\mu|X)=\mathrm{Beta}(\mu|\hat{a}=214, \hat{b}=359)が得られました(上式に代入するだけで求められます)。
実際に事前分布・事後分布を見てみましょう、下図の左側は事前分布、右側は事後分布を表しています。
f:id:sz_dr:20171125233749p:plain
事後分布は\mu=0.39あたりの部分でピークが見られます。

カテゴリ分布の学習

本の3.2.2節で紹介されているカテゴリ分布の学習と予測
UCI Machine Learning Repository: Adult Data Setに適用してみます。

このデータセットは人口調査のデータセットであり、人種情報(White, Asian-Pac-Islander, Amer-Indian-Eskimo, Other, Black)が与えられています。
データセットのサイズはN=32561であり、各人種のレコード数は以下の通りです。

属性 サイズ
White 27816
Black 3124
Asian-Pac-Islander 1039
Amer-Indian-Eskimo 311
Other 271

確率変数 \mathbf{s}を各属性に対応する1-hotベクトルと定義します。
例えば、"White"は\mathbf{s}=(1, 0, 0, 0, 0)と表し、"Black"は\mathbf{s}=(0, 1, 0, 0, 0)と表します。

5値を取る離散の事象を表すために、カテゴリ分布を考えます。

p(\mathbf{s}|\mathbf{\pi})=\mathrm{Cat}(\mathbf{s}|\mathbf{\pi})
ここで、\mathbf{\pi}はカテゴリ分布のパラメータであり、\mathbf{\pi}の分布を訓練データS=\{\mathbf{s}_1, \dots, \mathbf{s}_N \}から推論します。

カテゴリ分布の共役事前分布であるディリクレ分布を考えます。

p(\mathbf{\pi})=\mathrm{Dir}(\mathbf{\pi}|\mathbf{\alpha})
ここで、\mathbf{\alpha}はディリクレ分布のパラメータであり、事前に適当に与えることにします。

すると、事後分布p(\mathbf{\pi}|S)は以下で与えられます、詳しくは本のP82を参照してください。

p(\mathbf{\pi}|S)=\mathrm{Dir}(\mathbf{\pi}|\hat{\mathbf{\alpha}}) \\ \hat{\alpha}_k=\sum_{n=1}^{N}s_{n,k}+\alpha_k

それではカテゴリ分布の学習を実装してみます。事前分布のパラメータ\mathbf{\alpha}に各属性のレコード数を足すだけです。

def fit(S, alpha):
    # S: numpy array, shape: (N, K)
    # alpha: hyper parameter of the prior distribution
    hat_alpha = np.sum(S, axis=0) + alpha
    return hat_alpha

\mathbf{\alpha}=(1, 1, 1, 1, 1)とすると、p(\mathbf{\pi}|S)=\mathrm{Dir}(\mathbf{\pi}|\hat{\mathbf{\alpha}}=(27817, 3125, 1040, 312, 272))が得られました。
\mathbf{\pi}は5次元(和が1になるという制約条件を考えれば4次元ですが)なので、事前分布・事後分布の描画はちょっと難しいです。。。

ポアソン分布の学習

本の3.2.3節で紹介されているポアソン分布の学習と予測
UCI Machine Learning Repository: Wine Quality Data Set
に適用してみます。

このデータセットはワインの質に関するデータセットであり、ワインの質を0から10の離散値で表しています。
赤ワインと白ワインのデータセットがありますが、今回は白ワインのデータセットを使います。
データセットのサイズはN=4898であり、ワインの質に対応するレコード数は下表の通りです。

ワインの質 レコード数
3 20
4 163
5 1457
6 2198
7 880
8 175
9 5

確率変数xをワインの質を表す離散値で定義します。
この事象を表現するために、ポアソン分布を考えます。

p(x|\lambda)=\mathrm{Poi}(x|\lambda)
ここで、\lambdaポアソン分布のパラメータであり、\lambdaの分布を訓練データX=\{x_1, \dots, x_N\}から推論します。

ポアソン分布の共役事前分布であるガンマ分布を考えます。

p(\lambda)=\mathrm{Gam}(\lambda|a,b)
ここで、a, bはガンマ分布のパラメータであり、事前に適当に与えることにします。

すると、事後分布p(\lambda|X)は以下で与えられます、詳しくは本のP85を参照してください。

p(\lambda|X)=\mathrm{Gam}(\lambda|\hat{a},\hat{b}) \\ \hat{a}=\sum_{n=1}^{N}x_n + a \\ \hat{b}=N+b

それではポアソン分布の学習を実装してみます。

def fit(X, a, b):
    # X: numpy array [6, 5, 9, ...]
    # a, b: hyper parameter of the prior distribution
    N = X.shape[0]
    hat_a = np.sum(X) + a
    hat_b = N + b
    return hat_a, hat_b

a=2, b=0.5とすると、p(\lambda|X)=\mathrm{Gam}(\lambda|\hat{a}=28792, \hat{b}=4898.5)が得られました。
実際に事前分布・事後分布を見てみましょう、下図の左側は事前分布、右側は事後分布を表しています。
f:id:sz_dr:20171126024639p:plain
事後分布は\lambda=0.59あたりの部分でピークが見られます。
このことは、\lambda_{\mathrm{MAP}}=\mathrm{arg~max}_{\lambda}p(\lambda|X)=\frac{\hat{a}-1}{\hat{b}}\simeq5.88となることからも確かめられます。

せっかくなのでサンプリングもやってみます。
ここでは比較のため、訓練データのヒストグラム・MAP推定の結果得られた\lambda_{\mathrm{MAP}}を用いたp(x|\lambda_{\mathrm{MAP}})からサンプリング・ベイズ推論による予測分布p(x_{\ast})からサンプリングの3つの結果を示します。
ここで、p(x_{\ast})=\int p(x_{\ast}|\lambda)p(\lambda)d\lambda=\mathrm{NB}(x_{\ast}|\hat{a}, 1/\hat{b})で与えられます、詳しくは本のP86を参照してください。

訓練データのヒストグラムp(x|\lambda_{\mathrm{MAP}})からサンプリング・p(x_{\ast})からサンプリングした結果を下図に示します。
f:id:sz_dr:20171128225514p:plain
サンプリングの結果、0~2や9~13の値も結構出てきています、もう少し尖って欲しいが…?

1次元ガウス分布の学習と予測

本の3.2.3節で紹介されている平均・精度が未知の場合のガウス分布の学習と予測UCI Machine Learning Repository: Wine Quality Data Setに適用してみます。

ポアソン分布の学習と予測で用いたデータセットですが、このデータセットにはワインの質以外にも、ワインのpHが与えられています。
pHについてのヒストグラムを下図に示します。
f:id:sz_dr:20171128231355p:plain

確率変数xをワインのpHを表す実数値で定義します。
この事象を表現するために、ガウス分布を考えます。

p(x|\mu, \lambda)=\mathcal{N}(x|\mu, \lambda^{-1})
ここで、\mu, \lambdaガウス分布のパラメータであり、\mu, \lambdaの分布を訓練データX=\{x_1, \dots, x_N\}から推論します。

ガウス分布の共役事前分布であるガウス・ガンマ分布を考えます。

p(\mu, \lambda)=\mathrm{NG}(\mu, \lambda|m, \beta, a, b)=\mathcal{N}(\mu|m, (\beta \lambda)^{-1})\mathrm{Gam}(\lambda | a, b)
ここで、m, \beta, a, bガウス・ガンマ分布のパラメータであり、事前に適当に与えることにします。

すると\muに関する事後分布は以下で与えられます。

p(\mu|\lambda, X)=\mathcal{N}(\mu|\hat{m}, (\hat{b}\lambda)^{-1}) \\ \hat{b}=N + \beta \\ \hat{m}=\frac{1}{\hat{\beta}}(\sum_{n=1}^{N}x_n+\beta m)

また、\lambdaに関する事後分布は以下で与えられます。

p(\lambda|X)=\mathrm{Gam}(\lambda|\hat{a},\hat{b}) \\ \hat{a}=N/2+a \\ \hat{b}=\frac{1}{2}(\sum_{n=1}^{N}x_{n}^{2}+\beta m^2 - \hat{\beta}\hat{m}^2)+b

これらの導出は本のP95を参照してください。

それではガウス分布の学習を実装してみます。

def fit(X, m, beta, a, b):
    # X: numpy array
    # m, beta, a, b: hyper parameter of the prior distribution
    N = X.shape[0]
    hat_beta = N + beta
    hat_m = 1 / hat_beta * (np.sum(X) + beta * m)
    hat_a = N / 2 + a
    hat_b = 1 / 2 * (np.sum(X ** 2) + beta * (m ** 2) - hat_beta * (hat_m ** 2)) + b
    return hat_beta, hat_m, hat_a, hat_b

m=0, \beta=1, a=1, b=1とすると、\hat{m}=3.19, \hat{\beta}=4899, \hat{a}=2450, \hat{b}=61.91が得られました。
実際に事前分布・事後分布を見てみましょう、下図の左側は事前分布、右側は事後分布を表しています。
f:id:sz_dr:20171130220905p:plain

\mu=3.19, \lambda=39あたりの部分でピークが見られます。
事前分布はべちゃーってなっていますが、事後分布はピシッてなっているのが良いですね(は?)

その他

事後分布が解析的に求まるなら実装も簡単でした、次は4章の混合モデルと近似推論について実装していきます