ベイズ推定①　モンテカルロ法

こんにちは、テックタヌキです。
普段はデータ分析のお仕事をしていて、学んだことのアウトプットをしています。
本記事では、ベイズ推定について学んだことをまとめます。
何か間違っているところがあれば指摘してもらえると喜びます。

はじめに

前回の記事では、EMアルゴリズムを使って尤度推定をしましたが、本記事ではベイズ推定を行います。

尤度推定は、未知のモデルパラメータ$\theta$に依存する観測データ$D=(x_1, \cdots, x_n)$の条件付確率分布であるモデル分布$p(D|\theta)$に従って生成されているとして、その尤度$L(\theta)=\prod_n p(x_n| \theta)$が最も大きくなるような$\theta$を推定することでした。

ベイズ推定でも、同様に未知のパラメータに依存するモデル分布$p(D|\theta)$が与えられますが、そのパラメータ$\theta$に対する事前分布$p(\theta)$を与えます。また、事前分布に対するハイパーパラメータ$\alpha$を与えて$p(\theta|\alpha)$とすることもあります。

ベイズ学習の枠組みでは、観測データDと確率モデル$\{p(D|\theta), p(\theta| \alpha)\}$が与えられます。そして、パラメータに関する事後分布$p(\theta|D, \alpha)$を求めるのが目的になります。

また、尤度推定ではパラメータに対して点推定を行っていましたが、ベイズ推定では、事後分布をもとめるので、分布推定の問題になります。^[1]

事後分布は、ベイズの定理から、以下のように求めることができます。^[2]

$$p(\theta|D, \alpha) = \frac{p(D|\theta)p(\theta|\alpha)}{p(D, \alpha)} = \frac{p(D|\theta)p(\theta|\alpha)}{\int p(D, \alpha, \theta) d\theta} = \frac{p(D|\theta)p(\theta|\alpha)}{\int p(D|\theta)p(\theta|\alpha) d\theta}$$

１次元ガウス分布の場合

ものは試しということで、観測値が１次元ガウス分布に従う場合を例に考えてみます。単純のために、分散$\sigma^2=1$は既知として、平均$\mu$を推定します。

また、事前分布を平均$\mu_0=0$, 分散$\sigma^2_0=1$の正規分布として設定する確率モデルを考えます。以上の準備の下、平均$\mu$の事後分布を求めます。

$$\begin{aligned} p(\mu|D, \sigma_0^2=1, \mu_0=1) &= \frac{p(D|\mu)p(\mu|\mu_0=0, \sigma^2_0=1)}{\int_{-\infty}^{\infty} p(D|\mu)p(\mu|\mu_0=0, \sigma^2_0=1)} \\ \\ &= \frac{\frac{\prod_n \exp(-(x_n- \mu)^2/2)}{ (2 \pi)^{N/2} }\frac{\exp(-\mu^2/2)} { (2 \pi)^{1/2} }} {\int_{-\infty}^{\infty} \frac{\prod_n \exp(-(x_n- \mu)^2/2)}{ (2 \pi)^{N/2} }\frac{\exp(-\mu^2/2)} { (2 \pi)^{1/2} }} \\ \\ &= \frac{\exp(- \sum_n(x_n- \mu)^2/2 - \mu^2/2)} {\int_{-\infty}^{\infty} \exp(- \sum_n(x_n- \mu)^2/2 - \mu^2/2)d\mu} \end{aligned}$$

ここで、

$$\begin{aligned} &\exp(- \sum_n(x_n- \mu)^2/2 - \mu^2/2) \\ &= \exp \Big(-\frac{N+1}{2}(\mu - \frac{\sum_n x_n}{N+1})^2\Big) \exp\Big(\frac{(\sum x_n)^2}{2(N+1)}- \frac{\sum_n x_n^2}{2}\Big) \end{aligned}$$

と分解できます。^[3] 元の式の分母の積分の中身は分子と同じで、積分は$\mu$についてだから、$\mu$と関係ない$\exp\Big(\frac{(\sum x_n)^2}{2(N+1)}- \frac{\sum_n x_n^2}{2}\Big)$は分母分子でキャンセルできます。よって、

