「データ解析のための統計モデリング入門」をPythonで写経 Vol.18 ～ ９章「GLMのベイズモデル化と事後分布の推定」①ベイズ統計モデリング～ポアソン回帰モデル

９章「GLMのベイズモデル化と事後分布の推定」

書籍の著者　久保拓弥 先生

書籍「データ解析のための統計モデリング入門」９章「GLMのベイズモデル化と事後分布の推定」の Python写経活動記録 です。 

この記事は ポアソン回帰のGLM をベイズモデル化 します。
GLM モデル と ベイズモデル の２モデルを実践します。
そしてアディショナルでは「気になる２つのあれこれ」を気ままに書いてみました。

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

はじめに

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

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

準備

準備

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

■ 利用データ
テキスト・サポートサイトのデータファイルを引用しています。
▶️ サポートサイト

Jupyter Notebook ファイルと同一フォルダ内に「data」フォルダを用意して、data フォルダ配下の章別フォルダにデータファイルを格納しています。

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

# インポート

# 数値計算
import numpy as np

# 統計計算
import scipy.stats as stats
import statsmodels.api as sm
import statsmodels.formula.api as smf

# PyMC
import pymc as pm
import arviz as az

# Rデータセットの読み込み
import rdata

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

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

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

■ 統計モデル
ポアソン回帰の GLMモデルと、GLM のベイズ統計モデルです。

■ モデリング手続き

データの確認

データを読み込み、データの外観を眺めてから、統計モデルを選択するためのデータの特徴確認を行います。

◼️ RDataファイルの読み込み
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（説明変数）です。

🍀🍀🍀

◼️ データの確認
基本的な統計量やチャートでデータを概観します。

① 要約統計量の表示

# 要約統計量
data.describe().round(3)

【実行結果】
種子数 y の値は 3 ～ 12 個です。
体サイズの値については、最小値・四分位数・最大値が 3, 4, 5, 6, 7 とキリのよい連番になっていて、意味深です。

② 標本分散の表示

# 標本分散
data.var(ddof=1).rename('var').to_frame().T.round(3)

【実行結果】
ｙ の分散 5.8 は平均値 7.3 より少し小さな値です（過小分散？）

③ 相関係数の表示

# 相関係数
data.corr().round(3)

【実行結果】
葉重量 x と花重量 y の間には弱い正の相関があります。

④ ヒストグラムと散布図の描画
seaborn の pairplot() を利用します。
著者の先生には叱られるかもですが、回帰直線を添えます。

# ヒストグラムと回帰直線付き散布図のペアプロットの描画
sns.pairplot(data=data, kind='reg',
             plot_kws={'scatter_kws': {'s': 70, 'alpha': 0.6},
                       'line_kws': {'color': 'tab:red', 'alpha': 0.7}},
             diag_kws={'edgecolor': 'white', 'bins': 5});

【実行結果】
体サイズ x は一様分布に見えます。
種子数 y は峰が１つです。

【考察】

• 種子数 y のヒストグラムは峰が１つです。

• 体サイズ x のヒストグラムは一様分布を示しています。

• 横軸 x、縦軸 y の散布図（左下）は、体サイズ x が大きくなるにつれて、種子数 y が大きくなる傾向があります。

🍀🍀🍀

◼️ データの特徴まとめ

種子数のばらつきを「確率分布」で表現します。
種子数の特徴を表現できる確率分布は「ポアソン分布」です。

今回の統計モデル・ベイズ統計モデルは「個体ごとの平均種子数が体サイズにどう依存しているのかを調べる」（テキストp.194）ことを目的にします。

ポアソン回帰のGLM

統計モデルをデータに当てはめ（モデルの理解）

ポアソン回帰の GLM を観測データに当てはめます。

◼️ 確率分布、リンク関数、線形予測子
GLMの３要素である「確率分布」、「リンク関数」、「線形予測子」を定義します。

🔷 確率分布と確率質量関数
種子数 $${y_i}$$ はパラメータ $${\lambda_i}$$ のポアソン分布に従います。

🔷 リンク関数と線形予測子
リンク関数は「対数」、線形予測子は「$${\beta_1 + \beta_2 x_i}$$」です。
パラメータ $${\lambda_i}$$ は指数関数「$${\exp(線形予測子)}$$」で表現されます。

🍀🍀🍀

◼️ 真値 $${\beta_1=1.5, \beta_2=0.1}$$ で可視化する
テキスト p.195 図 9.1 で、パラメータの真値（例題データ生成時のパラメータ）が $${\beta_1=1.5, \beta_2=0.1}$$ であると示されています。

図 9.1 に相当する、真値を用いた平均種子数 $${\lambda}$$ と観測値の散布図チャートを描画します。

# (B)20個体の架空植物のサイズx_iと種子数y_1の関係 p.194~195 図9.1(B)

# x軸の値の設定
x_val = np.linspace(data.x.min(), data.x.max(), 101)

# 描画領域の設定
fig, ax = plt.subplots()
# x,yの散布図の描画
sns.scatterplot(data=data, x='x', y='y', s=80, alpha=0.7, label='観測値', ax=ax)
# λの真値の点線の描画
ax.plot(x_val, np.exp(1.5 + 0.1 * x_val), color='black', ls='--',
        label='平均$\\lambda$の真値:exp(1.5+0.1$x$)')
# 修飾
ax.set_xlabel('サイズ $x_i$', fontsize=12)
ax.set_ylabel('種子数 $y_i$', fontsize=12)
ax.set(xticks=[3, 4, 5, 6, 7])
ax.legend();

【実行結果】
平均種子数 $${\lambda}$$ は指数関数的に、緩やかな上昇カーブを描いています。

２つのモデルによるパラメータ推定値はこの真値に迫れるのでしょうか！？

統計モデルをデータに当てはめ（当てはめと評価）

◼️ 統計モデルをデータに当てはめ、の準備
GLMの３要素「確率分布」「リンク関数」「線形予測子」を用いる統計モデルをデータに当てはめします。
統計モデルの対数尤度 $${\log L}$$ が最大になるパラメータ $${\beta_1, \beta_2}$$ を推定します。

🔷 今回の統計モデルの対数尤度

🍀🍀🍀

◼️ Python ライブラリで統計モデルをデータに当てはめ
statsmodels の glm を利用して統計モデルを実装します。
引数 family で確率分布、link でリンク関数、formula で線形予測子を指定します。
当てはめ結果を変数 result に格納します。

# 切片と傾きの最尤推定値 p.194

# 設定: GLMの引数familyの値の設定（ポアソン分布）
family = sm.families.Poisson()

# ポアソン回帰の当てはめ
result_glm = smf.glm(formula='y ~ x', data=data, family=family).fit()

# 結果表示
result_glm.summary()

【実行結果】
最下２行の Intercept が $${\beta_1}$$、x が $${\beta_2}$$ に対応しており、テキスト p.50 下部の summary 関数の実行結果 Coefficients に相当します。

🍀🍀🍀

◼️ 当てはめ結果の分析

