第4章：ベイズ統計入門

イントロダクション

前章までで、記述統計と確率の基礎を学びました。本章では、その応用として近年ますます重要性を増しているベイズ統計（Bayesian Statistics）を扱います。ベイズ統計は、観測データが得られるたびに知識を更新していくという、人間の学習プロセスに近い考え方を数学的に定式化した枠組みです。

これまで前提としてきた頻度論統計（Frequentist Statistics）では、母集団のパラメータ（例えばコインの表が出る確率）は未知だが固定された値であるとみなし、観測データからその値を推定します。一方、ベイズ統計では、パラメータ自体を確率分布として扱い、データを観測するたびにその分布を更新していきます。この違いは単なる数学的な形式の違いにとどまらず、不確実性の捉え方そのものに関わる重要な視点の転換です。

1. ベイズの定理の深い理解

第1章では、ベイズの定理を事象（陽性/陰性、病気/健康など）の間の確率関係として学びました。本章では、これをパラメータ推定という視点から捉え直します。

データ$D$が与えられたときに、未知のパラメータ$\theta$についての信念を更新するというのが、ベイズ推定の基本的な考え方です。ベイズの定理をパラメータ推定の文脈で書くと、次のようになります。

$$P(\theta|D) = \frac{P(D|\theta)P(\theta)}{P(D)}$$

各項の名称と役割は次の通りです。

• $P(\theta)$: 事前分布（Prior Distribution）— データを観測する前に持っていた、パラメータ$\theta$についての信念
• $P(D|\theta)$: 尤度（Likelihood）— パラメータが$\theta$であるという仮定のもとで、実際に観測されたデータ$D$が得られる確率（密度）
• $P(\theta|D)$: 事後分布（Posterior Distribution）— データ$D$を観測した後に更新された、パラメータ$\theta$についての信念
• $P(D)$: 周辺尤度（Marginal Likelihood）またはエビデンス（Evidence）— あらゆる$\theta$の値にわたってデータ$D$が観測される確率（密度）。正規化定数として働く

$\theta$が連続値をとる場合、周辺尤度は次の積分で与えられます。

$$P(D) = \int P(D|\theta)P(\theta)\,d\theta$$

この積分は多くの場合解析的に解けないため、ベイズ統計の実践では大きな課題となります。本章の後半で扱うMCMCは、この課題を回避するための手法の一つです。

2. 事前分布、尤度、事後分布の関係

ベイズ推定の全体像を整理すると、次の3つの要素の相互作用として理解できます。

2.1 事前分布：観測前の知識

事前分布はデータを見る前の知識や仮定を表現します。過去の研究、専門家の意見、あるいは「何も分からない」という状態さえも、事前分布として定式化できます。

• 無情報事前分布（Uninformative Prior）: どの値も同程度にありうるとする、事前知識をほとんど入れない事前分布
• 弱情報事前分布（Weakly Informative Prior）: おおまかな範囲だけを制約する事前分布
• 情報事前分布（Informative Prior）: 過去のデータや専門知識に基づく、強い制約を持つ事前分布

2.2 尤度：データが語ること

尤度関数はデータの生成過程をモデル化したものです。例えば、コインを$n$回投げて$k$回表が出たという観測は、パラメータ$\theta$（表が出る確率）を使って二項分布でモデル化できます。

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

2.3 事後分布：更新された知識

事後分布は事前分布と尤度を掛け合わせ（て正規化し）たものです。データが多くなるほど尤度の影響が相対的に大きくなり、事前分布の影響は薄れていきます。逆にデータが少ないうちは、事前分布が事後分布の形を強く左右します。

3. 共役事前分布の利用

ベイズ推定の計算を解析的に行うために便利な概念が共役事前分布（Conjugate Prior）です。

代表的な共役事前分布の組み合わせを以下に示します。

3.1 ベータ分布とベルヌーイ試行の共役性