$$\begin{aligned} p(\mu|D, \sigma_0^2=1, \mu_0=1) &= \frac{\exp \Big(-\frac{N+1}{2}(\mu - \frac{\sum_n x_n}{N+1})^2\Big) } {\int_{-\infty}^{\infty} \exp \Big(-\frac{2}{N+1}(\mu - \frac{\sum_n x_n}{N+1})^2\Big) d\mu}\\ \\ &= \frac{\exp \Big(-\frac{N+1}{2}(\mu - \frac{\sum_n x_n}{N+1})^2\Big) } {\int_{-\infty}^{\infty} \exp \Big(-\frac{N+1}{2}\mu^2\Big) d\mu} \\ \\ &= \frac{\exp \Big(-\frac{N+1}{2}(\mu - \frac{\sum_n x_n}{N+1})^2\Big) } {\sqrt{\frac{2}{N+1}}\int_{-\infty}^{\infty} \exp \Big(-\frac{N+1}{2}\mu^2\Big) d\mu} \\ \\ &= \frac{\exp \Big(-\frac{N+1}{2}(\mu - \frac{\sum_n x_n}{N+1})^2\Big) } {\sqrt{\frac{2 \pi}{N+1}}} \end{aligned}$$

となります。^[4]
よって、この式をグッと睨むと、$\mu$の事後分布は平均が$\frac{\sum_n x_n}{N+1}$で分散が$\frac{1}{N+1}$の正規分布に従うことがわかります。

この事後分布における確率が最大となるパラメータを求めると（MAP推定といわれます）、事後分布は正規分布のためそれは平均$\frac{\sum x_n}{N+1}$となります。改めて式に書くと、

$$\hat{\mu} = \text{argmax}_\mu p(\mu|D, \sigma^2_0=1, \mu_0=0) = \text{argmax}_\mu p(D|\mu)p(\mu|\sigma^2_0=1, \mu_0=0)$$

となります。^[5] 尤度推定では$p(D|\mu)$の最大化でしたが、MAP推定では事前分布$p(\mu|\sigma^2_0=1, \mu_0=0)$も考慮するようになっています。

最尤推定の結果と比べてみましょう。最尤推定の場合、分散既知に対する正規分布の平均の最尤推定量は標本平均である$\frac{\sum_n x_n}{N}$でした。一方で、今回のMAP推定では、$\frac{\sum_n x_n}{N+1}$となり尤度推定の結果と比較してやや小さく評価するようになります。

この理由は、$\mu$の事前分布に平均0の正規分布を与えたことで、絶対的に小さく見積もるように設計したためです。また、事後分布の分散は$1/(N+1)$のため、データ数に応じて柔軟にパラメータが中心$\frac{\sum_n x_n}{N+1}$からどの程度の範囲内にあるかを見積もることができます。^[6]

ベイズ推定の課題

上記のベイズ推定では、やや複雑ではありますが解析的に計算をすることができました。ベイズ推定の枠組みを再掲すると、以下のようになります。

$$p(\theta|D, \alpha) = \frac{p(D|\theta)p(\theta|\alpha)}{p(D, \alpha)} = \frac{p(D|\theta)p(\theta|\alpha)}{\int p(D, \alpha, \theta) d\theta} = \frac{p(D|\theta)p(\theta|\alpha)}{\int p(D|\theta)p(\theta|\alpha) d\theta}$$

改めて考えてみると、分母の積分が解析的に計算できますかどうかわかりません。また、数値計算の場合でも、積分の計算は重いです。^[7]

種明かしをすると今回の場合には、事前分布には、共役事前分布という性質の良い分布を選んだために、解析的に計算をすることができました。

解析的に計算できない例として、事前分布として、以下の標準コーシー分布
を与えてみます。

$$p(\mu|x_0=0, \gamma=1) = \frac{1}{\pi(1+\mu^2)}$$

事後分布の式に代入すると、

$$\begin{aligned} p(\mu|D, \alpha) &= \frac{p(D|\mu)p(\mu)}{\int p(D|\mu)p(\mu) d\theta} \\ &= \frac{\exp(- \sum_n(x_n- \mu)^2/2)/\pi(1+\mu^2)} {\int_{-\infty}^{\infty} \exp(- \sum_n(x_n- \mu)^2/2)/\pi(1+\mu^2) d\mu} \end{aligned}$$

のようになります。$\exp(\mu\text{の2次式})/(1+\mu^2)$の形のため、解析的には計算できそうにありません。
このようなとき、一般的には以下の２つの方法で解決します。ただし、事後分布を求めることは難しいので、パラメータに関する関数の値（たとえば平均や分散）を求めます。^[8]