🔷 係数の推定値
Intercept と x の coef に注目します。
パラメータである係数の最尤推定値は 切片（Intercept）$${\beta_1 = 1.5661}$$、傾き（x）$${\beta_2 = 0.0833}$$ です。
真値に近い結果になりました。

この推定値を線形予測子に当てはめてみます。

🔷 標準誤差の推定値など
続いて標準誤差 std err です。
パラメータ $${\beta_1, \beta_2}$$ の推定値の標準偏差であり、それぞれ $${0.360, 0.068}$$ です。

この標準誤差を用いて計算する $${p}$$ 値は、$${\beta_1}$$ は 5% 水準で有意ですが、$${\beta_2}$$ は $${0.223 > 0.05}$$ となり 5% 水準で有意ではありません。

🍀🍀🍀

ところで、ポアソン分布を用いる GLM は「観測値の平均＝分散」を仮定していました。
ただし、例題データは平均よりも分散がやや小さくなっていました。

この分散が小ささが影響して標準誤差が不適切にならないか、気になって…
ChatGPT に標準誤差の評価を手伝ってもらうと…
壮大な物語（当社比）が返ってきました！
「アディショナルタイム」で詳細を記載いたします。

平均種子数の予測

図 9.1 の散布図＋パラメータ真値による予測曲線のチャートに、GLM による平均種子数 $${\lambda}$$ の予測曲線（95%区間付き）を重ねて描画します。

# 図9.1(B)にGLMによるλの予測値を追加して描画

## 設定と準備
# x軸の値の設定
x_val = np.linspace(data.x.min(), data.x.max(), 101)
# GLMのパラメータ推定値
b1_hat, b2_hat = result_glm.params
# GLMのλの予測値
pred = result_glm.get_prediction(dict(x=x_val))
lam_pred = pred.predicted_mean  # λの予測値
lam_confint = pred.conf_int()   # λの予測値の95%区間

# 描画領域の設定
fig, ax = plt.subplots()
# x,yの散布図の描画
sns.scatterplot(data=data, x='x', y='y', s=80, alpha=0.7, label='観測値', ax=ax)
# λの真値の点線の描画
ax.plot(x_val, np.exp(1.5 + 0.1 * x_val), color='black', ls='--',
        label='$\\lambda$の真値:exp(1.5+0.1$x$)')
# λの予測値の描画
ax.plot(x_val, lam_pred, color='tab:red',
        label=f'$\\lambda$の予測値:exp({b1_hat:.2f}+{b2_hat:.2f}$x$)')
# λの予測値の95%区間の描画
ax.fill_between(x_val, *lam_confint.T, color='lightpink', alpha=0.3,
        label=f'$\\lambda$の予測95%区間')

# 修飾
ax.set_xlabel('サイズ $x_i$', fontsize=12)
ax.set_ylabel('種子数 $y_i$', fontsize=12)
ax.set(xticks=[3, 4, 5, 6, 7])
ax.legend();

【実行結果】
GLM の予測曲線の傾きは真値よりも緩やかになりました。
95%区間は体サイズ 5 あたりが狭く、両端（3、7）は広くなっています。

🍀🍀🍀

おまたせしました！
初・ベイズ統計モデリングに進みましょう！

ベイズ統計モデル

ベイズ統計モデルの理解

◼️ ベイズ統計モデリングが始まる
ベイズ統計モデリングでは「パラメータの事後分布の推定」を中心に置いて動きます。
ベイズ統計モデルの事後分布は尤度と事前分布の積に比例します。

今回のベイズ統計モデリングで推定するパラメータは、線形予測子のパラメータ $${\beta_1, \beta_2}$$ です。
データを $${\bm Y}$$ とし、事後分布を数式化します。

事後分布の算出に必要な「尤度」と「事前分布」を確認します。

🍀🍀🍀

◼️ 尤度
個体 $${i}$$ の種子数 $${y_i}$$ はポアソン分布 $${p(y_i \mid \lambda_i)}$$ に従うとします。
平均パラメータ $${\lambda_i}$$ は線形予測子と対数リンク関数を用いて $${\lambda_i = \exp (\beta_1 + \beta_2 x_i)}$$ とします。
尤度は次の数式で表されます。

PyMC のコードに寄せて、次のように表現します。

🍀🍀🍀

◼️ 事前分布
事後分布の数式化では $${\beta_1}$$ の事前分布 $${p(\beta_1)}$$ と $${\beta_2}$$ の事前分布 $${p(\beta_2)}$$ が示されていました。
テキストは、今回のベイズ統計モデルの事前分布に「無情報事前分布」を用います。

今回のモデルの無情報事前分布が「線形予測子のパラメータの値は $${[-\infty, \infty]}$$ の範囲で好きな値をとってよい、といったことを表現する事前分布」である旨をテキストは説明しています。
また、以下の２種類の確率分布がよく使われていると説明されています。

今回のモデルでは、ひらべったい正規分布 $${\text{Normal}(0, 100^2)}$$ を用います。

ひらべったい正規分布を可視化で確認します。
テキスト p.197 図 9.2 に相当します。

# 切片β1と傾きβ2の無情報事前分布 p.197 図9.2

# x軸の値
x_val = np.linspace(-11, 11, 101)
# 描画領域の設定
fig, ax = plt.subplots()
# 標準正規分布の確率密度関数の点線の描画
ax.plot(x_val, stats.norm.pdf(x_val, loc=0, scale=1), color='gray', ls='--')
# 標準正規分布のテキストの表示
ax.text(x=-5.5, y=0.35, s='標準正規分布\n$N(0,1)$', fontsize=12)
# ひらべったい正規分布の確率密度関数の実線の描画
ax.plot(x_val, stats.norm.pdf(x_val, loc=0, scale=100), lw=3)
# ひらべったい正規分布のテキストの表示
ax.text(x=4, y=0.018, s='ひらべったい\n正規分布\n$N(0,100)$', color='tab:blue',
         fontsize=12, fontweight='bold')
# 修飾
ax.set_xlabel(r'$\beta_1$あるいは$\beta_2$の値', fontsize=12)
ax.set_ylabel(r'確率密度 $p(\beta_{*})$', fontsize=12)
ax.set(xlim=(-10.5, 10.5));

【実行結果】
正規分布 $${\text{Normal}(0, 100^2)}$$ はとても「ひらべったい」形状です。
これが「無情報っぽい」と呼ばれているんです！

今回のベイズ統計モデルの数式をまとめます。
ベイズ統計モデルを数式で表しておくと、ベイズ統計モデルのコード化の際に有用です（と思っています）。
なお、説明変数 x を中央化（中心化）するので、線形予測子では x から標本平均 $${\bar{x}}$$ を引いています。

$${\sim}$$ は左側の変数が右側の確率分布に「従う」ことを意味します。
$${=}$$ で表現された変数は、確率分布に従わず、等号で結ばれた数式どおりに特定の値が決定される「決定論的変数」です。

パラメータの事後分布の推定

