「データ解析のための統計モデリング入門」をPythonで写経 Vol.19 ～ ９章「GLMのベイズモデル化と事後分布の推定」②ベイズ統計モデリング～ギブスサンプリング
９章「GLMのベイズモデル化と事後分布の推定」
書籍の著者 久保拓弥 先生
書籍「データ解析のための統計モデリング入門」９章「GLMのベイズモデル化と事後分布の推定」の Python写経活動記録 です。
この記事は MCMCアルゴリズムの一種、ギブスサンプリング に取り組みます。
複数パラメータの MCMC サンプリング に関するロジックや動きを体感いたします。
アディショナルもあります！
では書籍を開いて統計モデリングの旅に出かけましょう🚀
はじめに
このブログシリーズは、書籍「データ解析のための統計モデリング入門 一般化線形モデル・階層ベイズモデル・MCMC」（岩波書店、「テキスト」と呼びます）の Python 写経を通じて得た「統計モデリングの楽しさ」をご紹介します。
テキストの紹介と引用表記はリンク先の記事に掲載しています。
準備
準備
■ 記事の範囲
この記事はテキスト９章の以下の節を取り扱います。
9.6 複数パラメーターのMCMCサンプリング
■ 利用データ
テキスト・サポートサイトのデータファイルを引用しています。
▶️ サポートサイト
Jupyter Notebook ファイルと同一フォルダ内に「data」フォルダを用意して、data フォルダ配下の章別フォルダにデータファイルを格納しています。
■ ライブラリのインポート
この記事で用いるライブラリをインポートします。
# インポート
# 数値計算
import numpy as np
import pandas as pd
# 統計計算
import scipy.stats as stats
# ベイズ統計
import arviz as az
# Rデータセットの読み込み
import rdata
# 可視化
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo' # または import japanize_matplotlib
データの準備
前回記事と同じ例題データを利用します。
d.RData ファイルを rdata ライブラリで変換して、pandas データフレームの data に読み込みます。
# 例題：種子数のポアソン回帰（個体差なし）のデータ p.194
# RDataファイルの読み込み
data = rdata.read_rda('./data/ch09/d.RData', default_encoding='ASCII')['d']
# 変数yを整数型に変換
data['y'] = data['y'].astype('int')
# csvファイルの出力
# data.to_csv('./data/ch09/data.csv')
# データフレームの表示
print('data.shape: ', data.shape)
data.head()【実行結果】
データの個数（標本サイズ）は 20 です。
20 個体の植物に関する仮想の観測データです。
【変数の説明】
植物の個体ごとの種子数 y（目的変数）と体サイズ x（説明変数）です。
$$
\begin{array}{clll}
変数 & 説明 & 値 \\
\hline
\\
y & 個体 i の種子数 & 0以上の実数 \\
x & 個体 i の体サイズ & 0以上の実数 \\
\end{array}
$$
統計モデリング・サマリー
■ 統計モデル
ポアソン回帰のベイズ統計モデルです。
前回記事と同じモデルです。
$$
\begin{array}{ll}
モデル & 特徴 \\
\hline
\\
ベイズ & ポアソン分布・対数リンク関数 \\
\end{array}
$$
以下のベイズ統計モデルを扱います。
$$
\begin{align*}y_i &\sim \text{Poisson}(\text{mu}=\lambda_i) \\
\lambda_i &= \exp(\beta_1 + \beta_2 (x_i - \bar{x})) \\
\beta_1 &\sim \text{Normal(\text{mu}=0, \text{sigma}=100)} \\\beta_2 &\sim \text{Normal(\text{mu}=0, \text{sigma}=100)} \\
\end{align*}
$$
$${\bar{x}}$$ は $${x}$$ の標本平均です。
$${x_i - \bar{x}}$$ で 体サイズ x の中心化を行っています。
■ モデリング手続き
ギブスサンプリングと呼ばれる MCMC アルゴリズムの理解に注力します。
ギブスサンプリング
複数パラメータのサンプリング
テキストは複数パラメータを 同時に更新する MCMC サンプリングは簡単ではない、としています。
その代わりに複数パラメータを 交互に更新する 方法を検討します。
$${\beta_1, \beta_2}$$ の２つのパラメータの場合は、次の 1、2 を交互に繰り返して更新するそうです。
1. （更新後の）$${\beta_2}$$ を定数とみなして $${\beta_1}$$ を更新
2. （更新後の）$${\beta_1}$$ を定数とみなして $${\beta_2}$$ を更新
部分的なサンプリングがＯＫな理由を、テキストは、
MCMC アルゴリズムによって定められる定常分布は、このような交互に更新する操作によっても変わらないためです
と説明しています。
ギブスサンプリングの登場
ギブスサンプリングは MCMC における値の更新に特徴があります。
新しい値の確率分布を作り、その確率分布からのランダムサンプルを新しい値にします。
新しい値の確率分布は、多変量確率分布から一つの変数を除いて、他の変数全てを定数とする一変量確率分布＝ 全条件付き分布 です。
パラメータに関する全条件付き分布を作ることが特徴です。
ギブスサンプリングの動き
テキスト p.215 ～ 217 の図 9.7 ～ 9.9 で示される MCMC の step 1 ～ 3 の動きチャートをお借りして、ギブスサンプリングの動きと数式を確認します。
一気にコード化します。
◼️ 作図処理
① 設定
# 設定と準備
# β1,β2の値
b1_vals = np.arange(1.4, 2.6, 0.002)
b2_vals = np.arange(-0.1, 0.1, 0.0005)
# (A)グラフのグラデーションに用いる範囲
gradation_vals = [0.95, 0.75, 0.55, 0.35, 0.15]
# (B)グラフの描画用パラメータ
plot_params = [
dict(xlim=((1.55, 2.45)), xlabel='$\\beta_1$'),
dict(xlim=((-0.025, 0.065)), xlabel='$\\beta_2$')
]② 関数定義：サポートサイトの配布コード「plot.gibbs.R」に準拠（引用）
# 関数の定義
# テキストのRコード plot.gibbs.R をpythonに変換しています
# (A)のグラデーション用のbの幅を算出する関数
def qb(vb, pb, range):
a = (1 - range) * 0.5
range_items = vb[(pb >= a) & (pb <= (1 - a))]
return range_items.min(), range_items.max()
# 条件付き事後確率密度に関する何かを算出する関数１（関数２の引数になる）
def get_prob(b1, b2): # b2*data.y はRコードのママ
return np.exp(
np.sum(stats.poisson.logpmf(k=data.y, mu=np.exp(b1 + b2*data.y)))
)
# 条件付き事後確率密度に関する何かを算出する関数２
def d2p(vd):
cs = np.cumsum(vd)
return dict(d=vd, p=cs / np.max(cs))
# β1のときに、条件付き事後確率密度に関する何かを算出する関数を呼び出す関数
def get_prob_b1(b1_vals, b2):
return d2p([get_prob(b1, b2) for b1 in b1_vals])
# β2のときに、条件付き事後確率密度に関する何かを算出する関数を呼び出す関数
def get_prob_b2(b1, b2_vals):
return d2p([get_prob(b1, b2) for b2 in b2_vals])
# β1のときに(A)グラフのグラデーションを描画する関数
def draw_envelope_b1(b1s, bb2, ax):
draw_envelope(b1s[0], bb2, b1s[1], bb2, ax)
# β2のときに(A)グラフのグラデーションを描画する関数
def draw_envelope_b2(b1, b2s, ax):
draw_envelope(b1, b2s[0], b1, b2s[1], ax)
# (A)グラフのグラデーションの描画本体関数
def draw_envelope(b1L, b2L, b1U, b2U, ax):
ax.fill_between(
data.x,
np.exp(b1L + b2L * data.x),
np.exp(b1U + b2U * data.x),
color='tomato',
alpha=0.1
)③ 関数定義：自作関係
# 関数の定義 オリジナル関数
# (B)の事条件付き後確率密度の曲線を描画する関数
def draw_density(x_vals, density, b_new, params, ax):
# 確率密度の曲線の描画
ax.plot(x_vals, density)
# サンプリングされた値の垂直点線の描画
ax.vlines(b_new, 0, density.max()*1.1, color='gray', ls='--')
# 修飾
ax.set(xlim=params['xlim'], ylim=(0), xlabel=params['xlabel'])
# 枠線の消去
ax.spines['left'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
# 目盛り・目盛りラベルの消去
ax.tick_params(left=False, labelleft=False)
# β1の変化のグラフを描画する関数
def plot_gibbs_b1(b1, b2, b1_new, step):
# 描画用データの算出 ※テキストのRコードをpythonに変換
pb1 = get_prob_b1(b1_vals, b2)
post_density1 = pb1['d'] / sum(pb1['d']) # 条件付き事後確率密度として描画する
# 描画領域の設定
fig, ax = plt.subplots(1, 3, figsize=(10, 4), tight_layout=True)
# (A)の散布図の描画
sns.scatterplot(data=data, x='x', y='y', ec='tab:blue', fc='white', s=70,
color='tab:blue', ax=ax[0])
# (A)のβ1の条件付き事後確率密度のグラデーションの描画
for g_val in gradation_vals:
draw_envelope_b1(qb(b1_vals, pb1['p'], g_val), b2, ax[0])
# (A)の修飾
ax[0].set(title=f'(A)$\\beta_2$={b2:.3f}\nという条件のもとで',
xlabel='植物の体サイズ $x$', ylabel='種子数 $y$')
# (B)の条件付き事後確率密度の曲線、サンプリングされた新しいβ1の垂直点線の描画
draw_density(b1_vals, post_density1, b1_new, plot_params[0], ax[1])
# (B)の修飾
ax[1].set(title=f'(B)$\\beta_1^{{new}}$={b1_new:.3f}\nがサンプリングされた')
# (C)の散布図の描画
sns.scatterplot(data=data, x='x', y='y', ec='tab:blue', fc='white', s=70,
color='tab:blue', ax=ax[2])
# (C)の変化前の平均λの点線の描画
ax[2].plot(data.x, np.exp(b1 + b2 * data.x), color='gray', ls='--')
# (C)の変化後の平均λの実線の描画
ax[2].plot(data.x, np.exp(b1_new + b2 * data.x), color='tab:blue', lw=2)
# (C)の修飾
ax[2].set(title=f'(C) 平均値の変化\n$\\lambda=exp(\\beta_1+\\beta_2x)$',
xlabel='植物の体サイズ $x$', ylabel='種子数 $y$')
# ステップの表示
fig.suptitle(f'STEP. {step} → $\\beta_1$', x=0.1, weight='bold', fontsize=16)
# β2の変化のグラフを描画する関数
def plot_gibbs_b2(b1, b2, b2_new, step):
# 描画用データの算出 ※テキストのRコードをpythonに変換
pb2 = get_prob_b2(b1, b2_vals)
post_density2 = pb2['d'] / sum(pb2['d']) # 条件付き事後確率密度として描画する
# 描画領域の設定
fig, ax = plt.subplots(1, 3, figsize=(10, 4), tight_layout=True)
# (A)の散布図の描画
sns.scatterplot(data=data, x='x', y='y', ec='tab:blue', fc='white', s=70,
color='tab:blue', ax=ax[0])
# (A)のβ2の条件付き事後確率密度のグラデーションの描画
for g_val in gradation_vals:
draw_envelope_b2(b1, qb(b2_vals, pb2['p'], g_val), ax[0])
# (A)の修飾
ax[0].set(title=f'(A)$\\beta_1$={b1:.3f}\nという条件のもとで',
xlabel='植物の体サイズ $x$', ylabel='種子数 $y$')
# (B)の条件付き事後確率密度の曲線、サンプリングされた新しいβ2の垂直点線の描画
draw_density(b2_vals, post_density2, b2_new, plot_params[1], ax[1])
# (B)の修飾
ax[1].set(title=f'(B)$\\beta_2^{{new}}$={b2_new:.3f}\nがサンプリングされた')
# (C)の散布図の描画
sns.scatterplot(data=data, x='x', y='y', ec='tab:blue', fc='white', s=70,
color='tab:blue', ax=ax[2])
# (C)の変化前の平均λの点線の描画
ax[2].plot(data.x, np.exp(b1 + b2 * data.x), color='gray', ls='--')
# (C)の変化後の平均λの実線の描画
ax[2].plot(data.x, np.exp(b1 + b2_new * data.x), color='tab:blue', lw=2)
# (C)の修飾
ax[2].set(title=f'(C) 平均値の変化\n$\\lambda=exp(\\beta_1+\\beta_2x)$',
xlabel='植物の体サイズ $x$', ylabel='種子数 $y$')
# ステップの表示
fig.suptitle(f'STEP. {step} → $\\beta_2$', x=0.1, weight='bold', fontsize=16)いや～、長いコードですね…
では、図 9.7、9.8、9.9 に相当するチャートを一気に作成するコードを実行します。
# β1とβ2のギブスサンプリングの可視化 p.215~217 図9.7, 9.8, 9.9
# 設定：β1,β2の値 [初期値, new1, new2, new3]
b1_samples = [1.5, 2.052, 2.017, 1.901]
b2_samples = [0.0, -0.016, 0.015, 0.017]
# step.3までβ1とβ2を交互にグラフ(A)～(C)の描画を繰り返し処理
for i in range(6):
# 位置の取得
q, mod = divmod(i, 2)
# iが偶数の時、β1の変化ステップを描画
if mod == 0:
plot_gibbs_b1(b1_samples[q], b2_samples[q], b1_samples[q+1], q+1)
# iが奇数の時、β2の変化ステップを描画
else:
plot_gibbs_b2(b1_samples[q+1], b2_samples[q], b2_samples[q+1], q+1)【実行結果】
図全体が描画されました。
図をステップごとの説明で利用したいので、ステップごとに分割して掲載します。
🍀🍀🍀
◼️ ギブスサンプリングの動き
① STEP. 1：最初の $${\beta_1}$$ のサンプリング
$${\beta_2=0}$$ とし、全条件付き分布（FCD ）$${p(\beta_1 \mid \bm Y, \beta_2)}$$ は次の式で表されます。
$$
\begin{align*}
p(\beta_1 \mid \bm Y, \beta_2 = 0.0) &\propto \prod_i^N \cfrac{\lambda_i^{y_i} \exp(\lambda_i)}{y_i!}\ p(\beta_1) \\
\lambda_i &= \exp(\beta_1 + 0)
\end{align*}
$$
FCD から $${\beta_1^{新} = 2.052}$$ を得ました。
② STEP. 1：最初の $${\beta_2}$$ のサンプリング
$${\beta_1=2.052}$$ とし、全条件付き分布（FCD ）$${p(\beta_2 \mid \bm Y, \beta_1)}$$ は次の式で表されます。
$$
\begin{align*}
p(\beta_2 \mid \bm Y, \beta_1 = 2.052) &\propto \prod_i^N \cfrac{\lambda_i^{y_i} \exp(\lambda_i)}{y_i!}\ p(\beta_2) \\
\lambda_i &= \exp(2.052 + \beta_2 x_i)
\end{align*}
$$
FCD から $${\beta_2^{新} = -0.016}$$ を得ました。
③ STEP. 2：$${\beta_1}$$ のサンプリング
$${\beta_2=-0.016}$$ とし、全条件付き分布（FCD ）$${p(\beta_1 \mid \bm Y, \beta_2)}$$ は次の式で表されます。
$$
\begin{align*}
p(\beta_1 \mid \bm Y, \beta_2 = -0.016) &\propto \prod_i^N \cfrac{\lambda_i^{y_i} \exp(\lambda_i)}{y_i!}\ p(\beta_1) \\
\lambda_i &= \exp(\beta_1 -0.016 x_i)
\end{align*}
$$
FCD から $${\beta_1^{新} = 2.017}$$ を得ました。
④ STEP. 2：$${\beta_2}$$ のサンプリング
$${\beta_1=2.017}$$ とし、全条件付き分布（FCD ）$${p(\beta_2 \mid \bm Y, \beta_1)}$$ は次の式で表されます。
$$
\begin{align*}
p(\beta_2 \mid \bm Y, \beta_1 = 2.017) &\propto \prod_i^N \cfrac{\lambda_i^{y_i} \exp(\lambda_i)}{y_i!}\ p(\beta_2) \\
\lambda_i &= \exp(2.017 + \beta_2 x_i)
\end{align*}
$$
FCD から $${\beta_2^{新} = 0.015}$$ を得ました。
⑤ STEP. 3：$${\beta_1}$$ のサンプリング
$${\beta_2=0.015}$$ とし、全条件付き分布（FCD ）$${p(\beta_1 \mid \bm Y, \beta_2)}$$ は次の式で表されます。
$$
\begin{align*}
p(\beta_1 \mid \bm Y, \beta_2 = 0.015) &\propto \prod_i^N \cfrac{\lambda_i^{y_i} \exp(\lambda_i)}{y_i!}\ p(\beta_1) \\
\lambda_i &= \exp(\beta_1 + 0.015 x_i)
\end{align*}
$$
FCD から $${\beta_1^{新} = 1.901}$$ を得ました。
⑥ STEP. 3：$${\beta_2}$$ のサンプリング
$${\beta_1=1.901}$$ とし、全条件付き分布（FCD ）$${p(\beta_2 \mid \bm Y, \beta_1)}$$ は次の式で表されます。
$$
\begin{align*}
p(\beta_2 \mid \bm Y, \beta_1 = 1.901) &\propto \prod_i^N \cfrac{\lambda_i^{y_i} \exp(\lambda_i)}{y_i!}\ p(\beta_2) \\
\lambda_i &= \exp(1.901 + \beta_2 x_i)
\end{align*}
$$
FCD から $${\beta_2^{新} = 0.017}$$ を得ました。
こんな感じでギブスサンプリングが「交互に」進んでいきます。
ギブスサンプリングのメリットをテキストが挙げています。
・各 MCMC ステップで更新前の値と更新後の値の相関が小さい
・MCMC サンプリングの詳細を指定しなくてよい
ちなみに、現在のベイズ統計ライブラリは別の MCMC アルゴリズムを活用しています。
PyMC や Stan はデフォルトで NUTS と呼ばれる ハミルトニアンモンテカルロ（HMC）法をベースにしたアルゴリズムを使っています。
（注意）
アディショナルは趣味的な深堀りとコードです。
ご興味ない方はスルーしてくださって大丈夫です。
アディショナル：ギブス的サンプリング
今回モデルとギブスサンプリングの相性
◼️ ギブスサンプリングが使える条件
テキスト p.216 以降で 「全条件付き確率（FCD）に従う乱数の発生方法」を検討しています。
パラメータの事前分布と事後分布（または FCD）が同じ種類の確率分布になる「共役事前分布」のときにギブスサンプリングを活用できるそうです。
詳細はぜひテキストをお読み下さい！
ところで残念なお知らせがあります…
今回の「ポアソン分布のGLMで線形予測子のパラメータ $${\beta_1, \beta_2}$$ が正規分布に従うモデル」は、残念ですが、ギブスサンプリングに適さないようです…
🍀🍀🍀
◼️ ギブスサンプリングが使える条件の具体例
ChatGPTに「ギブスサンプリングが使えるパラメータと共役事前分布の関係を整理して教えて下さい」と軽く質問したところ、大量の回答が出てきました！
なお回答内容の適否を確認をしていないので、この点を踏まえてお読み下さい。
1️⃣ まずおさえる：指数型分布族と共役事前の一般形
$$
p(y \mid \eta) = h(y)\ \exp{\eta^\top T(y)-A(\eta)}
$$
のとき、共役事前は
$$
p(\eta\mid \alpha,\beta) \propto \exp \{\eta^\top \alpha-\beta A(\eta) \}
$$
となる（同じ族を保つ）。データ観測後、ハイパーパラメータは
$$
\alpha' = \alpha + \sum_i T(y_i), \quad \beta' = \beta+n
$$
のように足し込みで更新（= ギブス向き）。
2️⃣ 代表ペア（単層モデル）
3️⃣ 線形回帰・階層化（ガウス系）
4️⃣ GLM とデータ拡張で「条件付き共役」にする例
5️⃣ 混合・クラスタリング系（ギブス鉄板セット）
ギブスサンプリングは多くのケースで使える「頼りになるヤツ」です！
Metropolis-within-Gibbs の登場
しかしながら、今回のモデルはギブスサンプリング適用が難しいのです。
そこで ChatGPT に今回モデルでギブスサンプリング的な MCMC サンプリングはできないだろうか、と相談したところ…
Metropolis-within-Gibbs の提案をもらいました！
ギブスサンプリングに内包されるメトロポリス・ヘイスティング（MH）法ですって！
詳細を ChatGPT に訊きました。
🔷 Metropolis-within-Gibbs（MWG）ってなに？
要するに：
“ギブスサンプリングの枠組み”でパラメータを順番に更新しつつ、その場で閉形式の条件付き分布から引けないパラメータだけを Metropolis–Hastings（MH）で更新する方法です。
→ 共役でないモデル（GLMや階層モデルなど）でも、ギブスの手軽さを保ったまま回せます。
🔷 どうして必要？
純粋なギブスは「各パラメータの完全条件付き分布が“既知の分布”として書ける（=共役）」ことが前提。
でも実務の多くは共役じゃない → そこで そのパラメータだけ メトロポリス・ヘイスティング法（MH法） で一歩進める＝MWG。
🔷 仕組み（1ステップ）
あるパラメータ $${\theta_j}$$ 以外を固定しておく。
$${\theta_j}$$ の候補 $${\theta_j'}$$ を提案分布（例：現在値まわりの正規）からサンプル。
次のパラメータへ。これを全パラメータで1巡=1イテレーション。
提案分布を対称（例：正規のランダムウォーク）にすれば $${q}$$ は打ち消し合い、「事後比」だけでOK。
🔷 対応表（ここがギブス／ここがMH法）
ギブスの骨格（Gibbs skeleton）
パラメータを 1つずつ順番に更新する（他は固定）。
そのときの対象は 条件付き分布 $${p(\theta_j \mid \theta_{-j}, y)}$$。
→ つまり「順番に、条件づけて」回す部分がギブス。
MH法で置き換える部分（Metropolis “within”）
条件付きから直接サンプルできないとき、
3. 提案 $${\theta_j' \sim q(,\cdot\mid \theta_j)}$$ を出す。
4. 受容確率
$$
\alpha=\min \left(1,\
\frac{p(\theta_j'\mid \theta_{-j},y)\ q(\theta_j\mid \theta_j')}
{p(\theta_j \mid \theta_{-j},y)\ q(\theta_j'\mid \theta_j)}
\right)
$$
で採択/棄却する。
→ この 提案＋受容判定 が MH法 の担当。
注：もし $${p(\theta_j \mid \theta_{-j}, y)}$$ がベータや正規など 閉形式で直接引けるなら、そこは純ギブスに置き換えられます（MWGは「引ける所はギブス、引けない所だけMH法」）。
🔷 擬似コード（ラベル付き）
初期値 θ ← θ(0)
for t = 1..T: # [Gibbs骨格] 反復ループ
for j in 1..p: # [Gibbs骨格] パラメータを順番に
# --- ここから、条件付きの取り扱い ---
# Case A: 条件付きが閉形式で引けるなら（純ギブス）
# θ_j ~ p(θ_j | θ_-j, y) # [Gibbsステップ: 直接サンプル]
# Case B: 引けないなら（MHで置換）
提案 θ_j' ~ q(· | θ_j) # [MH: 提案]
r = log p(θ_j' | θ_-j, y) - log p(θ_j | θ_-j, y)
+ log q(θ_j | θ_j') - log q(θ_j' | θ_j) # [MH: 受容比(対数)]
受容: log U < r なら θ_j ← θ_j' # [MH: 受容/棄却]（U~Uniform(0,1)）
θ(t) を保存対称提案（例：ランダムウォーク正規 $${q(\theta_j'|\theta_j)=\mathcal N(\theta_j,s_j^2)}$$）なら、$${q}$$ が打ち消し合い、事後の差だけでOK：
$$
\alpha=\min\bigl(1,\ \exp{\log p(\theta_j'|\theta_{-j},y)-\log p(\theta_j|\theta_{-j},y)}\bigr).
$$
交互に更新する部分がギブスサンプリング的。
新しい値の提案と受容判定が MH 法的。
ギブスサンプリングと MH 法のハイブリッドが Metropolis-within-Gibbs なのです！
Metropolis-within-Gibbs の実装
いよいよクライマックス！
今回のベイズ統計モデルで Metropolis-within-Gibbs による MCMC サンプルを生成します。
コード実装も ChatGPT に全面依拠いたします…
🍀🍀🍀
① Metropolis-within-Gibbs 実行関数 by ChatGPT
ChatGPTにコメントをたくさんつけてもらいました。
【プロンプト】
コードには「やさしくて詳しくて網羅性ある日本語のコメント」を付けて下さい！
ぜひコメントをお読みいただき、コードの意味合いを感じてください。
# ============================================
# ポアソン回帰（ベイズ）を
# Metropolis-within-Gibbs で推定する最小実装
# --------------------------------------------
# モデル:
# y_i ~ Poisson(λ_i)
# λ_i = exp(β1 + β2 * (x_i - x̄)) # 説明変数は中心化
# β1 ~ Normal(0, 100)
# β2 ~ Normal(0, 100)
#
# ポイント：
# - このモデルは「完全条件付き分布」が閉形式にならないため、
# “純粋なギブスサンプリング”は組めません。
# そこで「各パラメータをランダムウォークMHで順番に更新する」
# Metropolis-within-Gibbs（MWG）を使います。
# - 1イテレーション内で
# β1 を MH で更新 → β2 を MH で更新
# と回せば、それはギブス骨格の MCMC（MWG）です。
# - 受容率（acceptance rate）＝「提案を採択した割合」も
# 併せて計算します。目安は 0.3〜0.5 程度。
# ============================================
# import numpy as np
# ---- 事後対数密度 log p(β | y, x) を返す関数 -------------------------
def _logpost(beta, y, xc, sd_prior=100.0):
'''
事後対数密度（定数分は落としてOK）を返す補助関数。
Parameters
----------
beta : (2,) array-like
[β1, β2]
y : (N,) array-like
観測カウント
xc : (N,) array-like
中心化済みの説明変数（x - x.mean()）
sd_prior : float
事前分布 N(0, sd_prior^2) の標準偏差
Returns
-------
float
事後対数密度（定数ぶんを落としているので相対値）
'''
# パラメータの取り出し
b1, b2 = beta
# 線形予測子の算出：η = β1 + β2 * x_c
eta = b1 + b2 * xc
# 数値安定化のため、expの過大化を軽く抑制（任意）
# ※ 過度なクリップは推定を歪めるので、必要なければ外してください。
lam = np.exp(np.clip(eta, -40, 40)) # λ = exp(η)
# ポアソンの対数尤度の算出： sum( y*log λ - λ ) = sum( y*η - λ )
# （-log(y!)はβに依らないので捨象）
ll = np.dot(y, eta) - lam.sum()
# 事前分布の対数確率密度の和の算出：β1, β2 ~ N(0, sd_prior^2)
# → 対数は -0.5 * (β/sd)^2 の和
lp = -0.5 * ((b1 / sd_prior)**2 + (b2 / sd_prior)**2)
return ll + lp
# ---- Metropolis-within-Gibbs 本体 ------------------------------------
def sample_mwg(
y,
x,
draws=5000,
tune=1000,
thin=1,
prop_sd=(0.1, 0.1),
seed=0,
sd_prior=100.0,
):
'''
ポアソン回帰（正規事前）を Metropolis-within-Gibbs で推定。
Parameters
----------
y : (N,) array-like
観測カウント
x : (N,) array-like
説明変数
draws : int
保存するサンプル数（バーンイン後・間引き適用後）
tune : int
バーンイン反復回数（保存しない）
thin : int
間引き間隔。thin=5 なら5ステップごとに1回保存。
prop_sd : (float, float)
それぞれ (β1, β2) のランダムウォーク提案N(0, sd^2)の標準偏差。
受容率が 0.3〜0.5 くらいになるよう調整すると安定。
seed : int
乱数シード（再現性のため）
sd_prior : float
事前 N(0, sd_prior^2) の標準偏差
Returns
-------
dict
- 'beta1', 'beta2' : 事後サンプル（draws個）
- 'acc_rate_all' : バーンイン込みの受容率（各β）
- 'acc_rate_post' : バーンイン後のみの受容率（各β）
- 'x_centered' : 中心化後の説明変数（確認用）
'''
## 設定と準備
# 乱数生成器の作成
rng = np.random.default_rng(seed)
# データをnumpy配列化
x = np.asarray(x, dtype=float)
y = np.asarray(y, dtype=float)
# x を中心化（式どおり）
# 切片と勾配の推定が安定し、β1の解釈も「平均的xでのlog λ」に。
xc = x - x.mean()
# MCMCサンプルを格納するリスト
out_b1, out_b2 = [], []
# 提案に使用する正規分布の標準偏差の取り出し
sd1, sd2 = prop_sd
## MCMC
# パラメータ初期値の設定：
# - β1 は log(mean(y)) を初期値にすると入りやすい（λ の初期が平均yに近い）
# - β2 は 0 から開始
b1 = np.log(max(y.mean(), 1e-12))
b2 = 0.0
# 現在点の対数事後密度の算出
cur_lp = _logpost((b1, b2), y, xc, sd_prior)
# 受容カウンタの初期化（「提案回数 = 反復回数」なので分母は total / (post)）
acc1_all = acc2_all = 0
acc1_post = acc2_post = 0
## MCMC の実行
total = tune + draws * thin # これだけ回して、thin間隔で保存する
for t in range(total):
# Step.1) β1 を MH で更新（β2 は固定のまま）---
# ランダムウォーク提案 β1' = β1 + ε, ε~N(0, sd1^2)
cand = b1 + rng.normal(0.0, sd1)
lp_c = _logpost((cand, b2), y, xc, sd_prior)
# 受容判定：log u < log posterior 差
if np.log(rng.uniform()) < (lp_c - cur_lp):
b1 = cand
cur_lp = lp_c
acc1_all += 1
if t >= tune:
acc1_post += 1
# Step.2) β2 を MH で更新（β1 を固定して）---
# ランダムウォーク提案 β2' = β2 + ε, ε~N(0, sd2^2)
cand = b2 + rng.normal(0.0, sd2)
lp_c = _logpost((b1, cand), y, xc, sd_prior)
# 受容判定：log u < log posterior 差
if np.log(rng.uniform()) < (lp_c - cur_lp):
b2 = cand
cur_lp = lp_c
acc2_all += 1
if t >= tune:
acc2_post += 1
# Step.3) 保存（tune 以降、thin間隔ごとに）---
if t >= tune and ((t - tune) % thin == 0):
out_b1.append(b1)
out_b2.append(b2)
# 受容率を集計（各パラメータ、反復ごとに1回提案している）
n_all = total
n_post = total - tune
acc_all = {'beta1': acc1_all / n_all, 'beta2': acc2_all / n_all}
acc_post = {'beta1': acc1_post / n_post, 'beta2': acc2_post / n_post}
return {
'beta1': np.asarray(out_b1),
'beta2': np.asarray(out_b2),
'acc_rate_all': acc_all,
'acc_rate_post': acc_post,
'x_centered': xc,
}
# --- ユーティリティ 簡易なESS, MCSEを算出するヘルパー関数
def acf(x, max_lag=100):
x = np.asarray(x); x = x - x.mean()
n = len(x); var = np.dot(x, x)/n
ac = np.empty(max_lag+1); ac[0] = 1.0
for k in range(1, max_lag+1):
ac[k] = np.dot(x[:-k], x[k:])/(n-k)/var
return ac
def ess(x, max_lag=100):
r = acf(x, max_lag)
# 0 をまたいだところで打ち切る素朴版（十分実用的）
pos = r[1:]
cutoff = np.argmax(pos < 0) + 1 if np.any(pos < 0) else max_lag
tau = 1 + 2*np.sum(r[1:cutoff])
n = len(x)
return n / tau
def mcse_mean(x, max_lag=100):
return np.std(x, ddof=1) / np.sqrt(ess(x, max_lag))
# --- 事後予測（PPC）と軽い整合チェック
# 0.05 未満や0.95 超など極端なら、その統計量に関するモデルとデータの不整合のサイン
def ppc_and_checks(res, y, rng_seed=0):
# 準備
# 乱数生成器
rng = np.random.default_rng(rng_seed)
# 事後予測用の変数の取得
b1, b2 = res['beta1'], res['beta2']
xc = res['x_centered']
# 平均λの事後予測の算出
lam = np.exp(b1[:, None] + b2[:, None] * xc[None, :]) # (draws, N)
# 観測yの事後予測の算出
y_ppc = rng.poisson(lam)
# 基本統計の整合
mean_obs = y.mean()
mean_ppc = y_ppc.mean()
var_ratio = y_ppc.var() / y.var(ddof=1)
# 両側PPC p値（統計量：最大値・ゼロ個数）
def ppc_two_sided(sim_stat, obs_stat):
ge = (sim_stat >= obs_stat).mean()
le = (sim_stat <= obs_stat).mean()
return min(1.0, 2 * min(ge, le))
# 最大値の両側PPC p値の算出
p_max = ppc_two_sided(y_ppc.max(axis=1), y.max())
# ゼロ個数の両側PPC p値の算出
p_zero = ppc_two_sided((y_ppc == 0).sum(axis=1), (y == 0).sum())
return {
'mean_obs': mean_obs, # 観測値の平均
'mean_ppc': mean_ppc, # yの事後予測の平均
'var_ratio_ppc_over_obs': var_ratio, # 分散比（予測/観測）
'p_max': p_max, # 最大値の両側PPCのp値
'p_zero': p_zero, # ゼロ個数の両側PPCのp値
'y_ppc': y_ppc, # yの事後予測：チャート作成で使える
}🍀🍀🍀
② Metropolis-within-Gibbsの実行
MCMC サンプルを取得して、ざっくり概観を確認します。
chain 数は１、draw 数は 10000 です。
# Metropolis-within-Gibbsの実行
## データの準備
# データセットの作成
x = data.x.values.copy()
y = data.y.values.copy()
# xの中心化データの作成
xc = x - x.mean()
## MCMCの実行
result_mwg = sample_mwg(
y=y, x=x, draws=10000, tune=1000, thin=1, prop_sd=(0.25, 0.2),
seed=123, sd_prior=100.0,
)
## 実行結果の概要表示
# 事後平均の表示
b1_mean = result_mwg['beta1'].mean()
b2_mean = result_mwg['beta2'].mean()
print('posterior mean (β1, β2) =', (b1_mean, b2_mean))
# 受容率の表示：0.3〜0.5 くらいを目標に。prop_sdで調整
print('acceptance (all) =', result_mwg['acc_rate_all']) # バーンイン込み
print('acceptance (post) =', result_mwg['acc_rate_post']) # バーンイン後のみ【実行結果】
あっという間に処理完了です。
事後平均（posterior mean）は前回記事の PyMC・NUTS による推定値に近いです！（次の③に前回記事の推定値を掲載）
受容率（acceptance）は新しい値提案を受容した割合です。
ChatGPT によると今回モデルの場合 0.3 ～ 0.5 が目安だそう。
🍀🍀🍀
③ 事後分布の要約統計量の表示
PyMC の summary() みたいな表を作ります。
ただし ESS と MCSE は簡易的に計算しています（ChatGPT 談）。
# 事後分布の要約統計量の表示
def post_summary(result, hdi_prob=0.95):
# HDI区間の算出
hdi_b1 = az.hdi(result['beta1'], hdi_prob=hdi_prob)
hdi_b2 = az.hdi(result['beta2'], hdi_prob=hdi_prob)
# 戻り値：MCMCサンプルの要約統計量をデータフレーム化
return pd.DataFrame(
{
'mean': [result['beta1'].mean(), result['beta2'].mean()],
'sd': [result['beta1'].std(), result['beta2'].std()],
'median': [np.median(result['beta1']), np.median(result['beta2'])],
f'hdi_{(1 - hdi_prob)/2:.2%}': [hdi_b1[0], hdi_b2[0]],
f'hdi_{1 - (1 - hdi_prob)/2:.2%}': [hdi_b1[1], hdi_b2[1]],
'ESS（簡易）' :
[ess(result_mwg['beta1']), ess(result_mwg['beta2'])],
'MCSE（簡易）':
[mcse_mean(result_mwg['beta1']), mcse_mean(result_mwg['beta2'])],
}, index=['beta1', 'beta2'],
)
# 事後分布の要約統計量の表示
post_summary(result_mwg).round(4)【実行結果】
PyMC による事後分布の要約統計量の各値によく似ています！
（参考：前回記事のPyMC・NUTSによる推定値）
ここからは事後分布の概要を見つつ、MCMC サンプルの良し悪しを確認します。
🍀🍀🍀
③ 事後分布とトレースプロットの可視化
# トレースプロットの描画
# 設定と準備
# MCMCサンプルの取り出し
samples = [result_mwg['beta1'], result_mwg['beta2']]
# チャートに表示するパラメータ名
var_names = [r'$\beta_1$', r'$\beta_2$']
# 描画領域の設定
fig, ax = plt.subplots(2, 2, figsize=(10, 6), tight_layout=True)
# パラメータごとに事後分布のKDEとトレースプロットを描画
for i, (sample, var_name) in enumerate(zip(samples, var_names)):
# 事後分布のKDEの描画（左）
sns.kdeplot(sample, ax=ax[i, 0])
ax[i, 0].set(title=f'{var_name} の事後分布')
ax[i, 0].set_xlabel(var_name)
# トレースプロットの描画（右）
ax[i, 1].plot(sample, lw=0.5)
ax[i, 0].set_ylabel('確率密度')
ax[i, 1].set(title=f'{var_name} のトレースプロット')
ax[i, 1].set_ylabel(var_name)【実行結果】
左：ほぼ単峰で正規分布に似ています。
右：満遍なくゲジゲジとしています。
両方がいい感じです。
🍀🍀🍀
④ 自己相関プロット
arviz の plot_autocorr() で各パラメータのMCMC サンプルの自己相関を確認します。
# 自己相関プロット
# β1, β2のMCMCサンプルの取り出し（辞書形式）
result_betas = {k: v for k, v in result_mwg.items() if k in ['beta1', 'beta2']}
# 自己相関プロットの描画
az.plot_autocorr(result_betas, figsize=(10, 3));【実行結果】
ChatGPT 的には「いい感じ」とのこと。
🍀🍀🍀
⑤ ２つのパラメータのMCMC サンプルの散布図
$${\beta_1}$$ のサンプルと $${\beta_2}$$ のサンプルのサンプルに相関がないことを確認します。
# MCMC サンプルの散布図 ※ほぼ相関なし
# MCMC サンプルの取り出し（各変数に代入）
b1_samples, b2_samples = result_betas['beta1'], result_betas['beta2']
# 散布図の描画
sns.scatterplot(x=b1_samples, y=b2_samples, alpha=0.3)
# 修飾
plt.xlabel('$\\beta_1$', fontsize=12)
plt.ylabel('$\\beta_2$', fontsize=12)
plt.title(f'$\\beta_1$ と $\\beta_2$ の相関係数：'
f'{np.corrcoef(b1_samples, b2_samples)[0, 1]:.3f}');【実行結果】
相関はないです。
🍀🍀🍀
⑥ パラメータ推定値を解釈
$${\exp(\hat{\beta_2})}$$ で解釈します。
# 効果の解釈用：IRR = exp(β2)
# β2のMCMCサンプルを観測値yの単位に戻す（log⇒exp）
irr = np.exp(result_mwg['beta2'])
# 平均値の算出
irr_mean = irr.mean()
# 95%信用区間の算出
irr_lo, irr_hi = np.quantile(irr, [0.025, 0.975])
# 結果の表示
print(f'\nIRR (=exp(β2)) 平均={irr_mean:.3f}, '
f'95%HDI=[{irr_lo:.3f}, {irr_hi:.3f}]')【実行結果】
体サイズ $${x_i}$$ が１単位増加したときの種子数 $${y_i}$$ の増減効果です。
95% HDI はゼロを挟んでいて、0.95 倍に減少する～ 1.2 倍に増加する、といった感じ。
良さげでは無さそうです。
🍀🍀🍀
⑦ 種子数 y の観測値と予測値の整合確認
ChatGPT に言われたまま、実行します。
# 予測整合（簡易PPC）
ppc = ppc_and_checks(result_mwg, y)
print('[PPC 簡易チェック]')
print(f"観測平均 = {ppc['mean_obs']:.3f}, 事後予測平均 = {ppc['mean_ppc']:.3f}")
print(f"分散比(予測/観測) = {ppc['var_ratio_ppc_over_obs']:.2f}")
print(f"PPC p値（両側）：max(y) = {ppc['p_max']:.3f}, "
f"count(y==0) = {ppc['p_zero']:.3f}")【実行結果】
【評価】
平均の比較：中心傾向の整合は良好
分散比：モデルが実データよりも“ばらつきを大きく見積もっている”状態
PPC両側 p値（最大値）：極端さは見られず
PPC両側 p値（ゼロの個数）：少なくともゼロの出現頻度は合っている示唆
🍀🍀🍀
⑧ 種子数 y の観測値の分布と事後予測分布の比較
# yの観測値vs事後予測の重ね描き
# yの観測値のKDE曲線付きヒストグラムの描画
sns.histplot(
x=y, bins=range(y.max()+2), stat='density', kde=True,
edgecolor='white', alpha=0.6, label='observed')
# yの事後予測のKDE曲線付きヒストグラムの描画
sns.histplot(
ppc['y_ppc'].ravel(), bins=range(y.max()+2), stat='density',
kde=True, kde_kws={'bw_adjust': 3},
color='tab:red', edgecolor='white', alpha=0.3, label='posterior pred')
# 修飾
plt.legend()
plt.xlabel('y')
plt.ylabel('density');【実行結果】
KDE 曲線は似ています。
事後予測分布（赤）は観測値の分布（青）から大外れしていない状況。
以上が用意したコードのすべてです。
MCMC サンプリングの実行時間はほんの一瞬でしたが、ベイズの奥深い分析・評価の一端に触れることができました。
面白かったですね！
ChatGPT さん、いつもありがとう🍀
まとめ
◼️ 振り返り
MCMC アルゴリズムの一種、ギブスサンプリングをテキストに沿って学びました。
ただし、今回のポアソン回帰のベイズ統計モデルは、パラメータ $${\beta_1,\beta_2}$$ の完全条件付き分布（FCD）が共役事前分布に対応しない（非共役）ため、FCDから直接乱数を生成するのが難しく、ギブスサンプリングが適さないことが分かりました。
代わりに、ChatGPT を活用してギブスサンプリングとメトロポリス・ヘイスティング法のハイブリッド「Metropolis-within-Gibbs」を動かしました。
ところで、PyMC や Stan などの最近のベイズ統計モデリング・ライブラリは、 ハミルトニアン・モンテカルロ（HMC）法の一種、NUTS（No-U-Turn Sampler）という MCMC アルゴリズムをデフォルトで用いるようです。
MCMC アルゴリズム、たくさん出てきて、うわー！ってなりますね苦笑い
🍀🍀🍀
◼️ MCMC アルゴリズムの学びの続きを
そんなときには、こちらの資料が助けてくれるでしょう。
分寺杏介 先生によるベイズ授業教材です。
３つの MCMC アルゴリズムを楽しく学べる最高のドキュメントです！
メトロポリス・ヘイスティングス法
ギブスサンプリング
ハミルトニアンモンテカルロ法
こちらのサイトで pdf ファイルをダウンロードできます。
先生、ありがとうございます！
今回のブログは以上です。
次回は 階層ベイズモデル を学びます。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
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の教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
いいなと思ったら応援しよう！
応援ありがとうございます。これからもがんばって記事を作成します！