はじめに

「コインを10回投げて7回表が出た。このコインは偏っているか？」

頻度論（p値）では「帰無仮説のもとで…」という回りくどい手順が必要です。ベイズ統計なら「表が出る確率の分布」を直接得られます。

本記事では、ベイズ統計の核心的な数式とPythonによる実装を丁寧に解説します。

ベイズの定理：いつ使うか

いつ使うか

過去の知識（事前分布）とデータを組み合わせて推論したいとき
サンプルが少なく、頻度論の推定が不安定なとき
「〜の確率は何%か」という形で答えを直接出したいとき
A/Bテスト・パラメータ推定・分類モデルの確率校正など

数式

ベイズの定理：

$$P(\theta \mid \mathcal{D}) = \frac{P(\mathcal{D} \mid \theta) \cdot P(\theta)}{P(\mathcal{D})}$$

各項の意味：

項	名称	意味
$P(\theta \mid \mathcal{D})$	事後分布（Posterior）	データを見た後のパラメータの確率分布
$P(\mathcal{D} \mid \theta)$	尤度（Likelihood）	パラメータ $\theta$ のもとでデータが得られる確率
$P(\theta)$	事前分布（Prior）	データを見る前のパラメータの信念
$P(\mathcal{D})$	周辺尤度（Evidence）	正規化定数（定数なので比例式で扱うことが多い）

比例式で書くと：

$$P(\theta \mid \mathcal{D}) \propto P(\mathcal{D} \mid \theta) \cdot P(\theta)$$

事後分布 ∝ 尤度 × 事前分布

Beta分布と共役事前分布

いつ使うか

$\theta \in [0, 1]$（確率・割合）をパラメータとして推定するとき
二項分布の尤度と組み合わせると解析的に更新できる（共役事前分布）

数式

Beta分布の確率密度関数（PDF）：

$$f(\theta; \alpha, \beta) = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \theta^{\alpha-1}(1-\theta)^{\beta-1}, \quad \theta \in [0,1]$$

期待値と分散：

$$E[\theta] = \frac{\alpha}{\alpha + \beta}, \quad \text{Var}[\theta] = \frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}$$

二項尤度 × Beta事前分布 → Beta事後分布（共役性）：

データ：$n$ 回試行で $k$ 回成功

$$P(\mathcal{D} \mid \theta) = \binom{n}{k}\theta^k(1-\theta)^{n-k}$$

事前分布：$P(\theta) = \text{Beta}(\alpha_0, \beta_0)$

事後分布：

$$P(\theta \mid \mathcal{D}) = \text{Beta}(\alpha_0 + k,; \beta_0 + n - k)$$

→ 成功回数を $\alpha$ に、失敗回数を $\beta$ に足すだけで更新できる。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# コインを10回投げて7回表
n_trials, n_heads = 10, 7

# 事前分布: Beta(1, 1) = 一様分布（知識なし）
prior_a, prior_b = 1, 1

# ベイズ更新：足し算だけ
post_a = prior_a + n_heads
post_b = prior_b + (n_trials - n_heads)

prior_dist = stats.beta(prior_a, prior_b)
post_dist  = stats.beta(post_a, post_b)

print(f"事後分布: Beta({post_a}, {post_b})")
print(f"事後期待値 E[θ] = {post_a/(post_a+post_b):.3f}")
print(f"事後分散 Var[θ] = {post_a*post_b/((post_a+post_b)**2*(post_a+post_b+1)):.4f}")

theta = np.linspace(0, 1, 300)
plt.figure(figsize=(9, 4))
plt.plot(theta, prior_dist.pdf(theta), "--", color="gray", label="事前 Beta(1,1)")
plt.plot(theta, post_dist.pdf(theta), color="steelblue", lw=2,
         label=f"事後 Beta({post_a},{post_b})")
plt.axvline(0.5, color="red", linestyle=":", label="θ=0.5（公平）")
plt.xlabel("θ（表が出る確率）")
plt.ylabel("確率密度")
plt.legend(); plt.tight_layout(); plt.show()

信用区間（Credible Interval）

いつ使うか

「θ が X〜Y の範囲にある確率は 95%」と直接言いたいとき
頻度論の信頼区間（「この手順を繰り返すと 95% のケースで区間が真値を含む」）より解釈が直感的