PyMC ライブラリでベイズ統計モデリングを実装します。
テキストが用いるベイズ統計ツール「WinBUGS」の設定を解釈しつつ、PyMC コードに書き換えていきます。

早速、モデルの定義をします。

◼️ モデルの定義

# モデルの定義

# 設定と準備
mean_X = data.x.mean()          # 説明変数xの平均値、データの中央化で使用
coords = {'data': data.index}   # 座標ラベルの設定：データ行の識別子

# モデリング
with pm.Model(coords=coords) as model:
    
    # dataの定義
    # 目的変数 Y
    Y = pm.Data('Y', value=data['y'].values, dims='data')
    # 説明変数 X
    X = pm.Data('X', value=data['x'].values, dims='data')

    # 事前分布: 平均0, 標準偏差100の正規分布
    beta1 = pm.Normal('beta1', mu=0, sigma=100) 
    beta2 = pm.Normal('beta2', mu=0, sigma=100)

    # 線形予測子: λ=exp(線形予測子)。説明変数Xは中央化
    lam = pm.Deterministic(
        'lam', pm.math.exp(beta1 + beta2 * (X - mean_X)), dims='data')

    # 尤度: 平均λのポアソン分布
    obs = pm.Poisson('obs', mu=lam, observed=Y, dims='data')

【実行結果】なし

【コードの補足説明】
PyMC のベイズモデルを with 構文で定義します。
モデル名 model を設定します。

with pm.Model(coords=coords) as model:

dims（次元）を定義してモデルの可視化をわかりやすくします。
今回は、データの行の次元を data と命名し、data のラベルを coords で設定します。
data の coords の値は例題データのインデックスです（0～19）。

coords = {'data': data.index}

with 文内の定義を見ていきます。
pm.Data でデータを登録します。
Data の引数は、PyMC内部の変数名、value：値、dims：次元などです。
＝の左辺は Python の変数名であり、個人的には Python の変数名と PyMC 内部の変数名は合わせておきたいです。

Y = pm.Data('Y', value=data['y'].values, dims='data')
X = pm.Data('X', value=data['x'].values, dims='data')

パラメータの事前分布を設定します。
正規分布は確率分布である Normal クラスを利用します。
引数は PyMC内部の変数名、分布のパラメータ（平均 mu、標準偏差 sigma）、dims：次元などです。
先の数式どおり、$${\text{mu}=0, \text{sigma}=100}$$ を設定します。

beta1 = pm.Normal('beta1', mu=0, sigma=100) 
beta2 = pm.Normal('beta2', mu=0, sigma=100)

線形予測子を設定します。
Deterministic は決定論的変数などと呼ばれ、数式で $${\sim}$$（確率分布に従う）ではなく、$${=}$$ で示されるものであり、MCMC サンプルの対象にしたい場合に使います。
引数は PyMC内部の変数名、式、dims：次元などです。
指数関数 EXP は pm.math.exp() を利用します。

lam = pm.Deterministic(
        'lam', pm.math.exp(beta1 + beta2 * (X - mean_X)), dims='data')

尤度を設定します。
ポアソン分布は Poisson クラスで表現します。
引数は PyMC内部の変数名、分布のパラメータ（平均 mu）、observed：観測値、dims：次元などです。
観測値を設定する点が特徴的です。

obs = pm.Poisson('obs', mu=lam, observed=Y, dims='data')

PyMC の裏側では 多次元計算ライブラリ PyTensow が動きます。
ちなみに、PyMC モデルの定義時に配列の不整合などの問題が生じていると PyTensow に怒られます…

🍀🍀🍀

◼️ モデルの確認
モデルの数式と有向グラフを可視化します。
数式表示はモデル名 model で呼び出せます。

# モデルの表示
model

【実行結果】
数式でざっくりベイズ統計モデルを確認できます。

モデルの有向グラフは model_to_graphviz(モデル名) で描画できます。
graphviz ライブラリを事前にインストールする必要があります。

# モデルの可視化
pm.model_to_graphviz(model)

【実行結果】
変数のつながりやデータ数が分かりやすいですね（個人的な感想です）。
平均種子数を示す lam は個体数の 20 個に対応しています。
パラメータ beta1, beta2 はそれぞれ１個です。

ChatGPTによる「図の凡例」です。
不正確な部分がありましたらごめんなさい。

🍀🍀🍀

◼️ MCMC の実行
MCMC サンプリングを行います。
MCMC アルゴリズムは NUTS（No-U-Turn Sampler）です。

こちらも with 文を用います。
idata（推論データ）に pm.sample() で生成した MCMC サンプルなどのデータが代入されます。
pm.sample() の引数は次の通り。

%%time
# MCMCサンプリング
# thinning:未実施, パラメータの初期値:未設定, サンプリングアルゴリズム:NUTS

with model:
    idata = pm.sample(draws=1500, tune=100, chains=3, random_seed=42)

【実行結果】
処理時間は約 45 秒でした。

パラメータ推定値の確認

ここからは ChatGPT が教えてくれた「健全性チェック」で MCMC で得られた事後分布からのサンプルを確認します。

🔷 1. ダイバージェンス（Divergences）

• 意味：サンプラー（特に HMC/NUTS）が不安定になったかどうか

• 理想値：0（ダイバージェンスがない）

• 多い場合 → モデルの再パラメータ化や事前分布の調整が必要

▶ MCMC実行の図の「Divergences」はすべて０です！
個別に確認するコードは以下の通りです。

# ダイバージェンスの個数を数える
print('Divergenceの個数:', idata.sample_stats.diverging.data.sum())

【実行結果】

🔷 2. R-hat（Gelman-Rubin 診断）

• 意味：チェーン間の収束度合い（between-chain vs within-chain variance）

• 理想値：1.00（1.01 未満が望ましい、1.05 を超えると要注意）

• 値が大きい → 収束していない可能性

R-hat が 1.01 以下になっていることを確認します。
arviz の rhat() を利用します。

# r_hat>1.01の確認
# 設定
idata_in = idata         # idata名
threshold = 1.01         # しきい値

# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())

【実行結果】
パラメータ各行の右端の数字が「R-hat ＞1.01」の個数です。
３つのパラメータの「R-hat ＞1.01」の個数が０です。
▶ すべての R-hat が 1.01 以下です。

🔷 3. 有効サンプルサイズ（Effective Sample Size, ESS）

• 意味：自己相関を考慮した「実質的に独立なサンプル数」

• 2種類：

  • Bulk ESS：平均や分散の推定精度

  • Tail ESS：分布の端（信頼区間の端）の推定精度

• 理想値：数百以上（特に 400 以上あれば安心）

ESS が 1000 以上になっていることを確認します。
arviz の ess() を利用します。

# 有効サンプルサイズ しきい値を1000とした
print((az.ess(idata) < 1000).sum())

【実行結果】
パラメータ各行の右端の数字が「ESS < 1000」の個数です。
３つのパラメータの「ESS < 1000」の個数が０です。
▶ すべての ESS が 1000 以上です。

🔷 4. サンプルのトレースプロット（Traceplot）

• 意味：チェーンが混ざり合い、安定して分布を探索しているかを可視化

