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

書籍の著者　馬場真哉 先生

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

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

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

はじめに

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

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

準備

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

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

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

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

■ Stan フォルダ
Stan ファイルを格納するフォルダは、Jupyter Notebook ファイルと同一フォルダの「stan」フォルダです。

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

# インポート

# 数値計算
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の結果の評価【前編】

前編では、前回記事（第４章）と同じ「架空のビール売り上げデータ」を同じベイズ統計モデルでモデリングし、MCMCサンプルを詳しく操作します。

データの読み込み

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

# ファイルの読み込み
file_beer_sales_1 = pd.read_csv('./data/2-4-1-beer-sales-1.csv')
# 結果の表示
print('file_beer_sales_1.shape: ', file_beer_sales_1.shape)
file_beer_sales_1.head(3)

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

🔵🔵🔵

モデリング・MCMC実行

テキスト 5.2 節に相当します。
次のモデルを実装します。前回記事と同じモデルです。

Stan ファイル（Stan コード）を作成します。
テキストの Stan ファイル名 および Stan コードを引用いたします。
📑ファイル名：2-4-1-calc-mean-variance.stan
（前回記事と同じファイルです）

data {
    int N;                // サンプルサイズ
    vector[N] sales;      // データ
}
parameters {
    real mu;              // 平均
    real<lower=0> sigma;  // 標準偏差
}
model {
    // 平均mu, 標準偏差sigmaの正規分布に従ってデータが得られたと仮定
    for (i in 1:N) {
        sales[i] ~ normal(mu, sigma);
    }
}

Stan に渡す変数を定義します。

# データセットの準備
sample_size = len(file_beer_sales_1)
data_dict = {'N': sample_size,
             'sales': file_beer_sales_1['sales']}

【実行結果】なし

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

%%time
# モデルのコンパイル

# stanプログラムファイルのパス指定
current_dir = os.path.abspath(os.getcwd())
stan_file = os.path.join(current_dir, 'stan', '2-4-1-calc-mean-variance.stan')

# モデルオブジェクトの作成（exeの作成）
model = CmdStanModel(stan_file=stan_file)  # stanファイルを指定

【実行結果】
爆速で終了しています。

実は、前回記事と同じStanファイルを用いていて、かつ、すでにStanファイルや実行ファイル（exeファイル）が作成されているので、コンパイル実行時に「コンパイルをしないで」終了します。
つまり、すでにコンパイル済みの場合には、当該コンパイル済みファイルを活かす、ということです。

MCMCを実行します。
前回記事の sample 文から次の点を変更しています。

• sig_figs 引数の追加（出力の数値精度の設定）

%%time
fit = model.sample(
    data=data_dict,          # 対象データ
    seed=1,                  # 乱数の種
    chains=4,                # チェーン数
    iter_warmup=1000,        # バーンイン期間（１チェーンあたり）
    iter_sampling=1000,      # サンプリング数（１チェーンあたり）
    thin=1,                  # 間引き数（１なら間引き無し）
    sig_figs=18,             # 出力CSV等に適用する数値精度
)

【実行結果】
変数 fit にMCMC サンプル等のデータを格納しました。

🔵🔵🔵

MCMCサンプルを操作する

テキスト 5.3 節に相当します。
ArviZ ライブラリで MCMC サンプルの操作を行いたいので、 CmdStan 出力を InferenceData 形式に変換します。

# p.126 MCMCサンプルの抽出 # shape=(draw=1000, chain=4, [mu, sigma])
idata = az.from_cmdstanpy(posterior=fit, log_likelihood='lp__')
idata

【実行結果】

ちなみに CmdStan の出力を格納した fit の MCMC サンプルを pandas のデータフレームに変換したい場合は、次のコードで実現できます。

# 参考：データフレーム
fit.draws_pd()

【実行結果】
左の chain__ と iter__ が chain、draw の番号を持ち、右の mu、sigma がMCMCサンプルになっています。

では、テキストのコードを参考にして、InferenceData 形式の idata の調理法を学んでいきます。

