「データ解析のための統計モデリング入門」をPythonで写経 Vol.17 ～ ８章「マルコフ連鎖モンテカルロ（MCMC）法とベイズ統計モデル」 ベイズ統計モデリング～メトロポリス法でMCMCを体感

８章「マルコフ連鎖モンテカルロ（MCMC）法とベイズ統計モデル」

書籍の著者　久保拓弥 先生

書籍「データ解析のための統計モデリング入門」８章「マルコフ連鎖モンテカルロ（MCMC）法とベイズ統計モデル」の Python写経活動記録 です。 

書籍は第８章で ベイズ統計モデルのパラメータ推定方法である MCMC法 を学びます。

この記事は MCMC アルゴリズムのひとつ「メトロポリス法」の実装を通じて、「MCMC サンプルは確率分布からの乱数だ」を実感できるように取り組みます。

「最尤推定」から「ベイズ・MCMC」へとステップ・バイ・ステップで学びを進めます。
ChatGPT に頼りっきりの記事ですが、どうぞお読みください。

では書籍を開いて統計モデリングの旅に出かけましょう🚀

はじめに

このブログシリーズは、書籍「データ解析のための統計モデリング入門 一般化線形モデル・階層ベイズモデル・MCMC」（岩波書店、「テキスト」と呼びます）の Python 写経を通じて得た「統計モデリングの楽しさ」をご紹介します。

テキストの紹介と引用表記はリンク先の記事に掲載しています。

準備

準備

■ 記事の範囲
この記事はテキスト８章の以下の節を取り扱います。

■ ライブラリのインポート
この記事で用いるライブラリをインポートします。

# インポート

# 数値計算
import numpy as np

# 統計計算
import scipy.stats as stats

# 可視化
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib

統計モデリング・サマリー

この記事で扱う統計モデリングの概要です。

◼️統計モデル
二項分布に従う例題データの生存確率パラメータ $${q}$$ を４つの方法で推定します。
上から順に段階的に進めていきます。

データの準備

◼️データの登録
テキスト p.171 の例題データを numpy配列に登録します。

# 例題の架空データの登録 p.171
data = np.array([4, 3, 4, 5, 5, 2, 3, 1, 4, 0, 1, 5, 5, 6, 5, 4, 4, 5, 3, 4])
print('data.shape:', data.shape)

【実行結果】
データの個数（標本サイズ）は 20 です。
植物の個体（個体数 20）に関する仮想の観測データです。
個体 $${i}$$ の調査種子数 $${N=8}$$ のうち生存種子数 $${y_i}$$ が記録されています。

テキストによると、例題データは「真の」生存確率 $${q=0.45}$$ の二項分布から生成された乱数です。

◼️データの可視化
例題データのヒストグラムと「真の生存確率の二項分布」による生存種子数の折れ線グラフを重ね描きします。
テキスト p.172 図 8.1 (B) に相当します。

# 例題の架空データのヒストグラム p.172 図8.1(B)

# 設定
N = 8                 # 調査種子数
q_true = 0.45         # 真の生存確率
n = len(data)         # サンプルサイズ
x_val = np.arange(9)  # x軸の値

# 例題データのヒストグラムの描画
sns.histplot(data, bins=x_val - 0.5, edgecolor='white', label='観測値')

# q=0.45の二項分布に従う真の個体数分布の描画
plt.plot(x_val, stats.binom.pmf(k=x_val, n=N, p=q_true) * n, '--o', ms=8,
         color='tab:red', mfc='white', label='真の生存確率による分布')

# 修飾
plt.title(f'例題データのヒストグラムと $q={q_true}$ の二項分布')
plt.xlabel('生存種子数 $y_i$', fontsize=12)
plt.ylabel('個体数', fontsize=12)
plt.legend();

【実行結果】
生存種子数が二項分布に従っていると感じられるチャートです。

① 最尤推定・解析解

概要

テキスト p.173 の 二項分布 $${\text{binomial}(N, q)}$$ のパラメータ $${q}$$ の最尤推定値 $${\hat{q}}$$ の計算式

の最尤推定量の数式を導出します。

まず生存確率 $${q}$$ と対数尤度 $${\log L(q)}$$ のチャートを描画します。
最尤推定量の数式の導出「5. （確認）二次導関数と最大性」で確認する「上に凸」（凹関数）の見える化です。
テキスト p.172 図 8.2 に相当します。

# 生存確率qと対数尤度logL(q)の可視化 p.172 図8.2

## 設定
qs1 = np.arange(0.26, 0.64, 0.00001)  # x軸の生存確率qの値
log_Ls1 = np.zeros(len(qs1))          # 対数尤度を格納する配列の初期化

## 対数尤度の算出
# 生存確率qごとに対数尤度を算出
log_Ls1 = np.array([np.sum(stats.binom.logpmf(k=data, n=N, p=q)) for q in qs1])
# 最大対数尤度のインデックスを取得
max_idx = np.argmax(log_Ls1)

## 描画
# 描画領域の設定
fig, ax = plt.subplots()
# 対数尤度の曲線の描画
ax.plot(qs1, log_Ls1)
# 最大対数尤度のｘ点の描画
ax.plot(qs1[max_idx], log_Ls1[max_idx], 'x', color='tab:red', ms=10)
# 最大対数尤度の垂直点線の描画
ax.axvline(qs1[max_idx], color='black', ls='--')
# 修飾
ax.set_xlabel('生存確率 $q$', fontsize=12)
ax.set_ylabel('対数尤度', fontsize=12)
ax.set(xlim=(0.25, 0.65), ylim=(-55, -35),
       xticks=[0.3, 0.4, 0.5, 0.6], yticks=[-50, -45, -40],
       title='生存確率 $q$ と対数尤度 log $L(q)$\n'
             f'対数尤度が最大となる生存確率 $q$={qs1[max_idx]:.5f}');

【実行結果】
対数尤度の曲線は上に凸であり、傾き＝０の点（❌️点）が対数尤度最大値になっていることが分かります。

最尤推定量の数式の導出

二項分布のパラメータの最尤推定量の数式展開を ChatGPT が丁寧に解説してくれました。

二項分布 $${\mathrm{Bin}(N,q)}$$ の成功確率 $${q}$$ をデータ $${y_1,\dots,y_n}$$ から最尤推定する流れを、式変形を1行ずつ丁寧に示します。
（前提：各 $${y_i\in{0,\dots,N}}$$、観測は独立、$${n\ge1}$$、$${N\ge1}$$、微分では $${0<q<1}$$ を仮定します。）

1. 尤度関数 $${L(q)}$$ の展開

各観測の確率質量関数は

独立より標本全体の尤度は積：

積の中身を要素ごとに分けます（積の分配法則）：

指数法則 $${\prod_{i} a^{b_i}=a^{\sum_i b_i}}$$ を使うと

ここで

なので、記号 $${Y:=\sum_{i=1}^n y_i}$$ を導入すると

🍀🍀🍀

2. 対数尤度 $${\ell(q)=\log L(q)}$$ の展開

対数の性質 $${\log(ab)=\log a+\log b}$$、$${\log(a^b)=b\log a}$$ を順に適用します：

（$${\sum_{i=1}^n \log\binom{N}{y_i}}$$ は $${q}$$ に依存しない定数項。）

🍀🍀🍀