• 理想形：各チェーンが「毛糸玉」のように混ざり合い、ドリフトがないこと
  (ドリフトがないとは、MCMC の各チェーンが時間経過に従って少しずつずれていくような傾向（トレンド）を持たないこと)

• 収束が悪い場合 → 鎖ごとに異なる平均に固まる

トレースプロットを描画します。
PyMC（または arviz ）の plot_trace() を利用します。
テキスト p.209 図 9.5 に相当します。

# トレースプロットの表示 p.209 図9.5に相当
var_names = ['beta1', 'beta2']
pm.plot_trace(idata, var_names=var_names, figsize=(7, 5))
plt.tight_layout();

【実行結果】
右側がトレースプロットです。
▶ 各チェーンが「毛糸玉」のように混ざり合い、ドリフトがないです。
また、左側はパラメータ推定値の事後分布の確率密度（カーネル密度推定）です。
３本のチェーンはほぼ同じ形状をしています、ＯＫです。

🔷 5. エネルギー診断（Energy Bayesian Fraction of Missing Information, E-BFMI）

• 意味：HMC サンプラーの「飛び方」が十分多様かをチェック

• 理想値：0.3 以上（小さいとサンプラーの動きが制限されている）

• PyMC / ArviZ の `plot_energy()` で可視化できる

エネルギープロットを描画します。
PyMC（または arviz ）の plot_energy() を利用します。
凡例に BFMI 値が表示されます。

# エネルギープロット ※BFMI値を表示
pm.plot_energy(idata, fill_color=('tab:blue', 'tab:red'), fill_alpha=(0.7, 0.3));

【実行結果】
▶ ３本のチェーンのBFMIは 0.3 を超えています。

🔷 6. 事後分布の確認

• サマリー統計量（平均・標準偏差・HDI）を確認

• 事前分布や真値と比べて「あり得る範囲か」を検証

事後分布の要約統計量（サマリー統計量）を表示します。
テキスト p.211 の「事後分布の統計量」に相当します。
PyMC（または arviz ）の summary() を利用します。

# 推論データの要約統計情報の表示 p.211の図に相当
summary = pm.summary(idata, hdi_prob=0.95, var_names=var_names, round_to=3)
summary

【実行結果】
実は事後分布の要約統計量の右端に R-hat と ESS が表示されています。

【 事後分布の要約統計量の項目説明 】
表の各項目の意味と確認ポイントを ChatGPT にまとめてもらいました。
不正確な部分がありましたらごめんなさい。

🍀🍀🍀

◼️ 95% 区間を深堀りする
テキストは $${\beta_2}$$ の 95% 信用区間が０を含んでいることに言及しています。
この表の $${\beta_2}$$ の 95% 区間も０を含んでいます。
$${\beta_2}$$ が０（つまり傾きが無い）から充分に離れているとは言いにくい感じです。

① 可視化
パラメータ推定値の 95% 区間を可視化しましょう。
フォレストプロットを描画します。
PyMC（または arviz ）の plot_forest() を利用します。

# フォレストプロット
pm.plot_forest(idata, combined=True, var_names=['beta1', 'beta2'], hdi_prob=0.95)
plt.axvline(0, color='tab:red', ls='--')
plt.grid(axis='x', alpha=0.5);

【実行結果】
$${\beta_2}$$ の 95% 区間が０を含むことを直感的に把握できました。

② 推定値 $${\hat{\beta}_2}$$ の解釈
$${\exp(\hat{\beta}_2)}$$ を見ましょう。

# exp(beta2)の算出
(
    np.exp(summary.loc['beta2', ['hdi_2.5%', 'mean', 'hdi_97.5%']])
    .rename('exp(beta2)')
    .to_frame()
    .round(2)
    .T
)

【実行結果】
HDI の 2.5%点では 0.96、つまり、体サイズ $${x}$$ が１単位増えると平均種子数 $${\lambda}$$ は 0.96 倍になり、減少します。
HDI の 97.5%点では 1.25、つまり、体サイズ $${x}$$ が１単位増えると平均種子数 $${\lambda}$$ は 1.25 倍になり、増加します。

$${\beta_2 > 0}$$ となる確率を計算してみます。
テキスト p.212 ～ 213 の計算に相当します。

# β2>0となる確率 p.212~213
(idata.posterior.beta2.data > 0).mean()

【実行結果】
$${\beta_2 > 0}$$ となる確率は 0.89 です。

③ HDI（最高密度区間） とは
ところで、上の表で 95% 区間は「HDI 2.5%」「HDI 97.5%」の列名で表示されています。
hdi は「HDI（最高密度区間）」を差します。
95% HDI は 95% 最高密度区間であり、95% 信用区間とは異なる概念です。

PyMC・ArviZ 界隈は HDI を用いるようですので、この機会に HDI の理解を進めましょう。

【 95% HDI の直感的理解 】
下の図のイメージで考えてみます。

（注）
図の分布は、今回のベイズ統計モデルのパラメータ事後分布とは異なりますので、ご留意下さい。

95% 信用区間は、区間の外側の確率が等しく、信用区間の外に区間内の確率密度よりも高い点があります。

95% HDI は確率密度が最も高い区間なので、HDI 区間の外に区間内の確率密度よりも高い点がありません。
図では、 95% HDI の青いエリアの確率密度は赤い点線以上です。
95% HDI の外側の確率密度は赤い点線よりも小さいです。

さらに 95% HDI の区間は、95% に達する 最短の区間 だそうです。

テキストは信用区間を用いています。
この記事シリーズは HDI を用います。
区間の定義が相違しますので、ご留意下さい。

平均種子数の予測

◼️ テキストのチャート
テキスト p.210 図 9.6 (A) に相当する「平均種子数 $${\lambda}$$ の予測チャート」を描画します。
あわせて 同図 (B) に相当する「$${\beta_1, \beta_2}$$ のMCMCサンプルの散布図」を描画します。

# （A) 平均λの予測、（B)β1,β2の事後分布の描画 p.210 図9.6
# ★(B)のβ1,β2の相関は、中央化：テキスト同様の相関なさげ、非中央化：強い相関、となる

## 描画用データの取得
# x軸の値
x_vals = np.linspace(data.x.min(), data.x.max(), 101)
# idataをthinning（2個飛ばし）
thinned_idata = idata.sel(draw=slice(None, None, 3))
# β1, β2の事後分布からのMCMCサンプル(thinning後)を平坦化して取得
beta1s = thinned_idata.posterior.beta1.stack(sample=('chain', 'draw')).data
beta2s = thinned_idata.posterior.beta2.stack(sample=('chain', 'draw')).data
# β1, β2のサンプル(thinning前)の中央値を取得
b1_median = idata.posterior.beta1.median().data
b2_median = idata.posterior.beta2.median().data

## 描画設定
# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4), tight_layout=True)

## (A)の描画
# 観測値の散布図の描画
sns.scatterplot(data=data, x='x', y='y', s=70, ec='tab:blue', fc='white',
                label='観測値', ax=ax1)