• ①サンプル平均で近似する
• ②事後分布にできるだけ近い分布を選んで計算する

サンプル平均で近似する場合

$\mu$の事後分布を得ることはあきらめて、$\mu$の事後分布による期待値を求めます。しかし、事後分布はわからないのでベイズの定理を使って事前分布と尤度で計算します。すなわち、

$$\begin{aligned} \int \mu p(\mu|D) d\mu &= \int \mu \frac{p(D|\mu)p(\mu)}{p(D)}d\mu \\ &= \frac{\int \mu \exp\Big(- \sum_n(x_n- \mu)^2/2\Big)/\pi(1+\mu^2) d\mu}{\int \exp\Big(- \sum_n(x_n- \mu)^2/2\Big)/\pi(1+\mu^2) d\mu} \\ &= \frac{\int \frac{\mu}{1+\mu^2} \exp\Big(- \frac{N}{2}(\mu-\frac{\sum_n x_n}{N})^2\Big) d\mu}{\int \frac{1}{1+\mu^2} \exp\Big(- \frac{N}{2}(\mu-\frac{\sum_n x_n}{N})^2\Big) d\mu} \\ &= \frac{E_{\mu \sim N(\sum_n x_n /N, 1/N)}[\frac{\mu}{1+\mu^2}]}{E_{\mu \sim N(\sum_n x_n /N, 1/N)}[\frac{1}{1+\mu^2}]} \\ &\simeq \frac{\frac{1}{M}\sum_m \frac{\mu_m}{1+\mu_m^2}}{\frac{1}{M}\sum_m \frac{1}{1+\mu_m^2}} \end{aligned}$$

つまり、$\mu$の事後分布による期待値は、$N(\sum_n x_n /N, 1/N)$に従う確率変数による期待値の比であらわせます。^[9] 計算機上では期待値を求められないので、実際にはサンプルをМ回とったときの平均とします。これを、モンテカルロ平均と呼びます。

分散も同様にして、$\int (\mu- E[\mu])^2 p(\mu|D) d\mu$を求めればよいですが、分散公式を用いて求めたほうが簡単でしょう。そのために、$\int \mu^2 p(\mu|D) d\mu$は同様にして、以下のように計算できます。

$$\begin{aligned} \int \mu^2 p(\mu|D) d\mu &= \frac{E_{\mu \sim N(\sum_n x_n /N, 1/N)}[\frac{\mu^2}{1+\mu^2}]}{E_{\mu \sim N(\sum_n x_n /N, 1/N)}[\frac{1}{1+\mu^2}]} \\ &\simeq \frac{\frac{1}{M}\sum_m \frac{\mu_m^2}{1+\mu_m^2}}{\frac{1}{M}\sum_m \frac{1}{1+\mu_m^2}} \end{aligned}$$

分散は分散公式から、以下のように推定できます。

$$\sigma^2 = \int \mu^2 p(\mu|D) d\mu - \Big(\int \mu p(\mu|D) d\mu\Big)^2 \simeq \frac{\frac{1}{M}\sum_m \frac{\mu_m^2}{1+\mu_m^2}}{\frac{1}{M}\sum_m \frac{1}{1+\mu_m^2}} - \Big(\frac{\frac{1}{M}\sum_m \frac{\mu_m}{1+\mu_m^2}}{\frac{1}{M}\sum_m \frac{1}{1+\mu_m^2}}\Big)^2$$

シミュレーション

次の条件で、シミュレーションを行いました。

データ数:[5, 10, 50, 100, 200, 500, 1000, 2000]の８種類で生成
データ生成方法：未知の平均3、既知の分散が1の正規分布からデータを生成。
事前分布：標準コーシー分布
モンテカルロ平均に使うサンプルの数：1000

使用したコードは以下になります。

import numpy as np
import matplotlib.pyplot as plt
var_true = 1 # 既知の分散
mu_true = 3 # nuknownな平均

# 図の作成
fig, axs = plt.subplots(2, 4, figsize=(20, 10))
fig.suptitle('Monte Carlo Mean Simulation (1000 samples)', fontsize=16)

