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

書籍の著者　馬場真哉 先生

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

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

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

はじめに

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

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

準備

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

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

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

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

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

# インポート

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

# ベイズ統計モデリング
import pymc as pm                   # pymc
import arviz as az                  # 分析・可視化

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

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

PyMC のモデルを書きましょう。

# モデリング

# coordsの設定
coords = {'id': file_beer_sales_1.index.values}  # id: 0始まりの連番

# モデルの定義
with pm.Model(coords=coords) as model:
    
    ## dataの設定
    # 目的変数
    Y = pm.Data('Y', value=file_beer_sales_1['sales'].values, dims='id')

    ## 事前分布: 無情報事前分布的な分布
    mu = pm.Normal('mu', mu=0, sigma=1e5)
    sigma = pm.HalfNormal('sigma', sigma=1e5)

    ## 尤度関数
    obs = pm.Normal('obs', mu=mu, sigma=sigma, observed=Y, dims='id')

【実行結果】なし

モデルの外観を眺めます。

# モデルの表示
model_norm

【実行結果】

【補足】
変数 pred は、「モデルを用いた amimal_num の予測値」（事後予測分布からの乱数）に相当します。
PyMC の場合、事後予測分布からの乱数を取得するには通常、別の方法を使いますが、ここは一旦、テキストに寄せて pred を算出します。

モデルをグラフィカルモデル形式で可視化します。

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

【実行結果】

MCMCを実行しましょう。

%%time
# MCMCの実行
with model_norm:
    idata_norm = pm.sample(draws=1000, tune=1000, chains=4, random_seed=1,
                           nuts_sampler='nutpie')

【実行結果】
Divergences（ダイバージェンス）は０件です。

収束の確認をします。
MCMC サンプルの要約表を表示します。

# 要約統計量の表示
var_names = ['mu', 'sigma']
az.summary(idata_norm, var_names=var_names, hdi_prob=0.95)

【実行結果】
$${\widehat{R}}$$（r_hat）、有効サンプル数（ess_bulk、ess_tail）に問題はなさそうです。

トレースプロットを描画します。

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

【実行結果】
右側のチャートがゲジゲジしています。

以上のチェックに基づいて、収束していると考えましょう。

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

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

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

PyMC のモデルを書きましょう。

# モデリング

# coordsの設定
coords = {'id': animal_num.index.values}

# モデルの定義
with pm.Model(coords=coords) as model_pois:
    
    ## dataの定義
    # 目的変数
    Y = pm.Data('Y', value=animal_num['animal_num'].values, dims='id')

    ## 事前分布: sigmaの値は調整しました
    lam = pm.HalfNormal('lam', sigma=100)

    ## 尤度関数
    obs = pm.Poisson('obs', mu=lam, observed=Y, dims='id')

    ## 事後予測
    pred = pm.Poisson('pred', mu=lam, dims='id')

【実行結果】なし

モデルの外観を眺めます。

# モデルの表示
model_pois

【実行結果】

【補足】
正規分布モデルと同じ理由で変数 pred を付け足しています。

モデルをグラフィカルモデル形式で可視化します。

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

【実行結果】

MCMCを実行しましょう。

%%time
# MCMCの実行 ※rhat<=1.01になるように、draw数を増やしました
with model_pois:
    idata_pois = pm.sample(draws=2000, tune=1000, chains=4, random_seed=1)

【実行結果】
Divergences（ダイバージェンス）は０件です。

【コードの補足】
推定パラメータに離散分布を設定した場合、NUTS サンプラーを使うことができず、PyMC デフォルトでは メトロポリス・ヘイスティング法サンプラーが適用されます。
今回のモデルでは、正規分布をポアソン分布の平均パラメータ $${\lambda}$$ にしており、「あれ？」と感じるかも知れません。
実は、pred にポアソン分布を設定したことで、メトロポリス・ヘイスティング法サンプラーを適用しています（適用になってしまいました）。