# 500×3組のβ1,β2のサンプル列を使って平均λ=exp(β1+β2x)の予測値を赤い透過線で描画
for b1, b2 in zip(beta1s, beta2s):
    ax1.plot(x_vals, np.exp(b1 + b2 * (x_vals - mean_X)), lw=0.1, color='tab:red',
    alpha=0.1)
# 事後分布の中央値を使った予測値を赤い太線で描画
ax1.plot(x_vals, np.exp(b1_median + b2_median * (x_vals - mean_X)), lw=3,
         color='tab:red', label='予測値:中央値')
# 修飾
ax1.set_xlabel('植物の体サイズ $x_i$', fontsize=12)
ax1.set_ylabel('種子数 $y_i$', fontsize=12)
ax1.set_title(r'(A) $x$ に依存する平均 $\lambda$ の予測')
ax1.legend()

## (B)の描画
# β1, β2のサンプルの散布図を描画
sns.scatterplot(x=beta1s, y=beta2s, s=50, alpha=0.3, ax=ax2)
# 修飾
ax2.set_xlabel(r'切片 $\beta_1$', fontsize=12)
ax2.set_ylabel(r'傾き $\beta_2$', fontsize=12)
ax2.set_title(r'(B) 事後分布 $p(\beta_1, \beta_2 \mid \boldsymbol{Y})$');

【実行結果】

【チャートの作り方と解釈】

① 右の $${\lambda}$$ の予測
$${\beta_1, \beta_2}$$ のMCMCサンプル 1500 個から 500 個を抜き取って、各サンプルの予測値を「薄く」描画しています。
太い赤実線は、$${\beta_1, \beta_2}$$ のMCMCサンプルの中央値を使った予測値です。
$${\lambda}$$ の予測値は $${\exp(\beta_{1_{\text{sample}}} + \beta_{2_{\text{sample}}} (x - x_{\text{mean}}))}$$ で計算しました。

体サイズが大きくなるにつれて、平均種子数 $${\lambda}$$ が（ちょっとずつ）大きくなる様子が分かります。
MCMCサンプル単位の薄い線は、中央値の線の周りに縦方向の広がりをもって、引かれています。
信用区間やHDIを思わせますね。
あとで、予測値の 95% 信用区間・95% HDI を描画してみましょう。

② 左の $${\beta_1, \beta_2}$$ のMCMCサンプルの散布図
同時分布である $${\beta_1, \beta_2}$$ の事後分布の描画です。
テキストにならって「$${\beta_1, \beta_2}$$ の MCMC サンプル感には強い相関がないようだ」と思っておきます。

🍀🍀🍀

◼️ $${\lambda}$$ の予測値の 95% 区間
GLM の予測のときに描画した ↓ このチャート ↓ をベイズも同様に描画します。 

まず 95% 信用区間を描画します。
MCMC サンプルの 2.5% 点～ 97.5% 点を numpy の quantile() （分位数）で算出して、間を塗りつぶします。

# 図9.1(B)にベイズ推定によるλの予測値を追加して描画 95%CI版

## 設定と準備
# x軸の値の設定
x_val = np.linspace(data.x.min(), data.x.max(), 101)
# パラメータのMCMCサンプルを取得
b1s = idata.posterior.beta1.data.flatten()
b2s = idata.posterior.beta2.data.flatten()
# 1500個のλ予測値を算出
lam_preds = np.exp(b1s + np.outer((x_vals - mean_X), b2s)).T
# λ予測値の中央値の算出
lam_median = np.median(lam_preds, axis=0)
# λ予測値の95%区間の算出
lam_95ci = np.quantile(lam_preds, q=[0.025, 0.975], axis=0)

# 描画領域の設定
fig, ax = plt.subplots()
# x,yの散布図の描画
sns.scatterplot(data=data, x='x', y='y', s=80, alpha=0.7, label='観測値', ax=ax)
# λの真値の点線の描画
ax.plot(x_val, np.exp(1.5 + 0.1 * x_val), color='black', ls='--',
        label='$\\lambda$の真値:exp(1.5+0.1$x$)')
# 事後中央値を使った予測値を描画
ax.plot(x_vals, lam_median, color='tab:red',
        label='$\\lambda$の予測値:'
              f'exp({np.median(b1s):.2f}+{np.median(b2s):.2f}$x$)')
# λの予測値の95%区間の描画
ax.fill_between(x_val, * lam_95ci, color='lightpink', alpha=0.3,
                label=f'$\\lambda$の予測95%区間')

# 修飾
ax.set_xlabel('サイズ $x_i$', fontsize=12)
ax.set_ylabel('種子数 $y_i$', fontsize=12)
ax.set(xticks=[3, 4, 5, 6, 7])
ax.legend();

【実行結果】
GLM とほぼほぼ同じです。

続いて 95% HDI を塗ります！
PyMC（または arviz ）の plot_hdi() を利用します。

# 図9.1(B)にベイズ推定によるλの予測値を追加して描画 95%HDI版

## 設定と準備
# x軸の値の設定
x_val = np.linspace(data.x.min(), data.x.max(), 101)
# パラメータのMCMCサンプルを取得
b1s = idata.posterior.beta1.data.flatten()
b2s = idata.posterior.beta2.data.flatten()
# 1500個のλ予測値を算出
lam_preds = np.exp(b1s + np.outer((x_vals - mean_X), b2s)).T
# λ予測値の中央値の算出
lam_median = np.median(lam_preds, axis=0)

# 描画領域の設定
fig, ax = plt.subplots()
# x,yの散布図の描画
sns.scatterplot(data=data, x='x', y='y', s=80, alpha=0.7, label='観測値', ax=ax)
# λの真値の点線の描画
ax.plot(x_val, np.exp(1.5 + 0.1 * x_val), color='black', ls='--',
        label='$\\lambda$の真値:exp(1.5+0.1$x$)')
# 事後中央値を使った予測値を描画
ax.plot(x_vals, lam_median, color='tab:red',
        label='$\\lambda$の予測値:'
              f'exp({np.median(b1s):.2f}+{np.median(b2s):.2f}$x$)')
# λの予測値の95%区間の描画 ★plot_hdi 利用
pm.plot_hdi(x_val, lam_preds.reshape(3, 1500, 101), ax=ax,
            hdi_prob=0.95, color='lightpink',
            fill_kwargs={'alpha': 0.3},
            backend_kwargs={'label': f'$\\lambda$の予測95%区間'})

# 修飾
ax.set_xlabel('サイズ $x_i$', fontsize=12)
ax.set_ylabel('種子数 $y_i$', fontsize=12)
ax.set(xticks=[3, 4, 5, 6, 7])
ax.legend();

【実行結果】
95% 信用区間よりも（ほんの少し）幅が広い印象です。

ベイズの基本事項を目一杯詰め込みました。
お腹いっぱいですね (^ ^;)
休憩しましょう。

（注意）
アディショナルの２テーマは趣味的な深堀りとコードです。
ご興味ない方はスルーしてくださって大丈夫です。