ベータ分布（Beta Distribution）は、$[0,1]$区間の値をとる確率変数のための確率分布で、確率$\theta$そのものの不確実性を表現するのに適しています。確率密度関数は次の通りです。

$$f(\theta;\alpha,\beta) = \frac{\theta^{\alpha-1}(1-\theta)^{\beta-1}}{B(\alpha,\beta)}$$

ここで$B(\alpha,\beta)$はベータ関数で、正規化定数として働きます。

事前分布をBeta($\alpha$, $\beta$)とし、$n$回のベルヌーイ試行（Bernoulli Trial、成功/失敗の2値をとる試行）で$k$回成功したというデータを観測すると、事後分布は次のように解析的に求まります。

$$\theta \sim \text{Beta}(\alpha, \beta) \quad \Rightarrow \quad \theta|D \sim \text{Beta}(\alpha+k, \beta+n-k)$$

この更新式が示す通り、パラメータ$\alpha$は「成功の仮想回数」、$\beta$は「失敗の仮想回数」として解釈でき、観測データの成功・失敗回数をそれぞれに単純に加算するだけで事後分布が得られます。

3.2 Pythonでの実装：ベータ-ベルヌーイ更新

コインを繰り返し投げ、観測データが増えるにつれて事後分布がどのように変化するかを実装してみましょう。

import numpy as np
from scipy import stats

# 真の確率(未知としてシミュレーション用に設定)でコイン投げを生成
np.random.seed(42)
true_p = 0.7
n_flips = 50
flips = np.random.binomial(1, true_p, size=n_flips)  # 1=表, 0=裏

# 事前分布 Beta(2, 2): やや公平寄りだが弱い事前知識
alpha_prior, beta_prior = 2, 2

# 観測数を増やしながら事後分布のパラメータを更新
checkpoints = [0, 5, 20, 50]

print("=== ベータ-ベルヌーイ更新の推移 ===")
for n in checkpoints:
    successes = flips[:n].sum()
    failures = n - successes
    alpha_post = alpha_prior + successes
    beta_post = beta_prior + failures
    posterior_mean = alpha_post / (alpha_post + beta_post)
    print(f"n={n:2d}回観測後: 成功{successes}回, "
          f"Beta({alpha_post}, {beta_post}), 事後平均={posterior_mean:.3f}")

実行結果:

=== ベータ-ベルヌーイ更新の推移 ===
n= 0回観測後: 成功0回, Beta(2, 2), 事後平均=0.500
n= 5回観測後: 成功3回, Beta(5, 4), 事後平均=0.556
n=20回観測後: 成功14回, Beta(16, 8), 事後平均=0.667
n=50回観測後: 成功39回, Beta(41, 13), 事後平均=0.759

事前分布から観測データが増えるにつれて事後分布がどのように鋭くなっていくかを可視化してみましょう。

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

np.random.seed(42)
true_p = 0.7
n_flips = 50
flips = np.random.binomial(1, true_p, size=n_flips)

alpha_prior, beta_prior = 2, 2
checkpoints = [0, 5, 20, 50]
theta = np.linspace(0.001, 0.999, 200)

plt.figure(figsize=(10, 6))
for n in checkpoints:
    successes = flips[:n].sum()
    failures = n - successes
    alpha_post = alpha_prior + successes
    beta_post = beta_prior + failures
    pdf = stats.beta.pdf(theta, alpha_post, beta_post)
    plt.plot(theta, pdf, label=f'n={n} (成功{successes}回)', linewidth=2)

plt.axvline(true_p, color='black', linestyle=':', label=f'真の確率={true_p}')
plt.xlabel('θ（表が出る確率）', fontsize=12)
plt.ylabel('確率密度', fontsize=12)
plt.title('観測データの増加に伴う事後分布の変化', fontsize=14)
plt.legend()
plt.grid(alpha=0.3)
plt.show()

4. ベイズ推定の実装

事後分布が得られたら、そこから具体的な推定値や不確実性の指標を取り出す必要があります。

4.1 点推定：事後平均とMAP推定