① 型を調べる
idata の型を調べます。Python 標準の type 関数を利用します。

# p.126 クラス⇒変数の型を確認
type(idata)

【実行結果】
arviz の InferenceData 形式であることを確かめました。

② MCMC サンプルの次元を調べる
idata 配下のグループ「posterior」（ポステリア：事後分布）に MCMC サンプルが格納されているので、まずはこちらを可視化します。
ストレートに「idata.posterior」と打ち込みます。

# p.127 次元数,各々の名称
idata.posterior

【実行結果】
なんだか忙しく数字が並んでいます。
少々インタラクティブに動かせまして、下の図は「Data variables」の「mu」「sigma」の右のドラム缶アイコンをクリックして、MCMCサンプルの明細を表示しています。

データの次元は「Dimensions」にあるとおり、MCMC実行時に指定したチェーン４本 ⇒ chain: 4、ドロー 1000個 ⇒ draw 1000 となっています。
ちなみに書籍の第３次元「3」は、「Data variables」の３つの変数「mu」「sigma」「sigma_log__」と対応しています。

idata.posterior の型を確認します。

type(idata.posterior)

【実行結果】
先頭のライブラリ名が「xarray」に変わりました。
idata.posterior は Xarray ライブラリの Dataset 形式です。

Dataset 形式で idata.posterior の操作をしたい場合には、ぜひ、Xarray の公式サイトをご覧ください。

ちなみに Xarray.Dataset の流儀で次元数を確認してみましょう。

idata.posterior.dims

【実行結果】
chain 軸で 4、draw 軸 で 1000 という結果を確認できました。

③ パラメータmuの１回目の chain の最初の MCMC サンプルを調べる
事後分布グループ化 posterior の配下に変数 mu があります。
mu から sel メソッド（select）を使用して、最初のchain軸 = ラベル 0（Pythonは0始まり）と 最初のdraw軸 = ラベル 0 を抽出します。

# p.127 パラメータmuの１回目のチェーンで得られた最初のMCMCサンプル
idata.posterior.mu.sel(chain=0, draw=0)

【実行結果】
選んだ mu の MCMC サンプルの値は「103.90311849」です。
左のドラム缶アイコンのそばに書かれています。

④ パラメータmuの１回目の chain 全体の MCMC サンプルを調べる
パラメータ mu の最初の chain の 1000 個のMCMCサンプルを表示します。
③ と同様に sel メソッドを使います。

# p.127 パラメータmuの１回目のチェーンのMCMCサンプル shape=(1000,)
idata.posterior.mu.sel(chain=0)

【実行結果】
乱数感が半端ないですね！

⑤ パラメータmuの MCMC サンプル の個数（要素数）を調べる
最初の chain の要素数（= 1000）を表示します。

# p.128 パラメータmuの１回目のチェーンのMCMCサンプルの個数
idata.posterior.mu.sel(chain=0).size

【実行結果】

全 chain（４本の chain）の要素数（= 4000）を表示します。

# p.128 ４つのチェーンすべてのMCMCサンプルの個数
idata.posterior.mu.size

【実行結果】

MCMC サンプルの形状（shape）を表示します。

# p.128 ４つのチェーンがあるので、4Chain x 1000iter
idata.posterior.mu.shape

【実行結果】

🔵🔵🔵

事後分布の代表値を調べる

テキスト 5.4 節に相当します。
MCMC サンプルの平均値・中央値・2.5%点・97.5%点を調査します！

① 準備：平坦な numpy 配列に変換する
テキストは R 言語で「4000の要素をもつベクトル」化しています。
Python でベクトル型の変数を定義するのは大変そうなので、普段遣いの numpy 配列に変換します。
xarray.DataArray の stack メソッドで単一次元（つまり平坦）に変換し、さらに data 属性で配列化します。

# p.128 ベクトルにする⇒平坦なnumpy配列化
mu_mcmc_vec = idata.posterior.mu.stack(sample=('chain', 'draw')).data

numpy 配列になっていることを確認します。

type(mu_mcmc_vec)