アディショナル：GLM標準誤差の評価

GLMの標準誤差の評価を深める

ポアソン回帰の GLM のパラメータ推定値の標準誤差を考察します。
GLM のサマリー表を再掲します。

標準誤差は後に続く $${z}$$ 値、$${p}$$ 値、95%信頼区間の値に影響します。
分散がやや過小気味な点が気になって、ChatGPT に標準誤差の評価をしてもらいました。

🍀🍀🍀

🔷 まず、標準誤差（SE）って何？

• 係数（切片や x の係数）は「データをもう一度取り直したら、どのくらい揺れそうか」を表す“幅”がつきます。これが標準誤差（SE）です。

• SE が小さいほど係数は “はっきり” 推定できています。
  大きいほど係数の推定値は “ぼんやり” します。

🍀

🔷 GLM の結果の表は何を前提にしてる？

• モデルは ポアソン回帰（リンク＝対数）です。

• 表には「nonrobust」「Scale=1.0」と表示されています。
  この表示は「ポアソン分布の理想（分散＝平均）」が完全に成り立つと仮定してSEを計算していることを示しています。

例題データとこのモデルの仮定が合っていれば、nonrobust の SE は “理想的” です。
ずれていると、SE が大きすぎたり（過分散）／小さすぎたり（過小分散）します。

🍀

🔷 仮定が合っているかを、2つの数字でざっくり点検

表にある

• Pearson chi2 = 13.6

• Deviance = 14.171

• 残差自由度 df_resid = 18

から“分散のずれ具合”を測れます。

(1) Pearson χ² / df

(2) Deviance / df

• どちらも 1 に近ければOK。

• 1 より大きい → 過分散（SEは本来もっと大きいはず）。

• 1 より小さい → 過小分散（SEは本来もう少し小さいはず）。

今回は 0.76～0.79 と 1より少し小さいので、軽い過小分散（データが“やや静か”）傾向を示しています。
ただし n=20 と小標本ですし、この程度のズレは偶然でも起こりうる範囲ですので、有意におかしいとは言い切れません。

🍀

🔷 「スケール調整」してみるとどうなる？

スケール調整 SE を計算します。
“本当の分散係数” を $${\hat\phi \approx 0.756}$$（$${\hat\phi_\text{Pearson}}$$）とみなすと、SE は $${\sqrt{\hat\phi}}$$ 倍になります。

つまり、SE が約 13% 小さくなるだけです。

• 切片の SE：$${0.360 \times 0.869 \approx 0.313}$$ 

• x の SE：$${0.068 \times 0.869 \approx 0.059}$$

このときのスケール調整後 95%CI（ざっくり）は

• $${x : 0.0833 \pm 1.96 \times 0.059  \Rightarrow [-0.032,\ 0.199]}$$
  （元の CI は $${[-0.051,\ 0.217]}$$）
  → 結論「x のパラメータ推定値は有意でない」は変わりません。

🍀

🔷 念のための“安全確認”（初心者向けベストプラクティス）

1. 頑健（ロバスト）SEも見ておく

   • コード例：$${\mathtt{fit(cov\_type='HC3')}}$$（小標本でおすすめ）

   • モデルの分散仮定が少しズレていても壊れにくい“保険”です。

   • 今回のように $${\hat\phi<1}$$ なら、nonrobust と大きくは変わらないはず。

2. スケール調整 SE（準ポアソン的）を併記

   • $${\hat\phi=\text{Pearson χ²}/\text{df}}$$ を掛けて計算した SE と、nonrobust の SE を並べて比較。

   • 2つが近ければ「SEは安定している」と言いやすいです。

3. 外れ値・影響点の確認

   • デビアンス残差やレバレッジのプロットで、極端な点が SE を膨らませていないかを簡単にチェック。

🍀

🔷 最終まとめ（超やさしく）

• 表の SE は「ポアソンの理想」が成り立つ前提の値。

• データは 少しだけ過小分散っぽい
  → スケール調整 SE で 13% ほど小さくしても結論は同じ。

• 念のため 頑健SE を見ても、x の“有意でない”は変わらないはず。

• よって、この標準誤差は実用的に妥当 と判断して大きな問題はありません。
  （ただし n=20 と小さいので、ロバストSEの併記が親切です）

🍀🍀🍀

【小まとめ】
データは少々過小分散のようですが、スケール調整 SE を用いても x のパラメータ推定値は「有意ではありません」との判断に変わりはありません。
通常の（nonrobust）GLM の標準誤差は実用的に妥当とのことです。

ロバスト標準誤差

◼️ 頑健 SE ～ ロバスト標準誤差、現る
新たに「頑健SEを見ても～有意でないは変わらないはず」という言葉が出現しました…気になります…

ひとまず、頑健 SE ＝ロバスト標準誤差 のことをChatGPTにざっくり訊いてみました。

ロバスト標準誤差は、標準誤差にデータのばらつきを加味したもののようです。

🍀🍀🍀

◼️ ロバスト標準誤差を求める
statsmodels でサクッとロバスト標準誤差のGLMをモデリングします。
「HC3」（不均一分散頑健共分散行列推定v3）という方法で標準誤差を調整します。
$${\mathtt{fit(cov\_type='HC3')}}$$ と設定します。

# ロバスト標準誤差のGLM

# 設定: GLMの引数familyの値の設定（ポアソン分布）
family = sm.families.Poisson()

# ポアソン回帰の当てはめ
result_glm_hc3 = smf.glm(formula='y ~ x', data=data, family=family
                         ).fit(cov_type='HC3')

# 結果表示
result_glm_hc3.summary()

【実行結果】
パラメータの推定値は変わりませんが、標準誤差、$${z}$$ 値、$${p}$$ 値、95% 信頼区間の値が変化しました。

◼️ 結果の分析
最初の GLM の結果と、ロバスト標準誤差 HC3 のGLM の結果を比べてみます。
説明は ChatGPT に委ねます！

タイトル「GLMの結果とHC3の結果をどう評価する？（やさしく丁寧に）」

🔷 結論（先に一言）

• どちらも5%水準では有意ではありません。
  点推定では x を1増やすと発生率は +8.7%（IRR=1.087）ですが、95%CIは1をまたぐため「上がるとも下がるとも断言できない」結果です。

• HC3は“SEの作り方”だけを頑丈にしたもの。平均の形（リンクや説明変数）は同じまま、実データのばらつき方をうまく反映した 標準誤差（SE）・p値・CI を出します。

• 今回は軽い過小分散（Pearson χ²/df ≈ 0.76）があるため、通常SEより HC3のSEが小さくなり、p値は 0.22 → 0.084 に近づきました。
  それでも 5%では非有意という意思決定は 変わりません。

🍀

🔷 ステップ1：通常のGLM（nonrobust）の読み方

• 係数（リンク上）: $${\hat\beta = 0.0833}$$
  → IRR（現実の倍率）に戻すと $${\exp(0.0833)}$$= 1.087（= +8.7%）