• 事後平均（Posterior Mean）: 事後分布の期待値$E[\theta|D]$。最も一般的な点推定量
• MAP推定（Maximum A Posteriori estimation）: 事後分布の密度が最大となる$\theta$の値、すなわち事後分布の最頻値

Beta($\alpha$, $\beta$)分布については、それぞれ次の式で解析的に求められます。

$$E[\theta|D] = \frac{\alpha}{\alpha+\beta} \qquad \theta_{\text{MAP}} = \frac{\alpha-1}{\alpha+\beta-2}$$

（ただしMAP推定の式は$\alpha > 1$かつ$\beta > 1$のときに成り立ちます）

4.2 区間推定：信用区間

信用区間（Credible Interval）は、事後分布において特定の確率質量（例えば95%）を含む区間です。頻度論統計の信頼区間（Confidence Interval）とよく似た役割を果たしますが、解釈は根本的に異なります（この違いは第6節で詳しく扱います）。

4.3 Pythonでの実装：点推定と区間推定

import numpy as np
from scipy import stats

def posterior_summary(alpha_post, beta_post, cred_mass=0.95):
    """
    ベータ事後分布から点推定と信用区間を計算する

    Parameters:
    -----------
    alpha_post, beta_post : float
        事後分布 Beta(alpha_post, beta_post) のパラメータ
    cred_mass : float
        信用区間の確率質量（デフォルト95%）

    Returns:
    --------
    mean, map_estimate, (lower, upper)
    """
    mean = alpha_post / (alpha_post + beta_post)
    if alpha_post > 1 and beta_post > 1:
        map_estimate = (alpha_post - 1) / (alpha_post + beta_post - 2)
    else:
        map_estimate = None
    lower = stats.beta.ppf((1 - cred_mass) / 2, alpha_post, beta_post)
    upper = stats.beta.ppf(1 - (1 - cred_mass) / 2, alpha_post, beta_post)
    return mean, map_estimate, (lower, upper)

# 前節のコイン投げデータ(50回中39回成功、事前分布 Beta(2,2))
alpha_post, beta_post = 41, 13
mean, map_est, (lower, upper) = posterior_summary(alpha_post, beta_post)

print(f"事後平均: {mean:.4f}")
print(f"MAP推定値: {map_est:.4f}")
print(f"95%信用区間: [{lower:.4f}, {upper:.4f}]")

実行結果:

事後平均: 0.7593
MAP推定値: 0.7692
95%信用区間: [0.6379, 0.8624]

4.4 実践例：A/Bテストのベイズ的比較

ウェブサイトの2つのデザイン案（A案とB案）のクリック率を比較するA/Bテストを例に、ベイズ推定を実践してみましょう。

import numpy as np
from scipy import stats

# A案: 500人中48人がクリック、B案: 480人中63人がクリック
alpha_prior, beta_prior = 1, 1  # 無情報事前分布 Beta(1,1) = 一様分布

n_a, success_a = 500, 48
n_b, success_b = 480, 63

alpha_a = alpha_prior + success_a
beta_a = beta_prior + (n_a - success_a)
alpha_b = alpha_prior + success_b
beta_b = beta_prior + (n_b - success_b)

def posterior_summary(alpha_post, beta_post, cred_mass=0.95):
    mean = alpha_post / (alpha_post + beta_post)
    map_estimate = (alpha_post - 1) / (alpha_post + beta_post - 2)
    lower = stats.beta.ppf((1 - cred_mass) / 2, alpha_post, beta_post)
    upper = stats.beta.ppf(1 - (1 - cred_mass) / 2, alpha_post, beta_post)
    return mean, map_estimate, (lower, upper)

mean_a, map_a, ci_a = posterior_summary(alpha_a, beta_a)
mean_b, map_b, ci_b = posterior_summary(alpha_b, beta_b)

print(f"A案: 事後平均={mean_a:.4f}, 95%信用区間=[{ci_a[0]:.4f}, {ci_a[1]:.4f}]")
print(f"B案: 事後平均={mean_b:.4f}, 95%信用区間=[{ci_b[0]:.4f}, {ci_b[1]:.4f}]")