【実行結果】

② 事後中央値を算出する
MCMC サンプルは事後分布からサンプリングされています。
ですので、MCMC サンプルの中央値は事後分布の中央値＝事後中央値です。
numpy.median で計算します。

# p.128 事後中央値
np.median(mu_mcmc_vec)

【実行結果】
パラメータ mu の事後中央値は 102.2 です。

変換前の idata のまま中央値を計算します。
xarray.DataArray の median メソッドを使います。

idata.posterior.mu.median().data

【実行結果】

PyMC 利用時に通常用いる arviz の summary 関数で事後分布の要約統計量を計算して、先ほどの中央値と一致することを確認しましょう。

az.summary(idata, var_names=['mu'], stat_focus='median', kind='stats',
           hdi_prob=0.95, round_to=6)

【実行結果】
中央値（median）は一致しています。

③ 事後期待値を算出する
MCMC サンプルの平均値は事後分布の平均値＝事後期待値です。
numpy.mean で計算します。

# p.129 事後期待値
mu_mcmc_vec.mean()

【実行結果】
パラメータ mu の事後期待値は 102.2 です。

変換前の idata のまま平均値を計算します。
xarray.DataArray の mean メソッドを使います。

idata.posterior.mu.mean().data

【実行結果】

arviz の summary 関数で事後分布の要約統計量を計算して、先ほどの平均値と一致することを確認しましょう。

az.summary(idata, var_names=['mu'], stat_focus='mean', kind='stats',
           hdi_prob=0.95, round_to=6)

【実行結果】
平均値（mean）は一致しています。

④ 95% ベイズ信用区間を算出する
MCMC サンプルの 2.5% 点、97.5% 点は事後分布の 95% ベイズ信用区間です。
numpy.quantile で計算します。

# p.129 95%ベイズ信用区間 ★fit.summary()の95%区間と違う値になった…
np.quantile(mu_mcmc_vec, [0.025, 0.975])

【実行結果】
パラメータ mu の 95% ベイズ信用区間は [98.6, 105.7] です。

変換前の idata のまま 2.5% 点・ 97.5% 点を計算します。
xarray.DataArray の quantile メソッドを使います。

idata.posterior.mu.quantile(q=[0.025, 0.975]).data

【実行結果】

CmdStan の出力 fit に対する summary メソッドで MCMC サンプルの要約統計量を表示することができます。

fit.summary(percentiles=[2.5, 50, 97.5], sig_figs=18)

【実行結果】

実は不可解な点がありまして…
numpy や xarray の quantile 関数で算出した 95% 信用区間と、cmdstanpy の summary メソッドで算出した 95% 信用区間が合わないのです…
生成AIは「分位数の補間アルゴリズムが一致していないのでは？」と情報提供してくれましたが、アルゴリズムの検証はできておらず、執筆時点で不一致理由は不明です…

🔵🔵🔵

トレースプロットの描画

テキスト 5.5 節に相当します。
arviz の plot_trace 関数でトレースプロットを描画します。

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

【実行結果】
右側のチャートがテキスト 図 2.5.1 に相当します。

【コード補足】
引数「backend_kwargs={'tight_layout': True}」は Python の可視化ライブラリ matplotlib に向けたメッセージです。
１つのチャートに複数のグラフが含まれるとき、グラフ間の間隔を自動調整して「見やすくするため」のおまじないが「backend_kwargs={'tight_layout': True}」なのです。

🔵🔵🔵

arviz を使わないで事後分布を可視化

テキスト 5.6 節に相当します。
seaborn ライブラリの kdeplot で MCMC サンプルのカーネル密度推定⇒事後分布を描画します。

# p.130 図2.5.2 seabornで事後分布のKDEを描画
# sns.kdeplot(idata.posterior.mu.stack(sample=('chain', 'draw')));
sns.kdeplot(mu_mcmc_vec)
plt.xlabel('mu', fontsize=12);

【実行結果】

２行目のコメント化されたコードでも描画できます。
変換前の idata を使って描画できることが分かります。