• 標準誤差: SE = 0.068
  → Wald 統計量 $${z=\hat\beta/SE=1.23}$$ → p≈0.22（両側）
  → 95%CI（係数）: $${0.0833 \pm 1.96\times0.068 = [-0.051,\ 0.217]}$$
  　 95%CI（IRR）: $${\exp([-0.051,\ 0.217]) = [0.95,\ 1.24]}$$
  ⇒ IRRの区間が 1 をまたぐので、「効果あり」とは言い切れない。

• MDE（最小検出可能効果）（5%、検出力80%の親指ルール）
  $${\text{MDE} \approx 2.8\times SE = 2.8\times0.068 = 0.190}$$
  　倍率で: $${\exp(\pm0.190) \approx +21\%\ / -17\%}$$
  ⇒ この設計だと ±約20% 超の効果なら拾えるが、±10%前後は見分けにくい。

🍀

🔷 ステップ2：HC3（ロバスト標準誤差）とは？

• 何をする？
  係数の点推定は そのままに、SEの計算式だけを“実データの揺れ方”に合わせて頑丈に作り直します（サンドイッチ分散のHC3）。

• 何に強い？
  異分散、分散関数のズレ（過分散／過小分散）、一部の外れ値・高レバレッジによるSEの歪み、など。
  ※ 平均構造の誤指定（リンク不適合・欠落変数）による係数のバイアスそのものは直せません。

• 今回どうなる？
  軽い過小分散のため、HC3のSEは小さく出て、p値は境界に近づきます。

🍀

🔷 ステップ3：HC3で計算し直す（同じ係数、SEだけ変更）

• $${\hat\beta = 0.0833}$$（同じ）

• HC3-SE = 0.048
  → $${z = 0.0833/0.048 = 1.73}$$ → p≈0.084（両側）
  → 95%CI（係数）: $${0.0833 \pm 1.96\times0.048 = [-0.011,\ 0.177]}$$
  　 95%CI（IRR）: $${\exp([-0.011,\ 0.177]) = [0.989,\ 1.194]}$$
  ⇒ 区間は まだ1をまたいでいる（5%では非有意）。

• MDE（HC3）
  $${\text{MDE} \approx 2.8\times0.048 = 0.134}$$
  　倍率で $${\exp(\pm0.134) \approx +14\%\ / -13\%}$$
  ⇒ “検出できる最小効果”は ±約14% に改善（ただし結論は非有意のまま）。

🍀

🔷 ステップ4：2つを並べてどう判断する？

• 通常SE：p=0.22、95%CI（IRR）=[0.95, 1.24]、MDE≈±20%

• HC3：p=0.084、95%CI（IRR）=[0.989, 1.194]、MDE≈±14%

• 意思決定：5%基準ならどちらも非有意で一致。
  HC3で“境界”に近づいたが、区間はまだ1をまたぐ → 結論は変えないのが安全。

• なぜ一致を確認するの？
  通常SEとロバストSE（HC3）で結論が同じなら、SEの妥当性への不安が小さくなり、報告が安心して出せます。
  逆に結論が割れるなら、クラスタ相関（同一個体内・施設内）や自己相関（時系列）に強い手法（cluster-robust / HAC）や、分布の見直し（負の二項・ZIP/ZINB）を検討します。

🍀

🔷 実務に載せる「一文」テンプレ（今回の数値に合わせた完成版）

🍀

🔷 いっしょに計算してみよう（手順つき）

🔢 手計算のポイント

1. IRR = $${\exp(\beta)}$$

2. 95%CI（係数） = $${\beta \pm 1.96\times SE}$$
   → 95%CI（IRR）は範囲の両端に $${\exp(\cdot)}$$ をかける

3. Waldのp値（両側）:$${z=\beta/SE}$$ を標準正規分布で

4. MDE（5%・80%検出力の近似）: $${2.8\times SE}$$（係数スケール）
   → IRRの最小検出倍率 = $${\exp(\pm\,\text{MDE})}$$

🔢 具体例（今回の数値）

• 通常SE: $${\beta=0.0833,\ SE=0.068}$$
  IRR = $${\exp(0.0833)=1.087}$$
  95%CI（係数） $${0.0833\pm1.96\times0.068=[-0.051, 0.217]}$$
  　→ 95%CI（IRR） = [0.95,1.24][0.95, 1.24]
  MDE = $${2.8\times0.068=0.190}$$
  　→ exp⁡(±0.190)≈+21%, −17%\exp(\pm0.190) \approx +21\%,\ -17\%

• HC3: $${\beta=0.0833,\ SE=0.048}$$
  IRR = 同じ 1.087（係数自体は不変）
  95%CI（係数） = $${0.0833\pm1.96\times0.048=[-0.011, 0.177]}$$
  　→ 95%CI（IRR） = [0.989,1.194][0.989, 1.194]
  MDE = $${2.8\times0.048=0.134}$$
  　→ $${\exp(\pm0.134) \approx +14\%,\ -13\%}$$

🍀

🔷 よくある質問（超やさしく）

Q. ロバストSEにすると“有意にする”ための秘策みたいでズルでは？
A. いいえ。係数は変えず、SEの作りだけを現実のばらつきに合わせる調整です。過小評価や過大評価を避けるのが目的で、ズルではありません。

Q. HC3と通常SEで結論が違ったら？
A. データの構造（同一個体内・施設内のクラスタ相関、時系列の自己相関）を疑い、$${\mathtt{cov\_type='cluster'}}$$（groups=…）や $${\mathtt{cov\_type='hac-groupsum'}}$$ を検討します。
　必要なら分布自体を Poisson→負の二項、ゼロ過多なら ZIP/ZINB に見直します。

Q. ロバストSEでも“平均の形が間違っていたら”ダメなの？
A. はい。リンクの不適合や欠落変数による係数バイアスは直せません。
　残差図やスプライン、相互作用で平均構造の適合も必ず点検しましょう。

🍀

🔷 まとめ

• 目的：現実のばらつき方を反映した妥当なSEでCI・pを出すこと。

• 今回：HC3にしても5%では非有意という意思決定は同じ。

• 伝え方：IRR・95%CI・p・MDE・診断一言（Pearson χ²/df）を一文セットで出す。

• 次の一手：±10%級の効果を確かめたいなら、サンプル増や設計改善を検討。

【小まとめ】
ロバスト標準誤差（頑健SE）も x のパラメータ推定値は「有意ではありません」との判断に変わりはありません。

【まとめ】
データにやや過小気味の分散があるものの、通常の（nonrobust）GLM の標準誤差は実用的に妥当と考えられます。

アディショナル：PyMC「idata」の操作

謎な idata の謎

PyMC の sample() 関数の結果を「idata」という名で保存して、これが MCMC サンプルだ、として扱ってきました。
idata は InferenceData（推論データ）を短くもじったものです。

デフォルトでは ArviZ ライブラリの InferenceData 形式になっています。
ArviZ ライブラリはベイズモデルの可視化を分担する嬉しいライブラリであり、InferenceData 形式のベイズモデルデータを扱えます。

idata の中に posterior（事後分布）などのグループが分けられているのは、ArViZ 仕様なのです。

