正規分布は非常に多くのシーンで、最適な確率分布として利用されています。今回は、正規分布のパラメータ $\mu$ と $\sigma^2$ を、ベイズ推定の枠組みで学習して予測を行います。
過去にベルヌーイ分布とポアソン分布でベイズ推定した記事の続きとなります。正規分布の場合はパラメータが $\mu$ と $\sigma^2$ の2つあるので、少し難易度が高くなりますが、逐一丁寧に導出していきます。
本記事の内容
- 分散既知の場合の平均 $\mu$ のベイズ推定
- 事後分布の導出（正規分布 × 正規分布 → 正規分布）
- 予測分布の導出
- Python での実装と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
正規分布の確率密度関数
まず復習として、正規分布の確率密度関数を示します。
$$ \begin{equation} p(x|\mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right) \end{equation} $$
ここで $\mu$ は正規分布の平均、$\sigma^2$ は分散です。
問題設定
データ $\mathcal{D} = \{x_1, x_2, \dots, x_N\}$ が正規分布から生成されたと仮定し、ベイズ推定でパラメータを学習します。
今回は簡単のため、分散 $\sigma^2$ は既知として、平均 $\mu$ のみをベイズ推定します。
尤度関数
データ $\mathcal{D}$ の尤度関数は、各データ点が独立に生成されたと仮定すると、
$$ \begin{align} p(\mathcal{D}|\mu) &= \prod_{n=1}^{N} p(x_n|\mu, \sigma^2) \\ &= \prod_{n=1}^{N} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x_n – \mu)^2}{2\sigma^2}\right) \\ &= \left(\frac{1}{\sqrt{2\pi\sigma^2}}\right)^N \exp\left(-\frac{1}{2\sigma^2}\sum_{n=1}^{N}(x_n – \mu)^2\right) \end{align} $$
事前分布
$\mu$ の事前分布として正規分布を選びます。
$$ \begin{equation} p(\mu) = \mathcal{N}(\mu|\mu_0, \sigma_0^2) = \frac{1}{\sqrt{2\pi\sigma_0^2}} \exp\left(-\frac{(\mu – \mu_0)^2}{2\sigma_0^2}\right) \end{equation} $$
ここで $\mu_0$ は事前分布の平均、$\sigma_0^2$ は事前分布の分散です。
事後分布の導出
ベイズの定理より、
$$ p(\mu|\mathcal{D}) \propto p(\mathcal{D}|\mu) \cdot p(\mu) $$
指数部分に注目すると、
$$ \begin{align} \ln p(\mu|\mathcal{D}) &\propto -\frac{1}{2\sigma^2}\sum_{n=1}^{N}(x_n – \mu)^2 – \frac{(\mu – \mu_0)^2}{2\sigma_0^2} \end{align} $$
$\mu$ の2次式として整理します。$\sum_{n=1}^{N}(x_n – \mu)^2$ を展開すると、
$$ \sum_{n=1}^{N}(x_n – \mu)^2 = \sum_{n=1}^{N} x_n^2 – 2\mu \sum_{n=1}^{N} x_n + N\mu^2 $$
$\mu$ に関する項のみを取り出すと、
$$ \begin{align} -\frac{N\mu^2 – 2\mu N\bar{x}}{2\sigma^2} – \frac{\mu^2 – 2\mu\mu_0}{2\sigma_0^2} \end{align} $$
ここで $\bar{x} = \frac{1}{N}\sum_{n=1}^{N} x_n$ はデータの標本平均です。
$\mu^2$ の係数と $\mu$ の係数を整理すると、
$$ -\frac{1}{2}\left(\frac{N}{\sigma^2} + \frac{1}{\sigma_0^2}\right)\mu^2 + \left(\frac{N\bar{x}}{\sigma^2} + \frac{\mu_0}{\sigma_0^2}\right)\mu $$
これを正規分布の標準形 $-\frac{(\mu – \mu_N)^2}{2\sigma_N^2}$ と比較すると、事後分布は次のようになります。
$$ \begin{equation} p(\mu|\mathcal{D}) = \mathcal{N}(\mu|\mu_N, \sigma_N^2) \end{equation} $$
ここで、事後分布のパラメータは次のとおりです。
$$ \begin{align} \frac{1}{\sigma_N^2} &= \frac{N}{\sigma^2} + \frac{1}{\sigma_0^2} \\ \mu_N &= \sigma_N^2 \left(\frac{N\bar{x}}{\sigma^2} + \frac{\mu_0}{\sigma_0^2}\right) \end{align} $$
精度の解釈
精度（分散の逆数）$\lambda = 1/\sigma^2$ を使って表現すると、
$$ \lambda_N = N\lambda + \lambda_0 $$
事後精度 = データからの精度 + 事前精度、という加法的な関係が成り立ちます。
また、事後平均 $\mu_N$ は、データの標本平均 $\bar{x}$ と事前平均 $\mu_0$ の精度で重み付けした加重平均と解釈できます。
$$ \mu_N = \frac{N\lambda}{N\lambda + \lambda_0}\bar{x} + \frac{\lambda_0}{N\lambda + \lambda_0}\mu_0 $$
データが増えるほどデータの標本平均の重みが大きくなり、事前分布の影響が薄れていきます。
予測分布
新しいデータ点 $x_*$ の予測分布は、事後分布を使って周辺化することで得られます。
$$ \begin{align} p(x_*|\mathcal{D}) &= \int p(x_*|\mu) \cdot p(\mu|\mathcal{D}) \, d\mu \\ &= \mathcal{N}(x_*|\mu_N, \sigma^2 + \sigma_N^2) \end{align} $$
予測分布の分散 $\sigma^2 + \sigma_N^2$ は、データのノイズ $\sigma^2$ とパラメータの不確かさ $\sigma_N^2$ の和になっています。
Python での実装
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
# 真のパラメータ
true_mu = 3.0
true_sigma = 1.5
sigma_sq = true_sigma**2 # 既知とする
# 事前分布のパラメータ
mu_0 = 0.0
sigma_0_sq = 4.0
# データ生成
np.random.seed(42)
N_total = 100
data = np.random.normal(true_mu, true_sigma, N_total)
# 逐次ベイズ更新
mu_range = np.linspace(-4, 8, 500)
fig, axes = plt.subplots(2, 3, figsize=(15, 9))
axes = axes.flatten()
observations = [0, 1, 3, 10, 30, 100]
for idx, N in enumerate(observations):
ax = axes[idx]
obs = data[:N]
if N == 0:
mu_N = mu_0
sigma_N_sq = sigma_0_sq
else:
x_bar = obs.mean()
sigma_N_sq = 1.0 / (N / sigma_sq + 1.0 / sigma_0_sq)
mu_N = sigma_N_sq * (N * x_bar / sigma_sq + mu_0 / sigma_0_sq)
# 事後分布 p(mu|D)
posterior = norm.pdf(mu_range, mu_N, np.sqrt(sigma_N_sq))
ax.plot(mu_range, posterior, 'b-', linewidth=2, label=f'Posterior')
ax.fill_between(mu_range, posterior, alpha=0.15, color='blue')
ax.axvline(x=true_mu, color='r', linestyle='--', alpha=0.7, label=f'True mu={true_mu}')
ax.axvline(x=mu_N, color='b', linestyle=':', alpha=0.7, label=f'mu_N={mu_N:.2f}')
ax.set_title(f'N={N}, mu_N={mu_N:.2f}, sigma_N={np.sqrt(sigma_N_sq):.2f}', fontsize=11)
ax.set_xlabel('mu', fontsize=10)
ax.set_ylabel('p(mu|D)', fontsize=10)
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
plt.suptitle('Bayesian Estimation of Normal Distribution Mean', fontsize=14)
plt.tight_layout()
plt.show()
データ数が増えるにつれて、事後分布が真の平均 $\mu = 3.0$ 付近に収束し、分散が小さくなっていく様子が確認できます。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
# 予測分布の可視化
np.random.seed(42)
true_mu = 3.0
true_sigma = 1.5
sigma_sq = true_sigma**2
mu_0 = 0.0
sigma_0_sq = 4.0
data = np.random.normal(true_mu, true_sigma, 100)
x_range = np.linspace(-5, 10, 500)
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
ns = [5, 20, 100]
for idx, N in enumerate(ns):
ax = axes[idx]
obs = data[:N]
x_bar = obs.mean()
sigma_N_sq = 1.0 / (N / sigma_sq + 1.0 / sigma_0_sq)
mu_N = sigma_N_sq * (N * x_bar / sigma_sq + mu_0 / sigma_0_sq)
# 予測分布
pred_var = sigma_sq + sigma_N_sq
pred_pdf = norm.pdf(x_range, mu_N, np.sqrt(pred_var))
ax.plot(x_range, pred_pdf, 'b-', linewidth=2, label='Predictive dist.')
# 真の分布
true_pdf = norm.pdf(x_range, true_mu, true_sigma)
ax.plot(x_range, true_pdf, 'r--', linewidth=1.5, label='True dist.')
ax.hist(obs, bins=15, density=True, alpha=0.3, color='green', label='Data')
ax.set_title(f'N={N}', fontsize=12)
ax.set_xlabel('x', fontsize=10)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
plt.suptitle('Predictive Distribution', fontsize=14)
plt.tight_layout()
plt.show()
まとめ
本記事では、正規分布のパラメータをベイズ推定で学習する方法について解説しました。
- 分散既知の場合、平均 $\mu$ の事前分布に正規分布を使うと、事後分布も正規分布になる（共役事前分布）
- 事後精度は事前精度とデータからの精度の和であり、事後平均は精度で重み付けした加重平均
- データが増えるにつれて事後分布は真のパラメータに収束し、事前分布の影響は薄れる
- 予測分布の分散はデータのノイズとパラメータの不確かさの和になる
次のステップとして、以下の記事も参考にしてください。