「ベイズ統計モデリングによるデータ分析入門」をPythonとStanで写経～ Vol.6 MCMCの結果の評価【後編】

ネイピア DS

2026年2月27日 11:05

書籍の著者　馬場真哉先生

この記事は、書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」第２部第５章「MCMCの結果の評価」の Python 写経活動記録です。

CmdStanPyでMCMCを実行した結果得られる「MCMCサンプル」を用いて、「事後予測」の確認を行います。

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

はじめに

このブログシリーズは書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」（講談社、「テキスト」と呼びます）の Python 写経です。

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

準備

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

5.11 事後予測チェックの対象となるデータとモデル
5.13 事後予測チェックのためのMCMCの実行
5.14 bayesplotによる事後予測チェック

ベイズ統計モデリングの理論面が気になった場合には、ぜひテキストの第１部をご覧いただき、本記事との繋がりをご確認下さいませ。

■ コード記述法
Jupyter Notebook 形式でコードを記述します。

■ 利用データ
テキスト・サポートサイトのデータファイルを引用しています。
Jupyter Notebook ファイルと同一フォルダ内の「data」フォルダにデータファイルを格納しています。

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

# インポート

# 数値計算
import numpy as np
import pandas as pd

# ベイズ統計モデリング
from cmdstanpy import CmdStanModel  # stan
import arviz as az                  # 分析・可視化

# ユーティリティ
import os

# 可視化
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_theme()                     # ggplot風のスタイル
plt.rcParams['font.family'] = 'Meiryo'

第５章 MCMCの結果の評価【後編】

後編では「ある動物の個体発見数データ」に基づいて、２つのベイズ統計モデルをモデリングし、事後予測チェックでモデルの良し悪しを確認します。

データの読み込み

テキストの 5.11 節に相当します。
テキストのデータを引用いたします。
csv ファイルを pandas のデータフレーム形式で変数 animal_num に読み込みます。

# p.135 分析対象データの読み込み
animal_num = pd.read_csv('./data/2-5-1-animal-num.csv')
print('animal_num.shape :', animal_num.shape)
animal_num.head(3)

【実行結果】
標本サイズ 200、変数が「animal_num」のみのデータです。

データを可視化しましょう。
ヒストグラムです。
Python の代表的な可視化ライブラリ matplotlib を利用します。

# ヒストグラム（実態は棒グラフ）の描画
plt.hist(animal_num['animal_num'], bins=np.arange(-0.5, 5))
plt.xlabel('発見個体数')
plt.ylabel('確率');

【実行結果】
データが離散値であり、この図のような確率質量関数の形状は…
ポアソン分布が候補かもです！

少々実験をしましょう。
データを当てはめてポアソン分布の平均パラメータ $${\lambda}$$ を推定し、確率質量関数とヒストグラムを重ね書きしてみます。
scipy.stats の fit 関数で確率分布のパラメータ推定（最尤法）を行います！

# 追加インポート
import scipy.stats as stats

# データをポアソン分布に当てはめて平均パラメータを推定
res = stats.fit(stats.poisson, animal_num['animal_num'], bounds={'mu': (0, 10)})

# ヒストグラムの描画
plt.hist(animal_num['animal_num'], density=True, bins=np.arange(-0.5, 5))

# 当てはめしたポアソン分布の確率質量関数の描画
plt.plot(range(5), stats.poisson.pmf(range(5), *res.params), '-o',
         color='tab:red', label=f'$\\lambda$={res.params.mu:.2f}のポアソン分布')

# 修飾
plt.xlabel('発見個体数')
plt.ylabel('確率')
plt.legend();

【実行結果】
平均パラメータ $${1.01}$$ のポアソン分布の確率質量関数（赤い折れ線）を重ねています。
データはポアソン分布に従っているように感じます。

データの要約統計量を確認します。

# 要約統計量の表示
animal_num.describe().round(2)

【実行結果】
平均値は $${1.01}$$ です。
推定したポアソン分布の平均値パラメータと同じですね！

🔵🔵🔵

モデリング・MCMC実行（正規分布モデル）

テキスト 5.11 節、5.13 節に相当します。
テキストが比べる２つのモデルのうち、こちらでは「正規分布モデル」を実装します。
データがポアソン分布に従っているだろう、と先ほど確認しました。
ですので、この正規分布モデルは「良くないモデル」の役目を担います。