さらに InferenceData は Xarray ライブラリの xarray.Dataset をベースに作られています。

個人的な感覚では、ベイズ的な操作は ArviZ のお作法で、細かなデータの操作は Xarray のお作法を、それぞれ使い分けする感じです。

いろいろと覚えることが増えて、嫌になりますね…
ですが、わざわざ ArviZ や Xarray のお作法を必ず使わないといけないか、というと、そうでも無いです。
使い慣れた numpy 配列に変換して、データ操作は numpy でがんばるという選択肢があります。

numpy 化する手順を検討しましょう。

idata を少々いじってみる

Jupyter Notebook や VSCode でJupyter Notebook 形式で動かしてみます。

① idata を開いてみる
idata を表示します。

# 推論データの表示１
idata

【実行結果】
InferenceData 形式で４つのグループがあることが分かりました。

「▶ posterior」（事後分布）をクリックすると…

何か開きましたね！
posterior の下には「xarray.Dataset」と表示されています。
各グループのデータは Xarray の Dataset 形式になっています。

Data variables には beta1、beta2、lam の事後分布からのサンプル（MCMC サンプル）が格納されています。
Coordinates には 次元の名称である Chain（MCMC のchain）、Draw（MCMC のdraw）、data（PyMCモデルの次元 data）が見えています。

ちょっと親近感が湧いてきました！

posterior 指名で見てみましょう。

# 推論データの表示２：事後分布
idata.posterior

【実行結果】
posterior グループだけを表示しています。
上の図と同じ内容ですね！

beta1 を見てみましょう。

# 推論データの表示３：beta_1の事後分布サンプル
idata.posterior.beta1

【実行結果】

データ形式が Xarray の DataArray に変わりました。
ドラム缶の隣にたくさんの数字が並んでいます。
これが MCMC サンプルです！
行が chain、列が draw の２次元です。
shape は (3, 1500) です。

🍀🍀🍀

② numpy 配列化
beta1 を numpy 配列に変換しましょう。
おまじないは「to_numpy()」です。

# 推論データの表示４：beta_1の事後分布サンプルをnumpy配列化
beta1_samples = idata.posterior.beta1.to_numpy()
beta1_samples

【実行結果】
無事に numpy 化できました！

配列の形状を確認します。

# numpy 配列のshape
beta1_samples.shape

【実行結果】
chain=3, draw=1500 の２次元のままです。

１次元にしましょう。

# numpy 配列を平坦化
beta1_samples_1d = beta1_samples.flatten()  # または .ravel()

print('beta1_samples_1d.shape:', beta1_samples_1d.shape)
beta1_samples_1d

【実行結果】
扱いやすくなりました。

この numpy 配列を用いて、MCMC サンプルから事後分布の各種統計量を計算してみましょう。

# MCMC サンプルの numpy 配列を操作して事後平均等を計算
print(f'事後平均　: {beta1_samples_1d.mean()}')
print(f'標準偏差　: {beta1_samples_1d.std(ddof=1)}')
print(f'中央値　　: {np.median(beta1_samples_1d)}')
print(f'95%CI low: {np.quantile(beta1_samples_1d, 0.025)}')
print(f'95%CI upp: {np.quantile(beta1_samples_1d, 0.975)}')

【実行結果】
じゃーん✨️ 簡単ですね！

🍀🍀🍀

③ Xarray 風で平坦化
InferenceData 形式・Xarray 形式にしておきたいけど、次元が多いと操作が大変なので、平坦化したい…
お悩みあるあるですね！
stack メソッドを使って解決しましょう。
posterior（事後分布）グループ配下のデータを平坦化します。

# stackでchainとdrawを平坦化する
idata_stack = idata.posterior.stack(sample=('chain', 'draw'))
idata_stack

【実行結果】
chain & draw の２次元を sample という１つの次元に変換しました。

beta1 の MCMC サンプルを見てみます。

# 平坦化データのbeta_1にアクセス
idata_stack.beta1

【実行結果】
無事に１次元に変換されています。

もちろん、numpy 配列化もできます。

# 平坦化データのbeta_1をnumpy配列化
idata_stack.beta1.to_numpy()

【実行結果】
１次元で取り出せました。

🍀🍀🍀

④ ArviZ のツールで「抽出」する
ArviZ にも便利ツールが用意されています。
「extract」を使って平坦化しましょう。

# 推論データをextractで抜き取る：posterior階層が無くなり、chainとdrawが統合される
idata_extract = az.extract(idata)
idata_extract 

【実行結果】
posterior グループを抽出して、chain & draw の２次元が sample という名の１次元に変換されました！

extract したデータの beta1 を見ます。

# 抜き取りデータのbeta_1にアクセス
idata_extract.beta1

【実行結果】
１次元になっています！（sample: 4500）

もちろん numpy 配列化もできます。

# 抜き取りデータのbeta_1をnumpy配列化
idata_extract.beta1.to_numpy()

【実行結果】

🍀🍀🍀

⑤ おまけ：間引きする
テキストでは Win BUGS の MCMC サンプル生成時に、間引き「 $${\texttt{n.thin}}$$」を設定できることが書かれています（p.202）。
PyMC は MCMC サンプル時の間引き処理に対して消極的であり、サンプリング後に自由に間引いてください、的な印象です。

では サンプリング済みの idata に対する間引き処理を書いてみましょう。
２個飛ばしで 500 個のサンプルを残します。

# 参考：間引きthinningの方法
thinned_idata = idata.sel(draw=slice(None, None, 3))
thinned_idata.posterior

【実行結果】

draw が 500 に減りました（減少前：1500）。
chain 数は変化なしの３ですので、合計 1500 個の MCMC サンプルが残っています。
Coordinates の draw の右数字群をよく見ると 0 3 6 … のように、０から３つごとに番号が振られています。
draw のインデックス 0, 3, 6, … を（２個飛ばしで）選択した、という意味合いです。

【コード補足】
idata に対して $${\texttt{sel}}$$ メソッドでデータを選択（select）します。
$${\texttt{draw=slice(\cdot)}}$$ は draw 次元をスライス処理する司令です。
$${\texttt{(start=None, stop=None, step=3)}}$$ は「インデックスの頭（$${\texttt{start=None}}$$）から最後まで（$${\texttt{stop=None}}$$）を３個ごとに（$${\texttt{step=3}}$$）」みたいなことでしょう。

まとめ

今回のベイズ統計モデルをまとめます。

🔷 尤度
観測データはポアソン分布に従います。

🔷 リンク関数と線形予測子
ポアソン分布のパラメータ $${\lambda}$$ は線形予測子の指数関数です。
（対数リンク関数に相当）

🔷 事前分布
線形予測子内の２つのパラメータ $${\beta_1, \beta_2}$$ の事前分布は、平均０、標準誤差 100 の「ひらべったい正規分布」です。

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

次回は MCMC アルゴリズムの一種、ギブスサンプリングを検討します。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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の教科書です。
よかったらぜひ、お試しくださいませ。

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