# モンテカルロ法でP(B案の方が優れている)を推定
np.random.seed(0)
n_mc = 200000
samples_a = np.random.beta(alpha_a, beta_a, n_mc)
samples_b = np.random.beta(alpha_b, beta_b, n_mc)
prob_b_better = np.mean(samples_b > samples_a)

print(f"\nP(B案のクリック率 > A案のクリック率) = {prob_b_better:.4f}")

実行結果:

A案: 事後平均=0.0976, 95%信用区間=[0.0732, 0.1250]
B案: 事後平均=0.1328, 95%信用区間=[0.1040, 0.1644]

P(B案のクリック率 > A案のクリック率) = 0.9582

5. マルコフ連鎖モンテカルロ法（MCMC）入門

ベータ分布とベルヌーイ試行のように共役事前分布が使える場合は、事後分布を解析的に求められます。しかし現実の多くのベイズモデルでは、尤度と事前分布の組み合わせが共役にならず、事後分布の式を閉じた形で書き下すことができません。このような場合に威力を発揮するのがマルコフ連鎖モンテカルロ法（Markov Chain Monte Carlo, MCMC）です。

MCMCは、事後分布の正規化定数$P(D)$を計算せずに、事後分布に従うサンプル列を生成する手法群の総称です。得られたサンプルの集合を使えば、事後平均や信用区間などの統計量をヒストグラムや標本平均から近似的に求められます。

5.1 マルコフ連鎖とモンテカルロ法

モンテカルロ法（Monte Carlo Method）は、乱数を用いた試行を繰り返すことで、解析的に求めるのが難しい量（積分値や期待値など）を近似的に計算する手法の総称です。

マルコフ連鎖（Markov Chain）は、次の状態が現在の状態のみに依存し、それ以前の履歴には依存しないという性質（マルコフ性）を持つ確率過程です。MCMCは、目的の事後分布を定常分布（Stationary Distribution）として持つようなマルコフ連鎖を構築し、それに沿ってランダムウォークすることでサンプルを生成します。

5.2 メトロポリス法

MCMCの中でも直感的に理解しやすいアルゴリズムがメトロポリス法（Metropolis Algorithm）です。手順は次の通りです。

1. 初期値$\theta_0$を適当に設定する
2. 現在の状態$\theta_t$から、提案分布（Proposal Distribution、多くの場合は正規分布などの対称な分布）を使って新しい候補$\theta^*$を生成する
3. 受理確率（Acceptance Probability）$a = \min\left(1, \dfrac{P(\theta^*|D)}{P(\theta_t|D)}\right)$を計算する
4. 確率$a$で候補を受理し$\theta_{t+1}=\theta^*$とする。そうでなければ棄却し$\theta_{t+1}=\theta_t$とする
5. 手順2〜4を必要な回数だけ繰り返す

生成されたサンプル列の初期部分は、初期値$\theta_0$の影響を強く受けており、まだ定常分布に収束していない可能性があります。この初期部分をバーンイン（Burn-in）期間として捨て、残りのサンプルだけを事後分布の近似として使うのが一般的です。

5.3 Pythonでの実装：手書きメトロポリスサンプラー

ベータ-ベルヌーイモデルは解析的に事後分布が求まるため、MCMCで得られた結果を解析解と比較検証できます。これはMCMC実装が正しく動作しているかを確認する良い練習になります。

import numpy as np
from scipy import stats

def log_posterior(theta_val, successes, failures, alpha_prior, beta_prior):
    """正規化されていない事後分布の対数を計算(0<theta<1の範囲外は-infとする)"""
    if theta_val <= 0 or theta_val >= 1:
        return -np.inf
    log_lik = successes * np.log(theta_val) + failures * np.log(1 - theta_val)
    log_prior = (alpha_prior - 1) * np.log(theta_val) + (beta_prior - 1) * np.log(1 - theta_val)
    return log_lik + log_prior