3. $${q}$$ に関する微分 $${\ell'(q)}$$ の計算

微分公式と連鎖律を明示します：

これを項ごとに適用：

ゆえに

🍀🍀🍀

4. 一次条件 $${\ell'(q)=0}$$ の解法（最尤推定量）

右辺を移項して

両辺に $${q(1-q)\ (>0)}$$ を掛けます（不等号でないので符号は気にしなくてよい）：

左辺を展開し、右辺も分配則で展開：

両辺に $${Yq}$$ を加えて $${ -Yq}$$ を消します：

両辺を $${nN}$$ で割る（$${nN>0}$$）：

パラメータ推定量 $${\hat{q} = \frac{\sum_{i=1}^n y_i}{nN}}$$ を導出できました！

🍀🍀🍀

5. （確認）二次導関数と最大性

二次導関数も明示します：

従って $${\ell(q)}$$ は $${q}$$ に関して凹関数で、臨界点は最大です。
テキスト p.172 図 8.2 の対数尤度のチャートは「上に凸の形状」になっています。

🍀🍀🍀

6. 境界ケース（完全成功・完全失敗）

• $${Y=0}$$（全観測が失敗）のとき、対数尤度は $${q\downarrow 0}$$ で最大化され $${\hat q=0}$$。

• $${Y=nN}$$（全観測が成功）のとき、$${q\uparrow 1}$$ で最大化され $${\hat q=1}$$。

• $${0<Y<nN}$$ のとき、上の内部解 $${\hat q=Y/(nN)}$$ が成り立ちます。

以上が、尤度 ⇒ 対数尤度 ⇒ 微分 ⇒ 0 をおいて解くまでの完全展開です。

大ボリュームですね…
でも数式の展開を１つ１つ丁寧に行ってくれたので、理解しやすかったです。

二項分布 $${\text{binomial}(N, q)}$$ のパラメータ $${q}$$ の最尤推定値 $${\hat{q}}$$ は

となりました。
例題データは調査種子数 $${N=8}$$、標本サイズ $${n=20}$$ です。
最尤推定値 $${\hat{q}}$$ を計算しましょう。

# 二項分布のパラメータqの最尤推定値の算出 p.173 の数式より
# N = 8

q_hat = sum(data) / (N * len(data))
print('q の最尤推定値: ', q_hat)

【実行結果】
テキスト p.173 の計算結果と一致しました！

② 最尤推定・数値解

パラメータ推定の手順

テキスト 8.2 節のふらふら試行錯誤による最尤推定を python で実践します。

次の手順でパラメータ $${q}$$ を推定します。
なお、$${q}$$ を離散化している点にご留意下さい。

手順 (2) で $${q}$$ を増やす・減らす量（ステップ間隔）は 0.01 です。
0.01 ずつ $${q}$$ が動いて対数尤度が移り変わる様子をテキストの図 8.3、図 8.4 で確認します。

ｑを動かすイメージ

0.01 のステップ $${q}$$ で動かすイメージを掴みます。
0.01 刻みの $${q}$$ と対数尤度を可視化します。
テキスト p.174 図 8.8 に相当します。

二項分布の対数尤度は、データ点ごとに scipy.stats の binom.logpmf() で対数確率を求めて、足し合わせます。

# 生存確率qを離散化して各qにおける対数尤度をプロット p.174 図8.3

## 設定と準備
# x軸の生存確率qの値
qs2 = np.arange(0.01, 1.00, 0.01)
# 各qにおける対数尤度の算出
log_Ls2 = np.array([np.sum(stats.binom.logpmf(k=data, n=N, p=q)) for q in qs2])

## 描画
# 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 4))
# 生存確率と対数尤度のデータ点の描画
ax.plot(qs2, log_Ls2, 'o', ms=6, mfc='white')
# 修飾
ax.set_xlabel('生存確率 $q$', fontsize=12)
ax.set_ylabel('対数尤度', fontsize=12)
ax.set(xlim=(0.25, 0.65), ylim=(-54, -36),
       xticks=[0.3, 0.4, 0.5, 0.6], yticks=[-50, -45, -40]);

【実行結果】

パラメータ推定の手順では、ランダムに移動先候補を $${-0.01}$$ か $${0.01}$$ に決めて、移動先候補の対数尤度が現地点の対数尤度より大きい場合に移動をします。
ですので、１つのステップで「対数尤度が大きい⇒移動する」か「対数尤度が小さい⇒現地点に留まる」という動きをします。
ステップを繰り返すことで、対数尤度最大の $${q=0.46}$$ あたりに到達して、以後は留まり続けようとします。

続いてテキスト p.174 図 8.4 を描画します。
$${q}$$ が $${0.28}$$ ～ $${0.32}$$ のあたりをクローズアップして、「対数尤度の大きい $${q}$$」へ移動する様子を確認します。

# 図8.3のq=0.30付近を拡大した図 p.174 図8.4

## 設定
idx = np.where(np.isclose(qs2, 0.30))[0][0]  # q=0.30のインデックス
plot_qs = qs2[idx-2:idx+3]                   # 描画対象のq
plot_log_Ls = log_Ls2[idx-2:idx+3]           # 描画対象のlog L(q)
colors = ['lightblue', 'tab:blue']           # 矢印の色
adjust_x, adjust_y = 1.006, 0.995            # 矢印の始点の調整率

## 描画
# 描画領域の設定
fig, ax = plt.subplots()
# 矢印の描画
for i in range(4):
    ax.quiver(plot_qs[i] * adjust_x,
              plot_log_Ls[i] * adjust_y,
              plot_qs[i+1] - plot_qs[i],
              plot_log_Ls[i+1] - plot_log_Ls[i],
              angles='xy', scale_units='xy', scale=1.5, color=colors[i%2])
# データ点に添える対数尤度の数値の表示
for q, log_L in zip(plot_qs[1:-1], plot_log_Ls[1:-1]):
    ax.vlines(q, -50, log_L, color='black', ls='--', lw=0.5)
    ax.hlines(log_L, 0.275, q, color='black', ls='--', lw=0.5)
    ax.text(x=q-0.005, y=log_L+0.25, s=f'${log_L.round(2)}$')
# ５つのデータ点の描画
ax.plot(plot_qs, plot_log_Ls, 'o', ms=20, mfc='white')
# q=0.30のデータ点（青色で塗りつぶし）の描画
ax.plot(qs2[idx], log_Ls2[idx], 'o', ms=20, mec='tab:blue', mfc='lightblue')
# 修飾
ax.set_xlabel('生存確率 $q$', fontsize=12)
ax.set_ylabel('対数尤度', fontsize=12)
ax.set(xlim=(0.278, 0.322), ylim=(-49.2, -43.9),
       xticks=[0.28, 0.29, 0.30, 0.31, 0.32]);

【実行結果】
パラメータ推定の手順では「対数尤度が小さくなる移動はしない」ので、ステップを繰り返すことで、対数尤度の大きな $${q}$$ へ移動します。

最尤推定シミュレーションの実行

それではパラメータ推定の手順で $${q}$$ の最尤推定を行う「シミュレーション」を開始します！

パラメータ推定の手順で最尤推定値を探索する関数を定義します。
各ステップの動きの履歴を保存して後から参照できるようにします。

# 試行錯誤による対数尤度最大化にともなうqを算出する関数 p.175
def search_optim_prob(data, start_q, N=8, step=0.01, iter=100, seed=123):

    ## 設定と準備
    # 初期化
    result = np.zeros((iter + 1, 2))      # 結果を格納する配列
    rng = np.random.default_rng(seed)     # 乱数生成器
    # 初期設定
    current_q = start_q                   # スタート時の現在のq
    current_log_L = np.sum(               # スタート時の現在のqの対数尤度
        stats.binom.logpmf(k=data, n=N, p=current_q))
    result[0] = current_q, current_log_L  # スタート時点の結果を格納
    
    ## 試行錯誤をiter回繰り返し処理
    for i in range(iter):
        # となりの生存確率qをランダムに決定 ※q-0.01かq+0.01をランダムに選択
        next_q = rng.choice(a=[current_q - step, current_q + step])
        # となりの生存確率qの対数尤度を算出
        next_log_L = np.sum(stats.binom.logpmf(k=data, n=N, p=next_q))
        # となりの生存確率qの対数尤度のほうが大きい場合、現在のqをとなりのqに変更
        if current_log_L < next_log_L:
            current_q, current_log_L = next_q, next_log_L
        # このステップの結果を格納
        result[i + 1] = current_q, current_log_L
    
    # 戻り値：shape=(iter+1, 2), １列目：qの推移、２列目：qの対数尤度
    return result

【実行結果】なし

ではシミュレーションの開始です！
テキストにならって、$${q}$$ の開始点を $${0.30, 0.60}$$ の２パターンで実践します。
各ステップの動きについては、テキスト p.175 図 8.5 に相当する可視化で確認します。

# 試行錯誤による対数尤度最大化にともなうqの変化の可視化 p.175 図8.5

## q=0.30, 0.60の生存確率qデータを作成
result_030 = search_optim_prob(data, 0.30)
result_060 = search_optim_prob(data, 0.60)

## 描画
# 描画領域の設定
ax = plt.subplot()
# q=0.30の折れ線グラフ（青色）の描画
ax.plot(result_030[:, 0], label='$q=0.30$から開始した試行錯誤')
# q=0.60の折れ線グラフ（赤い点線）の描画
ax.plot(result_060[:, 0], color='tab:red', ls='--',
        label='$q=0.60$から開始した試行錯誤')
# 修飾
ax.set_xlabel('試行錯誤のステップ数', fontsize=12)
ax.set_ylabel('生存確率 $q$', fontsize=12)
ax.set(ylim=(0.26, 0.64), yticks=[0.3, 0.4, 0.5, 0.6])
ax.legend();

【実行結果】

【考察】
ステップの前半で最大対数尤度の $${q=0.46}$$ に到達して、以後は移動しない結果となりました。
$${q=0.60}$$ から始めた赤点線は、20ステップあたりで最大対数尤度に達したようです。
$${q=0.30}$$ から始めた青実線は、35ステップあたりで最大対数尤度に達したようです。

$${q=0.60}$$ から始めたケースの ステップ 20 ～ 29 の $${q}$$ の位置と対数尤度を表示しましょう。

# q=0.60の20～29ステップの生存確率と対数尤度
result_060[20:30]

【実行結果】
0.01 ずつ $${q}$$ の値が小さくなり（ときどき変化無しで）、最尤推定値 $${\hat{q}=0.46}$$ に到達しています。

推定したいパラメータ $${q}$$ をランダムに両隣へ「一定の幅」で移動を試みて、対数尤度が大きくなる方向に $${q}$$ を「一定の幅」で動かし、最大対数尤度の地点で足踏みさせる、こんなシミュレーションでした。

このような動きのイメージを押さえて、MCMC へ進みます！

③ 尤度・MCMC

パラメータ推定の手順

テキスト 8.3 節の MCMC アルゴリズムを python で実装します。
「メトロポリス法」と呼ばれるアルゴリズムで、最尤法ではない方法によるパラメータ推定を行います。

MCMC はマルコフ連鎖モンテカルロ法（Markov Chain Monte Carlo method）の略称です。
テキストの MCMC 関連用語の説明を引用いたします。

次の手順＝メトロポリス法でパラメータ $${q}$$ を推定します。
こちらも $${q}$$ を離散化します。
１回のステップ間隔は $${0.01}$$ です（先のシミュレーションと同じ）。
推定の結果、MCMC サンプルを得ます。

先の最尤推定・数値解とよく似た手順です。
太字が追加・変更点です。
大きな変更は (4) の手順が追加されたことです。

$${q^{新}}$$ の尤度が小さい場合も、確率 $${r}$$ で $${q^{新}}$$ に変更します。
つまり、尤度が悪化する方向の移動が追加されたのです。
悪化する方向へ移動する確率は $${r}$$ です。

📊 確率 $${r}$$ の公式
$${r}$$ は $${q}$$ と $${q^{新}}$$ の尤度比です。

$${q}$$ と $${q^{新}}$$ の対数尤度の差が小さいほど $${r}$$ は大きくなり、移動しやすくなります。

メトロポリス法の実行

それではメトロポリス法で $${q}$$ のパラメータ推定の「シミュレーション」を開始します！

メトロポリス法で MCMC サンプリングを行う関数を定義します。
各ステップの $${q}$$ の履歴を保存して後から参照できるようにします。
この履歴が「MCMC サンプル」です。

# メトロポリス法によるMCMCサンプリングの例 p.179 図8.8

# 二項分布の対数尤度算出関数の定義
def log_likelihood_binom(Y, T, q):
	"""二項分布の対数尤度（定数は省略してOK。ここでは分かりやすく q 依存部分のみ）
	L(q) ∝ q^Y * (1-q)^(T - Y) なので
	log L(q) = Y * log q + (T - Y) * log(1 - q)
	数値安定のため log(1-q) は np.log1p(-q) を使用
	"""
	if q <= 0.0 or q >= 1.0:
		return -np.inf
	# 戻り値：対数尤度
	return Y * np.log(q) + (T - Y) * np.log1p(-q)

# MCMCサンプリング関数の定義 ※テキストのRスクリプトとは異なる実装です
def exec_metropolis(data, start_q=None, N=8, interval=0.01, iter=100, seed=42):

	## 設定と準備
	Y = sum(data)
	T = len(data) * N
	# 初期化
	result = np.zeros((iter + 1, 2))                      # 結果を格納する配列
	rng = np.random.default_rng(seed)                     # 乱数生成器
	# 1.パラメータqの初期を決める
	if start_q is None:  # 引数未指定：開始時のqをランダムに設定
		current_q = round(rng.random(), 2) 
	else:                # 引数指定　：指定された開始時のqを設定
		current_q = start_q
	# その他の初期値設定
	current_log_L = log_likelihood_binom(Y, T, current_q) # 開始時のqの対数尤度
	result[0] = current_q, current_log_L                  # 開始時の結果を格納

	## 試行錯誤をiter回繰り返し処理
	for i in range(iter):
		# 2.qを増やすか減らすかをランダムに決定 ※q-interval,q+intervalをランダム選択
		new_q = rng.choice(a=[current_q - interval, current_q + interval])
		#   q新の対数尤度の算出
		new_log_L = log_likelihood_binom(Y, T, new_q)
		# 3.q新の対数尤度が大きい場合、現在のqをq新に更新 ※等しい場合もq新で置き換え
		if current_log_L <= new_log_L:
			current_q, current_log_L = new_q, new_log_L
		# 4.q新の対数尤度が小さい場合、確率rで現在のqをq新に更新するかどうか決定
		elif current_log_L > new_log_L:
			# 確率（尤度比）rの算出 ※テキストp.177「移動する確率」の計算式
			r = np.exp(new_log_L - current_log_L)
			# 確率rでq新を選択して、現在のqをq新に更新
			if rng.random() <= r:
				current_q, current_log_L = new_q, new_log_L
		# このステップの結果を格納
		result[i + 1] = current_q, current_log_L

	# 戻り値：shape=(iter+1, 2), １列目：qの推移、２列目：qの対数尤度
	return result

【実行結果】なし

【コード補足説明：二項分布の対数尤度算出関数】
ChatGPTがこの関数の原型を作りました。
二項分布の対数尤度関数を次の式で実装しています。
ここで、$${Y=\sum_i y_i, T = nN}$$ です。

では、MCMC サンプリングを行います。
サンプルサイズを 100, 1000, 100000 の３ケースで行います。

%%time
# 描画用データの作成 p.179 図8.8, p.181 図8.9, p.183 図8.10, p.184 図8.11

# 関数を用いて、q=0.30からスタートするMCMCサンプルを100, 1000, 100000個作成
result_100 = exec_metropolis(data, 0.30, iter=100)
result_1000 = exec_metropolis(data, 0.30, iter=1000)
result_100000 = exec_metropolis(data, 0.30, iter=100000)

# マルコフ連鎖の定常分布データの作成 ※テキストのRスクリプトのコードを引用
# 設定
Y = sum(data)                                       # 全生存種子数
T = len(data) * 8                                   # 全調査種子数
interval = 0.01                                     # 生存確率qの間隔
# 生存確率q（等間隔）の作成 ※部分区間を切り出し
vq = np.arange(0.25, 0.65 + interval, interval)
# 上記vqに対応する対数尤度の算出
# vlogL = [log_likelihood_binom(Y, T, q) for q in vq]
vlogL = [sum(stats.binom.logpmf(k=data, n=8, p=q)) for q in vq]  # 図8.10対策
# 尤度の確率 ※p.183のp(q|Y)の数式
density_vL = np.exp(vlogL) / sum(np.exp(vlogL))

【実行結果】
割と短い時間で処理が完了しました。

パラメータ $${q}$$ の MCMC サンプルを可視化します。
テキスト p.179 の図 8.8 に相当します。

# グラフの描画 p.170 図8.8

## 描画用の設定 ※３つのデータの設定のリスト化
# MCMC stepのx軸（step数）のデータ
x_vals = [range(len(result_100)), range(len(result_1000)),
          np.arange(0, len(result_100000), 100)]
# MCMC stepのy軸（生存確率）のデータ
y_vals = [result_100[:, 0], result_1000[:, 0], result_100000[0::100, 0]]
# サンプリングされたqのヒストグラムのデータ
hist_vals = [result_100[:, 0], result_1000[:, 0], result_100000[:, 0]]
# ヒストグラムのbin
binses = [20, 26, 32]
# グラフタイトルに使用する値
title1s = ['A', 'B', 'C']
title2s = [len(result_100)-1, len(result_1000)-1, len(result_100000)-1]

## 描画領域の指定
# 全体の設定
fig = plt.figure(figsize=(8, 7))
# Grid領域の設定
grid = plt.GridSpec(3, 6)

# ３つのデータの描画を繰り返し処理
for i, (x_val, y_val, hist_val, bins, title1, title2) in enumerate(
    zip(x_vals, y_vals, hist_vals, binses, title1s, title2s)):
    
    ## 描画領域の設定
    ax1 = fig.add_subplot(grid[i, 0:5])            # MCMC stepの描画領域
    ax2 = fig.add_subplot(grid[i, 5], sharey=ax1)  # ヒストグラムの描画領域

    ## プロット
    # MCMC stepの描画
    ax1.plot(x_val, y_val, lw=1)
    # ヒストグラムの描画
    ax2.hist(hist_val, bins=bins, density=True, orientation='horizontal',
             edgecolor='white')
    # マルコフ連鎖の定常分布の描画 ※×100でスケール合わせ
    ax2.plot(density_vL*100, vq, color='tab:red', lw=2.5)
    # 修飾1: 左のグラフのx・y軸の範囲、タイトル
    ax1.set(xlim=(0, title2), ylim=(0.25, 0.65),
            title=f'({title1}) {title2} MCMC stepまで')
    # 修飾2: グラフの行ごとに変える修飾
    match i:
        case 0:  # 1行目のヒストグラム等の文章の表示
            ax2.text(x=0.8, y=0.58, s='マルコフ連鎖の\n定常分布', color='tab:red')
            ax2.text(x=6, y=0.28, s='サンプルされた\n$q$ のヒストグラム')
        case 1:  # 2行目のグラフのylabelの表示
            ax1.set_ylabel('生存確率 $q$', fontsize=12)
        case 2:  # 3行目のグラフのxlabelの表示
            ax1.set_xlabel('MCMC step 数', fontsize=12)
    # 修飾3: ヒストグラムの枠線を消去
    ax2.axis('off')
# 修飾
plt.subplots_adjust(hspace=0.5, wspace=0.1)
plt.show()

【実行結果】
なんだか不思議な図が現れました…

【考察】
左側の枠内のギザギザ線は $${q}$$ のMCMC サンプルの軌跡になっています。
尤度が大きくなる方向だけでなく、尤度が小さくなる方向にも動いています。
また、最適なパラメータ値と考えられる位置に留まること無く、常にフラフラ動いています。

右側の枠外のヒストグラムは $${q}$$ のサンプルの度数分布です。
赤い線は「定常分布」です。
定常分布は、この例題の統計モデルとメトロポリス法によって決まるマルコフ連鎖の確率分布です。
MCMC ステップ数が大きくなるにつれて、$${q}$$ の度数分布が定常分布に近づく様子が分かります。

定常分布とMCMCサンプル

◼️定常分布
大切なお知らせがあります📣

テキストによると、観測データ $${\bm Y}$$ が得られたときの $${q}$$ の定常分布 $${p(q \mid \bm Y)}$$ は、

この例題に関しては、定常分布 $${p(q \mid \bm Y)}$$ は尤度 $${L(q)}$$ に比例する離散化した $${q}$$ の確率分布であり、次の式で示されます。

📈 定常分布の公式（この例題の場合）
$${q}$$ を離散化しているので $${\sum}$$ を使っています。

そしてメトロポリス法によって得られたMCMC サンプルは、定常分布からのランダムサンプルなのです。

🍀🍀🍀

◼️定常分布と MCMC サンプルのつながり
テキスト p.181 の図 8.9 で定常分布と MCMC サンプルの関連を確認します。
「７本」のマルコフ連鎖（チェーン）でサンプリングします。
１回のステップ間隔は $${0.01}$$ です（先のシミュレーションと同じ）。

# いろいろなqの値から開始した,メトロポリス法によるサンプリングの反復 p.181 図8.9

## 描画領域の設定
# 全体の設定
fig = plt.figure(figsize=(10, 4))
# Grid領域の設定
grid = plt.GridSpec(1, 10, wspace=0.1)         # GridSpecの設定
ax1 = fig.add_subplot(grid[0, 0:9])            # MCMC stepの描画領域
ax2 = fig.add_subplot(grid[0, 9], sharey=ax1)  # 定常分布の描画領域

## 描画処理
# 開始するqごとにMCMCサンプリング実行と折れ線グラフ描画を繰り返し処理
for start_q in [0.25, 0.32, 0.39, 0.46, 0.53, 0.6, 0.65]:
    # 開始するqを用いてMCMCサンプリングを実行
    samples = exec_metropolis(data, start_q, iter=500, seed=int(start_q*100))
    # 取得したMCMCサンプルデータの折れ線グラフを描画
    ax1.plot(samples[:, 0], lw=1, alpha=0.7)
# step100あたりの目印の描画
ax1.axvline(100, color='gray', ls='--')
# 定常分布の描画
ax2.plot(density_vL, vq, color='tab:red', lw=2)
# 修飾（左グラフ）
ax1.set(yticks=(0.3, 0.4, 0.5, 0.6), xlim=(0, 500))
ax1.set_xlabel('MCMC step 数', fontsize=12)
ax1.set_ylabel('生存確率 $q$', fontsize=12)
# 修飾（右グラフ）
ax2.text(x=0.01, y=0.63, s='定常分布', fontsize=12)
ax2.axis('off');

【実行結果】
ベイズ統計モデルのトレースプロットのようになりました。

step 100 に目印の点線を付けています。
各 MCMC サンプルは step 100 あたりまでは初期状態に影響されて「外れた値」になっていますが、step 100 あたりから定常分布に「収まっている」印象です。
定常分布目線では、step 100 あたりから「定常分布からのサンプル」のように見えます。

テキストによると「初期状態を捨てる、複数の MCMC サンプルを比較するなどで定常分布の推定の改善につながる」ようです。

実際にベイズ統計モデルでは、初期状態をバーンイン・ウォームアップ・ チューニングなどと呼んでいて、サンプルから除外します。
また、マルコフ連鎖の数を複数（４本など）にしてサンプルの偏りをなくす努力をします。

🍀🍀🍀

◼️尤度関数と定常分布のつながり
この例題の定常分布の公式のとおり、定常分布 $${p(q \mid \bm Y)}$$ が尤度 $${L(q)}$$ と比例関係にある様子を可視化します。
テキスト p.183 図 8.10 に相当します。

# 尤度関数と定常分布の関係 p.183 図8.10

# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 3))

# (A)qの尤度関数の描画
ax1.plot(vq, np.exp(vlogL), '-o', ms=7, mfc='white')
ax1.set_xlabel('生存確率 $q$', fontsize=12)
ax1.set(yticks=[0, 1.98*1e-17, 3.97*1e-17], title='(A) $q$ の尤度関数 $L(q)$')

# (B)qの定常分布の描画
ax2.plot(vq, density_vL, '-o', ms=7, mfc='white')
ax2.set_xlabel('生存確率 $q$', fontsize=12)
ax2.set(yticks=[0, 0.05, 0.10], title=r'(B) $q$ の定常分布 $p(q|\boldsymbol{Y})$');

【実行結果】
分布の形状が瓜二つです！

🍀🍀🍀

◼️さらに定常分布からのサンプリングを体感
MCMC サンプルは右の $${q}$$ の定常分布から生成された乱数ということのイメージをいっそう高めるための可視化に続きます。
定常分布とMCMCサンプルの分布（ヒストグラム）を重ね描きします。
100000 個のMCMCサンプルを使います。
テキスト p.184 図 8.11 に相当します。

# メトロポリス法によるデータへのあてはめによって得られたパラメータqの確率分布
# p.184 図8.11

# 描画領域の設定
fig, ax = plt.subplots()
# MCMCサンプリングの標本のヒストグラムの描画
sns.histplot(result_100000[:, 0], bins=32, stat='probability', ec='white',
             label='MCMCサンプル', ax=ax)
# 尤度L(q)に比例する定常分布p(q|Y)の折れ線グラフの描画
ax.plot(vq, density_vL, '-o', ms=7, color='tab:red', mfc='white',
        label='定常分布')
# 修飾
ax.set_xlabel('生存確率 $q$', fontsize=12)
ax.set(ylim=(-0.003, 0.11), xticks=(0.3, 0.4, 0.5, 0.6), yticks=(0, 0.05, 0.10))
ax.legend();

【実行結果】
MCMCサンプルの分布は定常分布とほぼほぼ致しています！
MCMCサンプルのヒストグラムはベイズ統計モデルのパラメータ事後分布プロットのようになりました。

MCMC サンプルは右の $${q}$$ の定常分布から生成された乱数といえそうですね！

🍀🍀🍀

◼️MCMC サンプルを活用するパラメータ推定
MCMC サンプリングによるパラメータ推定では、最尤推定のように点推定結果を直接得るのではなく、定常分布から得たサンプルからパラメータの点推定値や区間推定を行います。
たとえば、サンプルの平均（または中央値）＝パラメータ点推定値、サンプルの95%区間＝パラメータの区間推定値、といった感じです。

Python で体感してみます。
100000 個のMCMCサンプルのうち、最初の 1000 個を「初期状態なので捨てる」ことにして、残りの 99000 個で統計量を計算します。

# 100000個のMCMCサンプルからパラメータqの統計量を表示

# 設定
burnin = 1000
# 中央値、95%区間の値の取得
quantiles = np.quantile(result_100000[burnin:, 0], [0.5, 0.025, 0.975])

# データフレーム化
pd.DataFrame({
    'mean': [result_100000[burnin:, 0].mean()],
	'std': [result_100000[burnin:, 0].std()],
	'median': [quantiles[0]],
	'95%CI lower': [quantiles[1]],
	'95%CI upper': [quantiles[2]],
	}, index=pd.Series(['q'], name='パラメータ')
).round(4)

【実行結果】
PyMC のサマリーに似た表になりました。

【読み取り】
パラメータ $${q}$$ の点推定値は平均 $${0.4583}$$ で、95%区間は $${[0.38, 0.53]}$$ です。
なお、ステップ間隔 $${0.01}$$ でサンプルを生成したので、中央値や 2.5%点・97.5%点が小数第２位ぴったりの離散値になっています。

定常分布という「ほしい確率分布」から「MCMCサンプル」を得て、パラメータを推定する、といった流れを噛み締めながら、いよいよベイズ統計モデルの第一歩を踏み出します👣

④ ベイズ・MCMC

概要

ここからはベイズ統計の話になります。
ベイズ統計では「パラメータ $${q}$$ は確率変数」なんです！
堂々とパラメータ $${q}$$ の確率を語れます！

引き続き $${q}$$ を離散化します。
$${q}$$ を離散化しているので $${\sum}$$ を使えています。

ベイズ統計と言えば「ベイズの定理」です。
観測データ $${\bm Y}$$ が得られたときのパラメータ $${q}$$ の条件付き確率 $${p(q \mid \bm Y)}$$ を、ベイズの定理は次の式で表します。

📈 ベイズの定理（例題の場合）

各項には名前が付いています。
数式と名前をマッピングします。

$${p(\bm Y \mid q)}$$ は尤度なので $${L(q)}$$ と書き換えましょう。

この数式は「③ 尤度・MCMC」の「定常分布の公式」によく似ています。

テキストによるとこの「尤度・MCMC」のケースでは、「事前分布 $${p(q)}$$ が $${q}$$ の値によらずに定数になっていると、辻褄があっているように見えます」とのこと。

ベイズ統計モデルでは、MCMC サンプリングによって「尤度 × 事前分布」に比例するパラメータの「事後分布」を推定します。

🍀🍀🍀

さて、ここからは ChatGPT が教えてくれたベイズ統計モデリングの基礎知識に進みたいと思います。
事後分布の比例の数式を次のように書き換えて、スタートします！

（注意）
こちらは趣味的な深堀りとコードです。
ご興味ない方はスルーしてくださって大丈夫です。

ChatGPTと一緒に学ぶベイズ統計モデルの基礎知識

ChatGPTにベイズ統計モデルの基礎知識と Python コードを書いてもらいました。

最初にベイズ統計の基礎知識で、ベイズ統計モデルと メトロポリス法による MCMC サンプルの関係性などを確認します。

難易度高めですので、読み飛ばしていただいても大丈夫です。

1️⃣ 目的（まず直感）

• 点推定（いちばん“ありそう”な $${q}$$）ではなく、$${q}$$ がどのあたりにどれくらいの確率であり得るか＝分布（事後分布）を知りたいことが動機です。

• ベイズでは次の確率分布を扱います。

• この分布＝事後分布が複雑な場合は、メトロポリス法で“その分布から”乱数のサンプルを集めます。

🍀🍀🍀

2️⃣ モデルの検討

2-1. 尤度 $${L(q)}$$

• 各観測データが従う確率分布は $${y_i\sim\mathrm{Bin}(N,q)}$$（互いに独立）です。

• 積をとると $${L(q)\propto q^{Y}(1-q)^{T-Y}}$$ です。
  例題データの場合、$${Y=\sum_i y_i=73,\ T=nN=160}$$ です。
  $${N=8}$$ は 調査種子数（試行回数）、$${n=20}$$ は標本サイズです。

🍀

2-2. 事前分布 $${\text{prior}(q)}$$ の例
二項分布の成功確率パラメータ（ここでは $${q}$$）が従う事前分布にはベータ分布を用いることが多いです。
ベータ分布の確率変数が０から１の値をとるので、確率 $${q}$$ を表すのに適しています。
ベータ分布のパラメータは $${\alpha, \beta}$$ です。
$${\alpha, \beta}$$ に事前知識を使わない「無情報事前分布」には、次のような値を使うことが多いようです。

• 一様事前分布：$${q\sim\mathrm{Beta}(\alpha=1,\beta=1)}$$

• Jeffreys事前分布：$${q\sim\mathrm{Beta}\left(\alpha=\tfrac{1}{2},\beta=\tfrac{1}{2}\right)}$$

🍀

2-3. 事後分布 $${p(q\mid \text{data})}$$

• ベータ分布は二項分布の共役事前分布なので、事前分布がベータ分布、観測値が二項分布に従う場合、事後分布がベータ分布になります。

• ここでは 円滑な学びに結ぶように MCMC（メトロポリス法）でこの理論的な事後分布と一致するサンプルを作ります。
  理論的な事後分布と MCMC サンプルから推定する事後分布を比較できます。

🍀🍀🍀

3️⃣ メトロポリス法（ランダムウォーク・メトロポリス法）

3-1. 目標密度（Target Density）
推定したいパラメータ $${q}$$ の確率分布＝事後分布 $${\pi(q)}$$ を定義します。

🍀

3-2. 提案分布の作り方（対称な提案分布）
尤度・MCMCのときには、$${q_{新}}$$ の候補探しの際に、１回のステップ間隔を $${0.01}$$ としました。
このベイズ・MCMC では 正規分布 $${\text{Normal} (0, \sigma=0.03)}$$ の乱数 $${\varepsilon}$$ を１回のステップ間隔にします。
$${q_{新}}$$ への間隔が等間隔ではなくなります。

範囲外（$${q_{\text{新}}\notin(0,1)}$$）の場合は即棄却（据え置き）します。

🍀

3-3. 受理確率 $${r}$$ の算出と受理する／しないの判定

• 比を作る：
  尤度 $${\times}$$ 事前分布の比 $${r}$$ を作ります。

• ベータ事前分布を使うと目標密度 $${\pi(q)}$$ は次のように表されます。

• よって $${r}$$ は：

• 数値安定のため $${r}$$ を対数 にします。

• 一様乱数 $${u\sim \text{Uniform}(0,1)}$$ を生成して

（対称提案なので提案分布の比 $${J}$$ は打ち消し合います。）

🍀

3-4. （追補）ベイズの公式と「エビデンス（分母）」を無視できる理由

ここまでの計算では「ベイズの公式の分母：エビデンス」がなかったかのように進めてきました。
この追補では、エビデンスを無視できる理由を確認します！

まずベイズの公式をそのまま書きます：

ここでエビデンス（周辺尤度）$${p(\text{data})}$$

は $${q}$$ に依存しない定数（正規化定数）です。

メトロポリス法で使う 目標密度 $${\pi(q)}$$ を

と書き、エビデンスを $${Z :=p(\text{data})}$$ とおくと、受理比 $${r}$$（対称提案での単純形）は

同じ定数 $${1/Z}$$（エビデンスの逆数）が分子と分母に出て打ち消し合うので、エビデンス $${Z}$$ を計算せずに済みます。
（提案が非対称のときは、さらに提案密度 $${J}$$ の比 $${J(q\mid q_{\text{新}})/J(q_{\text{新}}\mid q)}$$ を掛けますが、やはり $${Z}$$ は消えます。）

数値安定のため 対数で書くと

で、ここにも $${\log Z}$$ は一切現れません（差分で消えているから）。

（おまけ：同じ理屈で消える“定数”の例）

• 二項の対数尤度に出てくる $${\sum_i \log \binom{N}{y_i}}$$ も $${q}$$ に依存しない定数なので、$${\log r}$$ の差分で打ち消し合い、無視してよい（計算を軽く・安定にできます）。

🍀🍀🍀

4️⃣ 数値で一歩だけやってみる（一様事前分布：$${\alpha=\beta=1}$$）

4-1. うまく行きそうな小さな移動

• 現在 $${q=0.45}$$、提案 $${q_{\text{新}}=0.47}$$ とする。

• 計算：

• $${u}$$ が 0.953 以下なら受理（約95%の確率で進む）。

🍀

4-2. 大きく上へ飛びすぎた移動

• 同じ現在 $${q=0.45}$$ から $${q_{\text{新}}=0.55}$$ に飛ぶと：

• 受理は約 6%（ほぼ棄却）。
  ⇒ step（$${\sigma}$$）が大きすぎると受理率が下がる直感がつかめます。

🍀🍀🍀

5️⃣ 全体の流れのまとめ（アルゴリズムを文章で）

1. 初期値 $${q^{(0)}\in(0,1)}$$ を決める（例：0.5）。

2. $${t=0,1,2,\dots}$$ と繰り返す：
   2-1. $${q_{\text{新}}=q^{(t)}+\mathcal N(0,\sigma^2)}$$ を提案。
   2-2. $${q_{\text{新}}\notin(0,1)}$$ なら棄却し $${q^{(t+1)}=q^{(t)}}$$。
   2-3. そうでなければ $${\log r}$$ を上の式で計算。
   2-4. $${u\sim \text{Uniform}(0,1)}$$ を生成して、
   　　$${\log u\le\log r}$$ なら受理（$${q^{(t+1)}=q_{\text{新}}}$$）、
   　　否なら棄却（$${q^{(t+1)}=q^{(t)}}$$）。

3. はじめの数千ステップをバーンインとして捨て、その後を MCMC サンプル（または事後サンプル）として使う。

4. サンプルの平均・分位点 → 推定値・信用区間にする。

🍀🍀🍀

6️⃣ ここで覚えておくコツ

• 対数にすると「積→和」「小数の掛け算→加算」で安定。

• 受理率（accept rate）は 0.2〜0.5 くらいが目安。
  大きすぎ→歩幅小さすぎ、小さすぎ→歩幅大きすぎ。

• 今回の問題はそもそも解析解（ベータ事後分布）があるので、MCMCは学習用の確認として最適。

🍀🍀🍀

7️⃣ 超ミニ要約（3行）

1. 事後分布 $${p(q\mid\text{data}) \propto q^{\alpha-1}(1-q)^{\beta-1}}$$ を目標にする。

2. $${q_{\text{新}}=q+\mathcal N(0,\sigma^2)}$$ を提案し、$${\log r=(\alpha-1)(\log q_{\text{新}}-\log q)+(\beta-1)(\log(1-q_{\text{新}})-\log(1-q))}$$。

3. $${\log u\le\log r}$$ なら受理。
   これを繰り返して MCMC サンプルを集め、平均や区間を出す。

長文をお読みいただき、お疲れ様でした。
一休みしましょう。

ChatGPTと一緒に動かすベイズ MCMC

いよいよクライマックスです。

二項分布・ベータ分布のベイズ統計モデルで MCMC サンプリングを行い、パラメータ $${q}$$ を推定します！

◼️ベイズ統計のための関数定義
まずは例題用のベイズ統計モデル＝関数を定義します。
前半は二項分布の対数尤度、ベータ分布の事前分布の対数確率密度、事後分布の対数確率密度を計算する関数群です。
後半はメトロポリス法でMCMCサンプリングを行う関数です。
ロジックはベイズ統計モデルの基礎知識の内容を踏襲しています。

# 関数定義

# ============================================================
#  ベイズ推定（qの事後分布）をメトロポリス法でサンプリングする最小実装
#  - 観測: y_i ~ Binomial(N, q)
#  - 事前: q ~ Beta(a, b)   ← ここを変えれば事前を変更できます
#  - 目的: 事後 p(q | data) ∝ L(q) * Beta(q | a, b) からMCMCサンプルを得る
#  - 出力: 事後平均, 信用区間, 受理率など
# ============================================================

# ------------------------------------------------------------
#  対数尤度・対数事前・対数事後の関数
#  （数値安定のため log(1-q) は np.log1p(-q) を使用）
# ------------------------------------------------------------
def log_likelihood_binom_(Y, T, q):
	"""二項分布の対数尤度（定数は省略してOKだが、ここでは分かりやすく q 依存部分のみ）
	L(q) ∝ q^Y * (1-q)^(T - Y) なので
	log L(q) = Y * log q + (T - Y) * log(1 - q)
	"""
	if q <= 0.0 or q >= 1.0:
		return -np.inf
	return Y * np.log(q) + (T - Y) * np.log1p(-q)

def log_prior_beta_(q, a, b):
	"""Beta(a, b) の対数密度の（q依存）部分
	Beta(a, b) ∝ q^(a-1) * (1-q)^(b-1) なので
	log prior = (a-1) * log q + (b-1) * log(1 - q)
	"""
	if q <= 0.0 or q >= 1.0:
		return -np.inf
	return (a - 1.0) * np.log(q) + (b - 1.0) * np.log1p(-q)

def log_posterior(q, a, b, Y, T):
	"""対数事後 log p(q | data) ∝ log L(q) + log prior(q)"""
	return log_likelihood_binom_(Y, T, q) + log_prior_beta_(q, a, b)

# ------------------------------------------------------------
#  メトロポリス法（ランダムウォーク・メトロポリス法）
#  - 提案: q_new = q + Normal(0, step^2)
#  - 範囲外(0,1)は即棄却（初学者向けにシンプルに）
# ------------------------------------------------------------
def metropolis_q(
	Y, T,              # 総生存種子数Y, 総調査種子数T
	a, b,              # 事前分布・ベータ分布のハイパーパラメータa, b
	n_samples=20000,   # 取得したいサンプル数（バーンイン後）
	burnin=2000,       # バーンイン（慣らし）ステップ数
	step=0.03,         # 提案分布の標準偏差（大きい→探索広いが棄却↑）
	q0=0.5,            # 初期値（0と1は避ける）
	seed=0,            # 乱数シード（再現性のため）
):
	rng = np.random.default_rng(seed)
	q = float(q0)
	lp = log_posterior(q, a, b, Y, T)  # 現在点の対数事後
	samples = []                       # サンプルの蓄積先（バーンイン後のみ入れる）
	n_accept = 0                       # 受理回数（受理率の確認用）

	total_steps = n_samples + burnin
	for t in range(total_steps):
		# --- 提案: 正規ノイズを足す（ランダムウォーク） ---
		q_new = q + rng.normal(0.0, step)

		# --- (0, 1) の外に出たら棄却（据え置き） ---
		if 0.0 < q_new < 1.0:
			# 新しい点の対数事後を計算
			lp_new = log_posterior(q_new, a, b, Y, T)

			# 受理確率 r = exp(lp_new - lp)
			# log u <= log r と比較する（数値安定のため対数空間で判定）
			logu = np.log(rng.uniform())
			logr = lp_new - lp

			# 受理するか判定
			if logu <= logr:
				q = q_new
				lp = lp_new
				n_accept += 1
		# 棄却なら q は据え置き

		# バーンインを過ぎたら保存
		if t >= burnin:
			samples.append(q)

	# MCMCサンプルをnumpy配列化
	samples = np.array(samples, dtype=float)
	# 受理率の算出（目安: 20〜50%くらいが扱いやすい）
	accept_rate = n_accept / float(total_steps)

	# 戻り値：MCMCサンプル、受理率
	return samples, accept_rate

【実行結果】なし

ここで、事前分布に用いるベータ分布のパラメータを確認します。
① $${\text{Beta}(\alpha=1, \beta=1)}$$ と ② $${\text{Beta}(\alpha=1/2, \beta=1/2)}$$ の確率密度を比べます。

# 事前分布プロット

# 描画
# qの事後分布のKDEプロット
x_vals2 = np.linspace(0, 1, 101)
plt.plot(x_vals2, stats.beta.pdf(x_vals2, a=1, b=1), label='Beta(a=1, b=1)')
plt.plot(x_vals2, stats.beta.pdf(x_vals2, a=0.5, b=0.5), color='tomato', ls='--',
         label='Beta(a=0.5, b=0.5)\nジェフリーの事前分布')
# 修飾
plt.xlabel('$q$', fontsize=12)
plt.ylabel('確率密度', fontsize=12)
plt.legend();

【実行結果】
① $${\text{Beta}(\alpha=1, \beta=1)}$$ は 1 で一定です。
② $${\text{Beta}(\alpha=1/2, \beta=1/2)}$$ はＵ字型です。
今回の ベイズ統計モデルでは ① $${\text{Beta}(\alpha=1, \beta=1)}$$ を使いました。

◼️MCMC の実行
MCMC サンプリングを実行して、パラメータ $${q}$$ を推定します。
あわせて、ベータ事後分布の理論値（解析解）によるパラメータ推定値も確認します。

# MCMCの実行

# ------------------------------------------------------------
#  設定と準備
# ------------------------------------------------------------
# ---- 観測データ ----
y = data.copy()             # 生存種子数データ
N = 8                       # 各観測の試行回数（調査種子数）
n = y.size                  # 観測の個数
Y = int(y.sum())            # 成功回数の合計 Y = Σ y_i
T = n * N                   # 総試行回数 T = n * N

# ---- 事前分布のハイパーパラメータ ----
# 既定: 一様事前 Beta(1, 1)
# Jeffreys事前にしたい場合は a = b = 0.5 に変更してください
a, b = 1.0, 1.0

# ---- 参考: 解析的な事後（共役）パラメータ（検算用） ----
# 目的変数が二項分布(n, q)に従い、パラメータqの共役事前分布がベータ分布(a, b)のとき
# ⇒事後分布はベータ分布(a+x, b+n-x)。事後分布の期待値は (a+x) / (a+x + b+n-x)
alpha_post = a + Y          # 事後のα
beta_post  = b + (T - Y)    # 事後のβ
post_mean_theory = alpha_post / (alpha_post + beta_post)  # 事後平均の理論値


# ------------------------------------------------------------
#  実行（既定設定でOK。受理率が高すぎ/低すぎなら step を調整）
# ------------------------------------------------------------
samples, acc = metropolis_q(
	Y, T, a, b,       # 二項分布・ベータ分布のパラメータ等
	n_samples=20000,  # サンプルサイズ（chain=1）
	burnin=2000,      # バーンイン数
	step=0.03,        # 受理率が高すぎる(>0.6)→大きく、低すぎる(<0.2)→小さく
	q0=0.5,           # 初期値
	seed=42           # 乱数シード
)

# ---- 事後要約統計（MCMC推定値）----
q_mean = float(samples.mean())                       # 事後平均（MCMC）
q_median = float(np.percentile(samples, 50))         # 事後中央値（MCMC）
q_ci_95 = np.percentile(samples, [2.5, 97.5])        # 95%信用区間（MCMC）

# ---- 参考: 解析的（理論）な事後平均（共役の検算）----
# Beta(alpha_post, beta_post) の平均 = alpha_post / (alpha_post + beta_post)
q_mean_theory = post_mean_theory

# ---- 表示 ----
print('--- 基本情報 ---')
print(f'n = {n}, N = {N}, Y = {Y}, T = nN = {T}')
print(f'事前分布: Beta(a={a}, b={b})')
print(f'事後分布 (共役): Beta(alpha=a+Y={alpha_post}, beta=b+T-Y={beta_post})')
print()

print('--- MCMC 結果（メトロポリス法）---')
print(f'受理率(accept rate):\t{acc:.3f}')
print(f'事後平均(MCMC):\t\t{q_mean:.6f}')
print(f'事後中央値 (MCMC):\t{q_median:.6f}')
print(f'95%信用区間 (MCMC):\t[{q_ci_95[0]:.6f}, {q_ci_95[1]:.6f}]')
print()

print('--- 検算（解析解との比較）---')
print(f'理論の事後平均 (Beta):\t{q_mean_theory:.6f}')
print('※ MCMCの平均が理論値に近ければ、サンプルが素直に取れている目安になります。')

# ============================================================
#  使い方メモ
#  - Jeffreys事前にしたい→ a=b=0.5 に変える
#  - 受理率が低すぎる(<0.2)→ step を小さく
#  - 受理率が高すぎる(>0.6)→ step を大きく
#  - サンプル数を増やしてCIを安定させたい→ n_samples を大きく
# ============================================================

【実行結果】
MCMC サンプルの平均（事後平均）によるとパラメータの推定値 $${\hat{q}=0.457}$$ です。
理論値（解析解）の $${\hat{q}=0.457}$$ とほぼ一致しています。
ベイズ統計モデルで良好な推定ができたようです。

◼️MCMC サンプルのチェック
MCMC サンプルの軌跡＝トレースプロットを確認します。

# トレースプロット

# 描画領域の設定
plt.figure(figsize=(10, 3))
# トレースプロットの描画
plt.plot(samples, lw=0.5)
# 事後平均の水平線の描画
plt.axhline(q_mean, ls='--', color='tab:red', label=f'事後平均:{q_mean:.4}')
# 修飾
plt.title('パラメータ $q$')
plt.xlabel('MCMC Step数', fontsize=12)
plt.ylabel('生存確率 $q$', fontsize=12)
plt.legend();

【実行結果】
サンプルたちは満遍なくギザギザになっており、偏りや一定の傾向が見られません。
いい感じにサンプリングできたようです。

◼️パラメータ推定値を探る
パラメータ $${q}$$ の事後分布を可視化しましょう。

# 事後分布プロット

# 設定と準備
# x軸の値
q_vals = np.linspace(0.3, 0.65, 101)
# KDEの算出
kernel = stats.gaussian_kde(samples)
kde_vals = kernel(q_vals)
# KDE最大値のqの値
q_max = q_vals[kde_vals.argmax()]

# 描画
# qの事後分布のKDEプロット
plt.plot(q_vals, kde_vals, lw=3, label='MCMCサンプルのKDE')
# 共役事前分布の事後分布（ベータ分布）のプロット
plt.plot(q_vals, stats.beta.pdf(q_vals, a=a+Y, b=b+T-Y), color='tomato',
         label=f'共役Beta(a={int(a+Y)}, b={int(b+T-Y)})')
# qの事後平均値の垂直線（赤点線）の描画
plt.axvline(q_mean, ls='--', color='tab:red', label=f'事後平均　:{q_mean:.4f}')
# qの事後最大値の垂直線（緑点線）の描画
plt.axvline(q_max, ls=':', color='tab:green', label=f'事後最大値:{q_max:.4f}')
# 修飾
plt.xlabel('$q$', fontsize=12)
plt.ylabel('確率密度', fontsize=12)
plt.legend();

【実行結果】

青太線が MCMC サンプルの KDE 曲線、つまり、パラメータ推定値の分布です。
赤線がベータ事後分布の理論値による確率密度です。
２つの曲線はほぼ一致しています。
ベイズ統計モデルの MCMC サンプルで、パラメータ $${q}$$ の確率密度を程よく推定できているようです。

最後にパラメータ $${q}$$ の統計量をまとめます。

# MCMCサンプルからパラメータqの統計量を表示

# 中央値、95%区間の値の取得
quantiles = np.quantile(samples, [0.5, 0.025, 0.975])

# データフレーム化
pd.DataFrame({
	'mean': [samples.mean()],
	'std': [samples.std()],
	'median': [quantiles[0]],
	'95%CI lower': [quantiles[1]],
	'95%CI upper': [quantiles[2]],
	}, index=pd.Series(['q'], name='パラメータ')
).round(4)

【実行結果】

📣 最後にもう一言 📣

ベイズ統計モデルの MCMC 生成で用いたアルゴリズムは独自関数で実装しています。

小難しいことやってるなぁ

とお思いのあなた！
その気持ち、当然な感覚だと思います！

そして、ご安心下さい。
Python にはベイズ統計ライブラリが整っているので、MCMC アルゴリズムのコードを自分で書く必要はありません！

次回のベイズ統計モデリングからは、PyMC ライブラリでサクッとモデリングする予定です。
よろしくお願いいたします！

まとめ

長い道のりをここまで辿ってくださり、ありがとうございます！
次回以降のベイズ統計モデルにつながるようにまとめます。

ベイズ統計モデリングでは、マルコフ連鎖の定常分布≒パラメータの事後分布を目標にして MCMC アルゴリズムでサンプルを生成し、サンプルを分析してパラメータを推定します。

ベイズ統計モデルの重要な数式は：

です。

MCMC法を用いるベイズ統計モデリングは：

と言えそうです（大雑把です）。

🍀🍀🍀

それからぜひテキストをお読みくださいね！
8.5 節「補足説明」では、「詳細釣り合いの条件」を紐解いて、メトロポリス法から得られるサンプルが定常分布からのランダムサンプルになることを説明しています。

📊 詳細釣り合いの条件の公式（テキストの例の場合）

詳細釣り合い条件は、パラメータの事後分布をマルコフ連鎖の定常分布にするための十分条件だそうです。
テキストによると「メトロポリス法に限らず、全ての MCMC アルゴリズムに共通する性質」とのこと。

今回のブログは以上です。

次回は GLM（ポアソン回帰）のベイズ統計モデリングを学びます。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。