🔵🔵🔵

arviz で事後分布を可視化

テキスト 5.7 節、5.8 節、5.10 節に相当します。
テキストの bayesplot のように、MCMC サンプルの可視化をサポートする Python のライブラリは「ArviZ」です。
MCMC の結果を格納した idata をそのまま適用できます。

① ヒストグラム
パラメータ $${\mu,\ \sigma}$$ のヒストグラムを plot_dist 関数で描画します。
引数 kind に 'hist'（ヒストグラム）を設定します。
１命令で１変数に対応するようですので、２変数を並べて描画するためにちょっと長めのコードです。

# p.131 ヒストグラム
fig, ax = plt.subplots(2, 1, figsize=(10, 5), tight_layout=True)
az.plot_dist(idata.posterior.mu, kind='hist', ax=ax[0])
az.plot_dist(idata.posterior.sigma, kind='hist', ax=ax[1])
ax[0].set(title='パラメータ $\\mu$', xlabel='$\\mu$', ylabel='density')
ax[1].set(title='パラメータ $\\sigma$', xlabel='$\\sigma$', ylabel='density');

【実行結果】

② カーネル密度推定プロット
plot_density 関数で複数パラメータを一気に描画します。

# p.131 カーネル密度推定
az.plot_density(idata, hdi_prob=1.00, shade=0.3);

【実行結果】

③ 事後分布プロット
事後分布のカーネル密度推定を可視化する場合、②の plot_density よりも、こちらの plot_posterior 関数を使う方が一般的だと思います。
早速描画しましょう。

# 事後分布プロットの描画
az.plot_posterior(idata, hdi_prob=0.95);

【実行結果】

④ トレースプロット
plot_trace 関数でトレースプロットを描画します。
左側に事後分布プロット（カーネル密度推定）、右側にトレースプロット（ゲジゲジ）を表示します。

# p.131, 132 図2.5.3 事後分布とトレースプロットの図示
az.plot_trace(idata, compact=True, figsize=(8, 4),
              backend_kwargs={'tight_layout': True});

【実行結果】
テキスト 図 2.5.3 相当の図になりました。

⑤ 事後分布の 95% HDI 区間
plot_forest 関数でフォレストプロットを描画します。
テキスト 図 2.5.4 に相当します。

# p.132 図2.5.4 範囲比較1 事後分布の範囲を比較
az.plot_forest(idata, combined=True, hdi_prob=0.95);

【実行結果】
横棒がパラメータの事後分布の 95% HDI の区間を示しています。
「多くのパラメータの比較を直感的に見たい」「95% HDI 区間が０を含んでいるか見たい」などのときに威力を発揮すると思います！

⑥ フォレストプロットとカーネル密度推定を重ね描き
テキスト 図 2.5.5 で紹介されているチャートです。
少々無理やり感は否めませんが、やってみました。

# p.133 図2.5.5 範囲比較2 カーネル密度を重ね描き
az.plot_forest(idata, kind='ridgeplot', combined=True, hdi_prob=0.6,
               ridgeplot_truncate=False,
               ridgeplot_quantiles=[.005, .5, .995],
               colors='lightblue',
               figsize=(6, 4));

【実行結果】

⑦ MCMC サンプルの自己相関
plot_autocorr 関数でMCMCサンプルの自己相関プロット（コレログラム）を描画します。
テキスト 図 2.5.6 に相当します。

# p.134 図2.5.6 コレログラム
az.plot_autocorr(idata, max_lag=20, figsize=(8, 5),
                 backend_kwargs={'tight_layout': True});

【実行結果】
青い縦棒が自己相関係数の値、グレイの帯が 95% 信頼区間です。
このチャートの場合、青い縦棒がグレイの帯の内側に包含されているので、自己相関は無いだろう、と見ます。

以上で MCMC サンプルを用いた事後分布の可視化を終わりにします。
面白かったですね。

５章のもう一つの論点「事後予測チェック」は次の記事で取り組みます。

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

シリーズの記事

次の記事

前の記事

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

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