def metropolis_sampler(log_post_fn, n_samples, init, proposal_std, seed=42):
    """
    メトロポリス法によるMCMCサンプラー

    Parameters:
    -----------
    log_post_fn : callable
        対数事後分布(正規化不要)を返す関数
    n_samples : int
        生成するサンプル数
    init : float
        初期値
    proposal_std : float
        提案分布(正規分布)の標準偏差

    Returns:
    --------
    samples : ndarray
        生成されたサンプル列
    acceptance_rate : float
        受理率
    """
    rng = np.random.default_rng(seed)
    samples = np.zeros(n_samples)
    current = init
    current_log_p = log_post_fn(current)
    n_accepted = 0

    for i in range(n_samples):
        proposal = current + rng.normal(0, proposal_std)
        proposal_log_p = log_post_fn(proposal)
        log_accept_ratio = proposal_log_p - current_log_p

        if np.log(rng.uniform()) < log_accept_ratio:
            current = proposal
            current_log_p = proposal_log_p
            n_accepted += 1
        samples[i] = current

    acceptance_rate = n_accepted / n_samples
    return samples, acceptance_rate

# 3節と同じコイン投げデータ(50回中39回成功)を使い、Beta(2,2)を事前分布とする
successes_total, failures_total = 39, 11
alpha_prior, beta_prior = 2, 2
target = lambda th: log_posterior(th, successes_total, failures_total, alpha_prior, beta_prior)

n_samples = 20000
samples, acc_rate = metropolis_sampler(target, n_samples, init=0.5, proposal_std=0.1)

burn_in = 2000
post_burn = samples[burn_in:]

analytic_mean = (alpha_prior + successes_total) / (alpha_prior + beta_prior + successes_total + failures_total)

print(f"受理率: {acc_rate:.3f}")
print(f"MCMCによる事後平均: {post_burn.mean():.4f}")
print(f"解析解(Beta分布)による事後平均: {analytic_mean:.4f}")
print(f"MCMCによる事後標準偏差: {post_burn.std():.4f}")

実行結果:

受理率: 0.548
MCMCによる事後平均: 0.7601
解析解(Beta分布)による事後平均: 0.7593
MCMCによる事後標準偏差: 0.0575

MCMCサンプルの挙動を可視化し、解析解と比較してみましょう。

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

def log_posterior(theta_val, successes, failures, alpha_prior, beta_prior):
    if theta_val <= 0 or theta_val >= 1:
        return -np.inf
    log_lik = successes * np.log(theta_val) + failures * np.log(1 - theta_val)
    log_prior = (alpha_prior - 1) * np.log(theta_val) + (beta_prior - 1) * np.log(1 - theta_val)
    return log_lik + log_prior

def metropolis_sampler(log_post_fn, n_samples, init, proposal_std, seed=42):
    rng = np.random.default_rng(seed)
    samples = np.zeros(n_samples)
    current = init
    current_log_p = log_post_fn(current)
    n_accepted = 0
    for i in range(n_samples):
        proposal = current + rng.normal(0, proposal_std)
        proposal_log_p = log_post_fn(proposal)
        if np.log(rng.uniform()) < proposal_log_p - current_log_p:
            current = proposal
            current_log_p = proposal_log_p
            n_accepted += 1
        samples[i] = current
    return samples, n_accepted / n_samples

successes_total, failures_total = 39, 11
alpha_prior, beta_prior = 2, 2
target = lambda th: log_posterior(th, successes_total, failures_total, alpha_prior, beta_prior)
samples, acc_rate = metropolis_sampler(target, 20000, init=0.5, proposal_std=0.1)
burn_in = 2000
post_burn = samples[burn_in:]

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# トレースプロット: サンプルが反復とともにどう推移したか
axes[0].plot(samples, linewidth=0.5, alpha=0.7)
axes[0].axvline(burn_in, color='red', linestyle='--', label='バーンイン境界')
axes[0].set_xlabel('反復回数')
axes[0].set_ylabel('θ')
axes[0].set_title('トレースプロット')
axes[0].legend()