次のモデルを実装します。

$$
\text{animal\_num} \sim \text{Normal}\ (\mu,\ \sigma^2)
$$

Stan ファイル（Stan コード）を作成します。
テキストの Stan ファイル名および Stan コードを引用いたします。
📑ファイル名：2-5-1-normal-dist.stan

data {
    int N;                  // サンプルサイズ
    vector[N] animal_num;   // データ
}

parameters {
    real<lower=0> mu;       // 平均
    real<lower=0> sigma;    // 標準偏差
}

model {
    // 平均 mu, 標準偏差 sigma の正規分布
    animal_num ~ normal(mu, sigma);
}

generated quantities {
    // 事後予測分布を得る
    vector[N] pred;
    for (i in 1:N) {
        pred[i] = normal_rng(mu, sigma);
    }
}

【コードの補足】
テキストのとおり、generated quantities ブロックに pred の乱数取得コードを書いています。

モデルのコンパイルを実行します。

%%time
# モデルのコンパイル：正規分布

# stanプログラムファイルのパス指定
current_dir = os.path.abspath(os.getcwd())
stan_file = os.path.join(current_dir, 'stan', '2-5-1-normal-dist.stan')
# モデルオブジェクトの作成（exeの作成）
model_norm = CmdStanModel(stan_file=stan_file)

【実行結果】（右側のファイルパスは記載省略）

MCMCを実行します。
sample メソッドの引数 $${\mathtt{sig\_figs}}$$ には MCMC サンプルデータを CSV ファイル出力する際の数値精度を指定します。
この CSV ファイルは arviz の idata 形式に変換する際に利用されます。

%%time
# p.137 MCMCの実行

## データセットのまとめ
# サンプルサイズ
sample_size = len(animal_num)
# 辞書にまとめる
data_dict = dict(N=sample_size, animal_num=animal_num['animal_num'])

## MCMCの実行：正規分布モデル
fit_norm = model_norm.sample(
    data=data_dict,          # 対象データ
    seed=1,                  # 乱数の種
    sig_figs=18,             # 出力CSV等に適用する数値精度
)

【実行結果】

英文が騒がしいです…

WARNING - Non-fatal error during sampling
> ワーニング：サンプリング中で致命的ではないエラーが発生、とのこと

Location parameter is inf, but must be finite!
> 位置パラメータ（データの中心位置を表す値）が無限大（inf）になっているが、有限な実数値でなければならない

Geminiによると「正規分布の対数確率密度関数 (normal_lpdf)を計算する際に、その分布の中心（ロケーションパラメータ、平均値）が無限大（inf）になってしまい、計算できないというエラー」が３回発生しています。

…
結果的に正規分布モデルは適切ではないので、このワーニングに目を瞑っておきます。

🍀🍀🍀

収束の確認をします。
診断サマリーを表示しましょう。

# MCMCの診断
print(fit_norm.diagnose())

【実行結果】
５つの項目をチェックして「no preblems …」（問題はない）という診断結果が表示されました。

要約統計量を確認しましょう。

# 要約統計量の表示
fit_norm.summary()

【実行結果】
上方に表示された mu と sigma の $${\hat{R}}$$ 値、および、有効サンプル数（N_Eff）は特段、問題はなさそうです。

トレースプロットを確認します。
まず arvis の from_cmdstanpy 関数を利用して fit を arviz の idata 形式に変換します。

# arvizのidataに変換
idata_norm = az.from_cmdstanpy(posterior=fit_norm, log_likelihood='lp__')

この idata 形式のデータは、事後予測チェックでも使用いたします。

ではトレースプロットを描画しましょう。

# トレースプロットの描画
var_names = ['mu', 'sigma']
az.plot_trace(idata_norm, var_names=var_names, compact=False, figsize=(8, 4),
              backend_kwargs={'tight_layout': True});

【実行結果】
右側のチャートは満遍なくゲジゲジしています。収束を感じます。
左側のチャートでは４本のチェーンがほぼ同じ分布になったことを確認できます。

🔵🔵🔵

モデリング・MCMC実行（ポアソンモデル）

テキスト 5.11 節、5.13 節に相当します。
テキストが比べる２つのモデルのうち、こちらでは「ポアソン分布モデル」を実装します。

次のモデルを実装します。