収束の確認をします。
MCMC サンプルの要約表を表示します。

# 要約統計量の表示
var_names = ['lam']
az.summary(idata_pois, var_names=var_names)

【実行結果】
$${\widehat{R}}$$（r_hat）、有効サンプル数（ess_bulk、ess_tail）に問題はなさそうです。

トレースプロットを描画します。

# トレースプロットの描画
az.plot_trace(idata_pois, var_names=var_names, compact=False, figsize=(8, 3),
              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, id=range(20)).data

【実行結果】
３次元目の軸名 id は、PyMCモデリングの際に coords で指定したアノ id です。

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

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

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

# ポアソン分布を仮定したモデル
print(idata_pois.posterior.pred.sel(chain=0, draw=0, id=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_pois.posterior.pred[:, :, :5], kind='hist', bins=15,
                  ec='white', round_to=3, ax=ax[1:]);

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

🚀🚀🚀

事後予測チェック２：PyMCっぽく

PyMC および ArviZ には事後予測分布からの乱数生成や事後予測チェックのための機能が備わっているので、わざわざ pred をモデルに組み込んで MCMC サンプルを生成する必要がありません。

① 正規分布モデルで事後予測分布からの乱数を生成する
pymc の sample_posterior_predictive 関数で事後予測サンプリングを行います。

# 事後予測チェック 正規分布モデル

# 事後予測サンプリング
with model_norm:
    pm.sample_posterior_predictive(
        idata_norm, extend_inferencedata=True, random_seed=123
    )

【実行結果】
進捗バーが 100% になり、サンプリングが完了しました。

【コード補足】
with 構文を用いて 正規分布モデル model_norm の下でサンプリングを行います。
引数は idata 本体、idata に事後予測サンプリング結果を格納するか、乱数シードです。

idata のグループ posterior_predictive（事後予測）にサンプルが追加されています。
確認してみましょう。

idata_norm.posterior_predictive

【実行結果】

Data variables 配下に「obs」の名で事後予測サンプルが格納されています。
次元は dimensions に記載のとおり (chain: 4, draw: 1000, id: 200) です。
観測値の id １つ１つに $${4 \times 1000 = 4000}$$ の事後予測サンプルデータがあるということです。

事後予測サンプルを使って事後予測チェックプロットを描画します。
arviz の plot_ppc 関数（Posterior Predictive/Prior Check）を利用します。

# 事後予測チェックプロットの描画
az.plot_ppc(idata_norm, num_pp_samples=100);

【実行結果】

黒線が観測値の分布、オレンジ点線が事後予測平均の分布、細い青線のかたまりは観測値１件ごとの事後予測の分布です。
観測値の分布と事後予測の分布は全然違います！
今回のデータに正規分布モデルは合っていないのです！

【コード補足】
引数 num_pp_samples=100 は細い青線で描画する事後予測の個数を 100 に制限しています。
このモデルでは事後予測のデータの個数は 4000 です。
4000 全部を描画するとチャート内が騒がしくなるので、描画数を制限（間引き）しています。

② ポアソン分布モデルで事後予測分布からの乱数を生成する
ポアソン分布モデルの事後予測チェックプロットを一気に描画します。

# 事後予測チェック ポアソン分布モデル

# 事後予測サンプリング
with model_pois:
    pm.sample_posterior_predictive(
        idata_pois, extend_inferencedata=True, random_seed=123
    )

# 事後予測チェックプロットの描画
az.plot_ppc(idata_pois, num_pp_samples=100);

【実行結果】

ポアソン分布モデルでは黒線の観測値の分布とオレンジ点線の事後予測平均の分布が似ています。
ポアソン分布モデルのほうが正規分布モデルよりも適していることを、直感的に確認できました。

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

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

シリーズの記事

次の記事

前の記事

Stan版

目次

ブログの紹介

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

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