# シミュレーション
for iter, size in enumerate([5, 10, 50, 100, 200, 500, 1000, 2000]):
    # 疑似データ
    data = np.random.normal(loc=mu_true, scale=np.sqrt(var_true), size=size)
    mean = np.mean(data)

    mean_list = []
    var_lsit = []
    # モンテカルロ平均
    sample = np.random.normal(loc=mean, scale=1/np.sqrt(size), size=1000) # 乱数を取得
    for i in range(1000):
        _sample = sample[:i] # 推移を見たいのでサンプルから部分列を取得
        
        # モンテカルロ平均を求める式（累積和などを使えばもう少し効率がよさそうですが
        # 本質ではないので愚直にやります）
        montecalro_mean = np.mean(_sample/(1+_sample**2))/np.mean(1/(1+_sample**2))
        montecalro_meansq = np.mean(_sample**2/(1+_sample**2))/np.mean(1/(1+_sample**2))
        montecalro_var = montecalro_meansq - montecalro_mean**2

        # 履歴を格納
        mean_list.append(montecalro_mean)
        var_lsit.append(montecalro_var)

    # 可視化
    row = iter // 4
    col = iter % 4
    axs[row, col].plot(range(1,1000 + 1), mean_list, color = 'black', label='estimate of mean of μ')
    axs[row, col].fill_between(range(1,1000 + 1), np.array(mean_list)-3*np.sqrt(var_lsit), np.array(mean_list)+3*np.sqrt(var_lsit), alpha=0.3)
    axs[row, col].axhline(y=mu_true, color='r', linestyle='--', label='μ_true')
    axs[row, col].set_ylim(0, 5)
    axs[row, col].set_title(f'Data size: {size}')
    axs[row, col].set_xlabel('iter')
    axs[row, col].set_ylabel('val')
    axs[row, col].legend(loc='lower left')

次の図が、各データに対する未知パラメータ$\mu$に対する事後平均および事後分散をモンテカルロ平均によって求めた結果です。横軸がサンプル数、縦軸が値になります。青く塗っている領域は、求めた事後平均±3事後標準偏差になります。

図を見ると、いずれのデータ数でも平均、分散ともに収束していることがわかりますが、データ数が少ない場合には、真の平均3に対して、精度が悪い結果になりました。しかし、データ数が少ない場合にはその分、分散も大きくなっており、推定結果に自信がないことを表していることが考えられます。

ただし、データ数がある程度ある場合でも、取得したデータによって、真の値からずれる場合もあるようです。

まとめ

本記事では、分散既知の１次元ガウス分布の平均推定を行いました。
方法としては、共役事前分布を設定して解析的に事後分布を求めた場合と、共役事前分布ではない分布を事前分布にした場合にモンテカルロ平均によって、パラメータの期待値と分散を求めました。
次回の記事では、記事の中で言及していた変分ベイズの場合はどうなるかについて考えてみようと思います。

脚注

1. 尤度推定では、真のパラメータがあるという立場をとります。いわゆる頻度論です。ギリシャ哲学的にいえば、神様が決めた真のパラメータで定まる確率モデル（イデア）からデータが生成（イデアのコピー）されていると考えるようです。オブジェクト指向のクラスとインスタンスの関係に似ているかもしれません。一方でベイズでは、真のパラメータなどないという立場をとり、パラメータも分布として現れます。

2. $p(D, \alpha, \theta)= p(D|\alpha, \theta) p(\theta|\alpha)$として扱う場合もあるようです。本記事では$D, \alpha$は独立として扱います。

3. $\mu$について平方完成しただけです。

4. 最後の等式変形は、ガウス積分の結果です。

5. 分母の$p(D, \sigma_0^2=1, \mu_0=0)$は$\mu$を含まない定数のため、最大化の意味では無視できます。

6. 最尤推定の場合は真のパラメータは一つという立場をとっているので、区間推定の結果はこのようなシンプルな解釈にはなりません。統計WEB
   から引用すると、頻度論では以下のような直感的ではない解釈になります。「母集団から標本を取ってきて、その標本平均から母平均の95%信頼区間を求める、という作業を100回やったときに、95回はその区間の中に母平均が含まれる」ということを意味します。」

7. 1次元の場合はあまり重くないかもしれませんが、多次元の場合に区分求積法で求めた場合、分割数が急激に増えるため計算量が非常に重くなります。

8. デルタ法などを活用すると、信用区間も求められますが、今回は割愛します。

9. 本記事では尤度関数の分布に従う確率変数に対して期待値を取りましたが、事前分布に従う確率変数で期待値を取る場合もあるようです。どちらがいいかはわかりません。有識者の方は教えてください。