# ヒストグラム vs 解析解
theta_grid = np.linspace(0.001, 0.999, 200)
analytic_pdf = stats.beta.pdf(theta_grid, alpha_prior + successes_total, beta_prior + failures_total)
axes[1].hist(post_burn, bins=50, density=True, alpha=0.6, label='MCMCサンプル')
axes[1].plot(theta_grid, analytic_pdf, 'r-', linewidth=2, label='解析解 Beta(41, 13)')
axes[1].set_xlabel('θ')
axes[1].set_ylabel('確率密度')
axes[1].set_title('MCMCサンプルと解析解の比較')
axes[1].legend()

plt.tight_layout()
plt.show()

6. ベイズ統計と頻度論統計の比較

ここまで学んだベイズ統計と、これまで前提としてきた頻度論統計の違いを整理します。両者は同じデータから異なる哲学で推論を行うため、それぞれの長所と限界を理解しておくことが重要です。

6.1 信頼区間と信用区間の解釈の違い

特に混同されやすいのが、頻度論の信頼区間と、ベイズの信用区間の解釈の違いです。

• 95%信頼区間: 「同じ手続きでデータ収集と区間計算を100回繰り返したとき、そのうち平均して95回は真のパラメータ値を含む区間になる」という、手続き全体についての性質。個々の区間が真の値を含むかどうかは0か1であり、確率的な言明はできない
• 95%信用区間: 「観測されたこのデータのもとで、パラメータがこの区間に入る確率は95%である」という、パラメータそのものについての直接的な確率的言明

6.2 小標本での挙動の違い：実装比較

標本サイズが小さい場合、両者の違いがより顕著になることがあります。10回中8回成功したという少ないデータで、頻度論のWald信頼区間とベイズの信用区間を比較してみましょう。

import numpy as np
from scipy import stats

n = 10
k = 8

# 頻度論: 最尤推定値とWald信頼区間
p_hat = k / n  # 最尤推定量(Maximum Likelihood Estimate, MLE)
se = np.sqrt(p_hat * (1 - p_hat) / n)
z = 1.96
freq_ci = (p_hat - z * se, p_hat + z * se)

# ベイズ: 無情報事前分布 Beta(1,1) を用いた事後分布
alpha_prior, beta_prior = 1, 1
alpha_post = alpha_prior + k
beta_post = beta_prior + (n - k)
bayes_mean = alpha_post / (alpha_post + beta_post)
bayes_ci = stats.beta.ppf([0.025, 0.975], alpha_post, beta_post)

print(f"最尤推定値: {p_hat:.4f}")
print(f"頻度論 95%信頼区間(Wald法): [{freq_ci[0]:.4f}, {freq_ci[1]:.4f}]")
print(f"\nベイズ事後平均: {bayes_mean:.4f}")
print(f"ベイズ 95%信用区間: [{bayes_ci[0]:.4f}, {bayes_ci[1]:.4f}]")

実行結果:

最尤推定値: 0.8000
頻度論 95%信頼区間(Wald法): [0.5521, 1.0479]

ベイズ事後平均: 0.7500
ベイズ 95%信用区間: [0.4822, 0.9398]

この例は、「どちらの手法が絶対的に正しいか」という話ではなく、それぞれの手法が異なる仮定と近似のもとで動いていることを示しています。標本サイズが大きくなり、弱い事前分布を用いる限り、頻度論とベイズの結果は多くの場合近い値に収束していきます。状況に応じてどちらの枠組みが適切かを判断する視点を持つことが重要です。

7. まとめと次のステップ

この章では、ベイズ統計の考え方を、理論からPythonでの実装まで一貫して学びました。

学習目標の振り返り

この章の冒頭で示した学習目標を振り返ってみましょう。