$$
\text{animal\_num} \sim \text{Poisson}\ (\lambda)
$$

Stan ファイル（Stan コード）を作成します。
テキストの Stan ファイル名および Stan コードを引用いたします。
📑ファイル名：2-5-2-poisson-dist.stan

data {
    int N;                    // サンプルサイズ
    array[N] int animal_num;  // int animal_num[N]の書き方はStan2.33.0で廃止された
}

parameters {
    real<lower=0> lambda;     // 平均
}

model {
    // 強度lambdaのポアソン分布
    animal_num ~ poisson(lambda);
}

generated quantities {
    // 事後予測分布を得る
    array[N] int pred;        // int pred[N]の書き方はStan2.33.0で廃止された
    for (i in 1:N) {
        pred[i] = poisson_rng(lambda);
    }
}

【コードの補足】
テキストのコードではエラーが発生します。
Stan 2.33.0 から、要素数 N の配列の定義で $${\mathtt{int\ animal\_num[N];}}$$ のように $${\mathtt{int}}$$ 始まりで定義することができなくなりました。
配列としての定義 $${\mathtt{array[N]\ int\ animal\_num;}}$$ のように書く必要があります。

モデルのコンパイルを実行します。

%%time
# モデルのコンパイル：ポアソン分布 ★stanコードで廃止エラーが発生

# stanプログラムファイルのパス指定
current_dir = os.path.abspath(os.getcwd())
stan_file2 = os.path.join(current_dir, 'stan', '2-5-2-poisson-dist.stan')
# モデルオブジェクトの作成（exeの作成）
model_poi = CmdStanModel(stan_file=stan_file2)

【実行結果】（右側のファイルパスは記載省略）

MCMCを実行します。

## MCMCの実行：ポアソン分布モデル
fit_poi = model_poi.sample(
    data=data_dict,          # 対象データ
    seed=1,                  # 乱数の種
    sig_figs=18,             # 出力CSV等に適用する数値精度
)

【実行結果】
ポアソン分布モデルはワーニングもなく、無事に実行できました！

🍀🍀🍀

収束の確認をします。
診断サマリーを表示しましょう。

# MCMCの診断
print(fit_poi.diagnose())

【実行結果】
５つの項目をチェックして「no preblems …」（問題はない）という診断結果が表示されました。

要約統計量を確認しましょう。

# 要約統計量の表示
fit_poi.summary()

【実行結果】
上方に表示された lambda の $${\hat{R}}$$ 値、および、有効サンプル数（N_Eff）は特段、問題はなさそうです。

トレースプロットを確認します。
まず arvis の from_cmdstanpy 関数を利用して fit を arviz の idata 形式に変換します。

# arvizのidataに変換
idata_poi = az.from_cmdstanpy(posterior=fit_poi, log_likelihood='lp__')

この idata 形式のデータは、事後予測チェックでも使用いたします。

ではトレースプロットを描画しましょう。

# トレースプロットの描画
var_names = ['lambda']
az.plot_trace(idata_poi, var_names=var_names, compact=False, figsize=(8, 2.5),
              backend_kwargs={'tight_layout': True});

🔵🔵🔵

事後予測チェック

テキスト 5.14 節に相当します。
MCMC サンプルの中の pred を利用して、テキストにならって事後予測チェックを行います。

① pred のサンプルサイズ
正規分布モデルについて、idata → posterior → pred でアクセスして形状 shape を調べます。

# p.138 MCMCサンプルの個数 shape=(chains=4, draws=1000, サンプルサイズ=200)
idata_norm.posterior.pred.shape

【実行結果】
chains=4, draws=1000, サンプルサイズ=200 の３次元データです。

正規分布モデルの最初のチェーン・最初のドロー（MCMCサンプル）・サンプルサイズの冒頭 20 個の pred の値を調べてみます。

idata_norm.posterior.pred.sel(chain=0, draw=0, pred_dim_0=range(20)).data

【実行結果】
３次元目の軸 pred_dim_0 は、観測データに対応するものであり、名前は変換時に自動設定されたようです。

② pred のMCMC サンプルを１セット取り出す
２つのモデルについて、最初のチェーン・最初のドロー（MCMCサンプル）・サンプルサイズの冒頭 20 個の pred の値を調べてみます。
sel メソッド（select）を利用します。

