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

この記事は何？
ベルヌーイ分布の学習
カテゴリ分布の学習
ポアソン分布の学習
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