数式

95% 等裾信用区間（Equal-tail Credible Interval）：

$$\text{CI}{0.95} = \left[\theta{0.025},; \theta_{0.975}\right]$$

$\theta_{q}$ は事後分布の $q$ 分位点。Beta分布なら：

$$\theta_q = F^{-1}_{\text{Beta}(\alpha, \beta)}(q)$$

lower = post_dist.ppf(0.025)
upper = post_dist.ppf(0.975)
print(f"95%信用区間: [{lower:.3f}, {upper:.3f}]")

# 「θ > 0.5 の確率」を直接計算
prob_biased = 1 - post_dist.cdf(0.5)
print(f"表が出やすい（θ>0.5）確率: {prob_biased:.1%}")

ベイズA/Bテスト

いつ使うか

「BがAより優れている確率は何%か」をビジネス判断に使いたいとき
サンプルが少なく途中経過でも意思決定が必要なとき

数式

A群・B群それぞれの事後分布：

$$\theta_A \mid \mathcal{D}_A \sim \text{Beta}(1+x_A,; 1+n_A-x_A)$$
$$\theta_B \mid \mathcal{D}_B \sim \text{Beta}(1+x_B,; 1+n_B-x_B)$$

「BがAより優れている確率」：

$$P(\theta_B > \theta_A \mid \mathcal{D}) = \int_0^1 P(\theta_B > \theta_A \mid \theta_A) , p(\theta_A), d\theta_A$$

解析解もあるが、モンテカルロ積分が実用的：

$$\hat{P}(\theta_B > \theta_A) \approx \frac{1}{S}\sum_{s=1}^{S} \mathbf{1}[\theta_B^{(s)} > \theta_A^{(s)}]$$

n_a, x_a = 1000, 50
n_b, x_b = 1000, 65

post_a = stats.beta(1 + x_a, 1 + n_a - x_a)
post_b = stats.beta(1 + x_b, 1 + n_b - x_b)

# モンテカルロ積分
S = 200_000
samp_a = post_a.rvs(S)
samp_b = post_b.rvs(S)

prob_b_better = np.mean(samp_b > samp_a)
expected_lift = np.mean((samp_b - samp_a) / samp_a)

print(f"E[θ_A] = {post_a.mean():.4f}")
print(f"E[θ_B] = {post_b.mean():.4f}")
print(f"P(θ_B > θ_A | D) = {prob_b_better:.1%}")
print(f"期待改善率: {expected_lift:.1%}")

データが増えると事後分布はどう変わるか

事前知識の影響は $n \to \infty$ で消え、事後平均は最尤推定量（MLE）に収束する。

$$E[\theta \mid \mathcal{D}] = \frac{\alpha_0 + k}{\alpha_0 + \beta_0 + n} \xrightarrow{n\to\infty} \frac{k}{n} = \hat{\theta}_{\text{MLE}}$$

experiments = [(0,0,"データなし"), (3,2,"3回中2表"), (10,7,"10回中7表"), (100,70,"100回中70表")]
fig, axes = plt.subplots(1, 4, figsize=(16, 4))
for ax, (n, k, label) in zip(axes, experiments):
    dist = stats.beta(1+k, 1+(n-k))
    theta = np.linspace(0, 1, 300)
    ax.plot(theta, dist.pdf(theta), color="steelblue", lw=2)
    ax.axvline(0.7, color="red", ls="--", alpha=0.5, label="真値")
    ax.set_title(label); ax.set_xlabel("θ"); ax.set_xlim(0, 1)
plt.suptitle("データが増えると事後分布が真値（0.7）に収束する", y=1.02)
plt.tight_layout(); plt.show()

まとめ

概念	数式	意味
ベイズの定理	$P(\theta\|D) \propto P(D\|\theta)P(\theta)$	事後 ∝ 尤度 × 事前
Beta分布の期待値	$E[\theta]=\alpha/(\alpha+\beta)$	更新後の割合の推定値
ベイズ更新（Beta-Binomial）	$\text{Beta}(\alpha_0+k,;\beta_0+n-k)$	成功・失敗を足すだけ
信用区間	$[F^{-1}(0.025), F^{-1}(0.975)]$	事後分布の分位点