# p.138 予測分布のサンプルを確認 chain0、draw0、sampleの先頭20を表示

# 正規分布を仮定したモデル
print(idata_norm.posterior.pred.sel(chain=0, draw=0, pred_dim_0=range(20)).data)

# ポアソン分布を仮定したモデル
print(idata_poi.posterior.pred.sel(chain=0, draw=0, pred_dim_0=range(20)).data)

【実行結果】
１つ目の正規分布モデルは連続値、２つめのポアソン分布モデルは離散値になっています。

③ 正規分布モデルの事後予測分布を可視化する
arviz の plot_posterior（事後分布プロット）を利用して、animal_num の観測値と事後予測サンプル（pred）の分布をヒストグラム形式で描画します。
テキスト図 2.5.7 に相当します。

# p.139 図2.5.7 正規分布を仮定したモデルの事後予測チェック by arviz

# 描画領域の設定
fig, ax = plt.subplots(2, 3, figsize=(10, 6), tight_layout=True)
# axesを平坦な配列に変換
ax = ax.flat
# 観測データの事後予測プロット ※animal_numを辞書型に変換して与える
az.plot_posterior(dict(y=animal_num), kind='hist', bins=10, ec='white',
                  color='green', round_to=3, ax=ax[0])
# 事後予測分布データの事後予測プロット
az.plot_posterior(idata_norm.posterior.pred[:, :, :5], kind='hist', bins=20,
                  ec='white', round_to=3, ax=ax[1:]);

【実行結果】

左上の緑のヒストグラムは観測値です。データ確認時に見たヒストグラムと同様に飛び飛びの値＝離散値になっています。
その他の青いヒストグラムは pred の冒頭 5 つの id 軸の MCMC サンプルです。
観測値の分布と事後予測分布が合っておらず、正規分布モデルが適切でないことが理解できます。
テキストの記載どおり「このモデルを信じて将来予測などを行うのは問題」だと分かります。

【コードの補足】
plot_posterior 関数に観測値 Y と MCMCサンプル pred を同時に与えて描画する方法が見つからなかったので、「観測値の描画」と「MCMCサンプルの描画」を別々に書きました。
グラフの描画領域は matplotlib の axes で 2 行 3 列に設定しており、最初の axes に観測値のヒストグラム、残りの axes にMCMCサンプルのヒストグラムを配置しています。

④ ポアソン分布モデルの事後予測分布を可視化する
正規分布モデルと同様に、arviz の plot_posterior（事後分布プロット）を利用して、animal_num の観測値と事後予測サンプル（pred）の分布をヒストグラム形式で描画します。
テキスト図 2.5.8 に相当します。

# p.140 図2.5.8 ポアソン分布を仮定したモデルの事後予測チェック by arviz

# 描画領域の設定
fig, ax = plt.subplots(2, 3, figsize=(10, 6), tight_layout=True)
# axesを平坦な配列に変換
ax = ax.flat
# 観測データの事後予測プロット ※animal_numを辞書型に変換して与える
az.plot_posterior(dict(y=animal_num), kind='hist', bins=10, ec='white',
                  color='green', round_to=3, ax=ax[0])
# 事後予測分布データの事後予測プロット
az.plot_posterior(idata_poi.posterior.pred[:, :, :5], kind='hist', bins=15,
                  ec='white', round_to=3, ax=ax[1:]);

【実行結果】
観測値も事後予測も離散分布になっています！
今回のデータにはポアソン分布モデルの方が適している感じがします。

以上で事後予測分布の可視化を終わりにします。
面白かったですね。

今回の記事は以上です。
楽しかったですね！

シリーズの記事

次の記事

前の記事

PyMC版

ブログの紹介

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

１．のんびり統計

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

２．統計・データ分析とつながる

シリーズ「統計・データ分析とつながる」は、統計・データ分析との「つながり」を発掘して、コラム風に仕立てたブログシリーズです。
生成 AI の力を借りながら、統計・データ分析の入り口をイメージして、自由気ままに書きました。
たとえば…
・日常生活と統計のつながり
・統計検定２級からその先へのつながり
気楽にお読みいただけたら嬉しいです🍀

３．実験！たのしいベイズモデリング１＆２を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の教科書です。
よかったらぜひ、お試しくださいませ。

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

いいなと思ったら応援しよう！

応援ありがとうございます。これからもがんばって記事を作成します！