• ✅ ベイズ統計の考え方を理解する — 事前分布・尤度・事後分布の関係として整理した
• ✅ 事前分布・尤度・事後分布の関係を説明できる — 第1節・第2節で数式と具体例を通じて学んだ
• ✅ 共役事前分布を用いたベイズ推定を実装できる — 第3節・第4節でベータ-ベルヌーイモデルを実装した
• ✅ MCMCの基本的な考え方を理解する — 第5節で手書きのメトロポリスサンプラーを実装し、解析解と照合した
• ✅ ベイズ統計と頻度論統計の違いを説明できる — 第6節で信頼区間と信用区間の解釈の違いを具体例で確認した

次のステップ

次章では、本章で扱った推定の枠組みを踏まえ、仮説検定など統計的推論のさらなるトピックへと進みます。ベイズ統計の考え方は、今後機械学習で頻繁に登場するベイズ最適化やベイズニューラルネットワークなど、より高度な手法の土台にもなります。

練習問題

from scipy import stats

alpha_prior, beta_prior = 2, 2
n, k = 10, 7

alpha_post = alpha_prior + k
beta_post = beta_prior + (n - k)

mean = alpha_post / (alpha_post + beta_post)
map_estimate = (alpha_post - 1) / (alpha_post + beta_post - 2)

print(f"事後分布: Beta({alpha_post}, {beta_post})")
print(f"事後平均: {mean:.4f}")
print(f"MAP推定値: {map_estimate:.4f}")

# 出力:
# 事後分布: Beta(9, 5)
# 事後平均: 0.6429
# MAP推定値: 0.6667

import numpy as np
from scipy import stats

n, k = 10, 7
alpha_post, beta_post = 9, 5

# ベイズ信用区間
bayes_ci = stats.beta.ppf([0.025, 0.975], alpha_post, beta_post)

# 頻度論Wald信頼区間
p_hat = k / n
se = np.sqrt(p_hat * (1 - p_hat) / n)
freq_ci = (p_hat - 1.96 * se, p_hat + 1.96 * se)

print(f"ベイズ95%信用区間: [{bayes_ci[0]:.4f}, {bayes_ci[1]:.4f}]")
print(f"頻度論95%信頼区間(Wald法): [{freq_ci[0]:.4f}, {freq_ci[1]:.4f}]")

# 出力:
# ベイズ95%信用区間: [0.3857, 0.8614]
# 頻度論95%信頼区間(Wald法): [0.4160, 0.9840]

import numpy as np

def log_posterior(theta_val, successes, failures, alpha_prior, beta_prior):
    if theta_val <= 0 or theta_val >= 1:
        return -np.inf
    log_lik = successes * np.log(theta_val) + failures * np.log(1 - theta_val)
    log_prior = (alpha_prior - 1) * np.log(theta_val) + (beta_prior - 1) * np.log(1 - theta_val)
    return log_lik + log_prior

def metropolis_sampler(log_post_fn, n_samples, init, proposal_std, seed=42):
    rng = np.random.default_rng(seed)
    samples = np.zeros(n_samples)
    current = init
    current_log_p = log_post_fn(current)
    n_accepted = 0
    for i in range(n_samples):
        proposal = current + rng.normal(0, proposal_std)
        proposal_log_p = log_post_fn(proposal)
        if np.log(rng.uniform()) < proposal_log_p - current_log_p:
            current = proposal
            current_log_p = proposal_log_p
            n_accepted += 1
        samples[i] = current
    return samples, n_accepted / n_samples

target = lambda th: log_posterior(th, successes=7, failures=3, alpha_prior=2, beta_prior=2)

for std in [0.01, 1.0]:
    samples, acc_rate = metropolis_sampler(target, 5000, init=0.5, proposal_std=std)
    print(f"proposal_std={std}: 受理率={acc_rate:.3f}, サンプル標準偏差={samples[500:].std():.4f}")

# 出力:
# proposal_std=0.01: 受理率=0.975, サンプル標準偏差=0.1299
# proposal_std=1.0: 受理率=0.156, サンプル標準偏差=0.1208