ベイズ推定は、機械学習分野でも広く利用されている手法であり、近年も多くの論文等が登場しています。
参考書等では、コイン投げの具体例を通してベルヌーイ分布のベイズ推定を行う例題が頻出しますが、今回取り上げる多変量ガウス分布のパラメータをベイズ推定する問題もよく登場します。
今回は共分散行列 $\bm{\Sigma}$ が既知の場合に、平均ベクトル $\bm{\mu}$ をベイズ推定します。
本記事の内容
- 問題設定と仮定
- 事後分布の解析的導出
- 精度行列を使った表現
- Python での実装と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
問題設定
$D$ 次元のデータ $\bm{x}_1, \bm{x}_2, \dots, \bm{x}_N$ が多変量ガウス分布に従うとします。
$$ \begin{equation} \bm{x}_n \sim \mathcal{N}(\bm{x}_n|\bm{\mu}, \bm{\Sigma}) \end{equation} $$
共分散行列 $\bm{\Sigma}$ は既知とし、平均ベクトル $\bm{\mu}$ をベイズ推定します。
尤度関数
データ $\mathcal{D} = \{\bm{x}_1, \dots, \bm{x}_N\}$ が独立に生成されたとき、
$$ \begin{align} p(\mathcal{D}|\bm{\mu}) &= \prod_{n=1}^{N} \mathcal{N}(\bm{x}_n|\bm{\mu}, \bm{\Sigma}) \\ &\propto \exp\left(-\frac{1}{2}\sum_{n=1}^{N}(\bm{x}_n – \bm{\mu})^T \bm{\Sigma}^{-1} (\bm{x}_n – \bm{\mu})\right) \end{align} $$
$\bm{\mu}$ に関する項を整理すると、
$$ \begin{align} \sum_{n=1}^{N}(\bm{x}_n – \bm{\mu})^T \bm{\Sigma}^{-1} (\bm{x}_n – \bm{\mu}) &= N(\bm{\mu} – \bar{\bm{x}})^T \bm{\Sigma}^{-1} (\bm{\mu} – \bar{\bm{x}}) + \text{const} \end{align} $$
ここで $\bar{\bm{x}} = \frac{1}{N}\sum_{n=1}^{N}\bm{x}_n$ は標本平均ベクトルです。
事前分布
$\bm{\mu}$ の事前分布として多変量ガウス分布を選びます。
$$ \begin{equation} p(\bm{\mu}) = \mathcal{N}(\bm{\mu}|\bm{\mu}_0, \bm{\Sigma}_0) \end{equation} $$
事後分布の導出
ベイズの定理より、
$$ p(\bm{\mu}|\mathcal{D}) \propto p(\mathcal{D}|\bm{\mu}) \cdot p(\bm{\mu}) $$
指数部分を $\bm{\mu}$ について整理すると（1次元の場合と同様に二次形式を完成させる）、
$$ \begin{equation} p(\bm{\mu}|\mathcal{D}) = \mathcal{N}(\bm{\mu}|\bm{\mu}_N, \bm{\Sigma}_N) \end{equation} $$
精度行列（共分散行列の逆行列）$\bm{\Lambda} = \bm{\Sigma}^{-1}$、$\bm{\Lambda}_0 = \bm{\Sigma}_0^{-1}$ を使って表現すると、
$$ \begin{align} \bm{\Sigma}_N^{-1} &= N\bm{\Lambda} + \bm{\Lambda}_0 \\ \bm{\mu}_N &= \bm{\Sigma}_N (N\bm{\Lambda}\bar{\bm{x}} + \bm{\Lambda}_0\bm{\mu}_0) \end{align} $$
解釈
(7) 式: 事後精度 = データからの精度 + 事前精度（1次元の場合の自然な拡張）
(8) 式: 事後平均は、標本平均 $\bar{\bm{x}}$ と事前平均 $\bm{\mu}_0$ の精度による重み付き平均
データ数 $N$ が増えるにつれて、$N\bm{\Lambda}$ の項が支配的になり、事後分布はデータの標本平均の周りに集中していきます。
Python での実装
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
from matplotlib.patches import Ellipse
def plot_confidence_ellipse(mu, Sigma, ax, n_std=2, **kwargs):
"""共分散行列からn_std標準偏差の楕円を描画"""
eigenvalues, eigenvectors = np.linalg.eigh(Sigma)
angle = np.degrees(np.arctan2(eigenvectors[1, 1], eigenvectors[0, 1]))
ell = Ellipse(xy=mu, width=2*n_std*np.sqrt(eigenvalues[1]),
height=2*n_std*np.sqrt(eigenvalues[0]),
angle=angle, **kwargs)
ax.add_patch(ell)
# 真のパラメータ
true_mu = np.array([3.0, 2.0])
Sigma_known = np.array([[1.0, 0.3],
[0.3, 0.5]])
Lambda = np.linalg.inv(Sigma_known)
# 事前分布のパラメータ
mu_0 = np.array([0.0, 0.0])
Sigma_0 = np.array([[4.0, 0.0],
[0.0, 4.0]])
Lambda_0 = np.linalg.inv(Sigma_0)
# データ生成
np.random.seed(42)
N_total = 100
data = np.random.multivariate_normal(true_mu, Sigma_known, N_total)
fig, axes = plt.subplots(2, 3, figsize=(16, 10))
axes = axes.flatten()
observations = [0, 1, 5, 15, 50, 100]
for idx, N in enumerate(observations):
ax = axes[idx]
if N == 0:
mu_N = mu_0
Sigma_N = Sigma_0
else:
x_bar = data[:N].mean(axis=0)
Sigma_N_inv = N * Lambda + Lambda_0
Sigma_N = np.linalg.inv(Sigma_N_inv)
mu_N = Sigma_N @ (N * Lambda @ x_bar + Lambda_0 @ mu_0)
# データ点
if N > 0:
ax.scatter(data[:N, 0], data[:N, 1], s=15, alpha=0.4, color='gray', label='Data')
# 真の平均
ax.scatter(*true_mu, s=150, marker='*', color='red', zorder=5, label='True mu')
# 事後分布の平均
ax.scatter(*mu_N, s=100, marker='D', color='blue', zorder=5, label='Posterior mu')
# 事後分布の楕円
for n_std in [1, 2]:
plot_confidence_ellipse(mu_N, Sigma_N, ax, n_std=n_std,
fill=False, color='blue', linewidth=1.5, linestyle='--')
ax.set_xlabel('x1', fontsize=10)
ax.set_ylabel('x2', fontsize=10)
ax.set_title(f'N = {N}', fontsize=12)
ax.set_xlim(-3, 6)
ax.set_ylim(-3, 5)
ax.legend(fontsize=8, loc='upper left')
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
plt.suptitle('Bayesian Estimation of Multivariate Gaussian Mean', fontsize=14)
plt.tight_layout()
plt.show()
データ数が増えるにつれて、事後分布の楕円が小さくなり（不確かさが減少）、事後平均が真のパラメータに収束していく様子が確認できます。
まとめ
本記事では、多変量ガウス分布の平均ベクトルのベイズ推定について解説しました。
- 共分散行列既知の場合、平均ベクトルの共役事前分布は多変量ガウス分布
- 事後分布も多変量ガウス分布になり、パラメータは精度行列を使って簡潔に表現できる
- 事後精度はデータからの精度と事前精度の和であり、1次元の場合の自然な拡張
- データ数が増えるほど事後分布は真のパラメータの周りに集中する
次のステップとして、以下の記事も参考にしてください。