「ベイズ統計モデリングによるデータ分析入門」をPythonとPyMCで写経 ～ Vol.7 PyMCコーディングの詳細

書籍の著者　馬場真哉 先生

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

PyMCのモデリングコードを詳しく見ていきます。
統計の典型例「平均値の差」のベイズ流推論も実践します！

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

はじめに

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

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

準備

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

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

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

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

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

# インポート

# 数値計算
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'

第６章 PyMCコーディングの詳細

ところどころテキストの 6.2 節～ 6.7 節を参照しながら、PyMC モデリングのコードを細かく見ていきます。

以下のコードは基本的にテキストの 6.8 節（最後の節）に沿って実装しています。

データの読み込みと外観の確認

テキストのデータを引用いたします。
２種類の銘柄 A、B のビールの売り上げデータです。
csv ファイルを pandas のデータフレーム形式で変数 file_beer_sales_ab に読み込みます。

# p.147 分析対象データの読み込み

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

【実行結果】
標本サイズ 200、変数 sales は 売り上げ量、beer_name は銘柄 A または B を示しています。

データを可視化しましょう。
テキスト 図 2.6.1 に相当するヒストグラムとKDE曲線の重ね描きチャートです。
seaborn ライブラリを利用します。
引数 hue に beer_name を指定して、ビール A とビール B 別に描きます。

# p.148 図2.6.1 2種類のビールの売上のヒストグラム

# ヒストグラムの描画
sns.histplot(data=file_beer_sales_ab, x='sales', hue='beer_name', bins=18,
             stat='density', edgecolor='white');
# KDEプロットの描画（曲線下を塗りつぶし）
sns.kdeplot(data=file_beer_sales_ab, x='sales', hue='beer_name', fill=True);

【実行結果】
ビール B の方が売り上げ量が大きいように見えます。
２つの銘柄の平均値には差があるのでは！？

銘柄別の要約統計量を確認します。

# A,B別の要約統計量

# Aの要約統計量
beer_a = (
    file_beer_sales_ab[file_beer_sales_ab['beer_name']=='A']['sales']
    .describe().rename('sales A')
)
# Bの要約統計量
beer_b = (
    file_beer_sales_ab[file_beer_sales_ab['beer_name']=='B']['sales']
    .describe().rename('sales B')
)
# データフレームに統合
pd.concat([beer_a, beer_b], axis=1).round(2)

【実行結果】
ビール B は平均値が大きいとともに、ばらつき（標準偏差）も大きいです。

では２変数の平均の差の $${t}$$ 検定へ！
ではなくて、ベイズ流の平均値差の推論へ進みます。

🚀🚀🚀

モデリング

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

変数 $${Y_A,\ Y_B}$$ はそれぞれビール A、ビール B の売り上げ量に関する尤度であり、正規分布に従うと仮定しています。
２つの正規分布の平均パラメータ変数 $${\mu_A,\ \mu_B}$$ は平均０、標準偏差 100000 の正規分布に、標準偏差パラメータ変数 $${\sigma_A,\ \sigma_B}$$ は標準偏差 100000 の半正規分布に従うとします。
ここまでの変数は確率分布と関連する「確率変数」です。

一方で、最後の変数 $${\text{diff}}$$ は確率変数ではなく、計算で求められる「決定論的変数」です。
$${\text{diff}}$$ はビール B の売り上げ量の期待値とビール A の売り上げ量の期待値の差であり、「平均値差」を表す変数です。

② PyMC のモデル設定
ビールの銘柄別に尤度を設定します。
テキストと同様に銘柄別にデータを分割します。

# p.148 データセットの準備

# 売上データをビールA、ビールBに分ける
beer_a = file_beer_sales_ab[file_beer_sales_ab['beer_name']=='A']['sales'].values
beer_b = file_beer_sales_ab[file_beer_sales_ab['beer_name']=='B']['sales'].values

【実行結果】なし

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

# モデリング

# coordsの設定
coords = {'id_a': range(len(beer_a)), 'id_b': range(len(beer_b))}

# モデルの定義
with pm.Model(coords=coords) as model:
    
    ## dataの設定
    # 目的変数
    sales_a = pm.Data('sales_a', value=beer_a, dims='id_a')
    sales_b = pm.Data('sales_b', value=beer_b, dims='id_b')

    ## 事前分布: 無情報事前分布的な分布
    mu_a = pm.Normal('mu_a', mu=0, sigma=1e5)            # ビールAの平均
    mu_b = pm.Normal('mu_b', mu=0, sigma=1e5)            # ビールBの平均
    sigma_a = pm.HalfNormal('sigma_a', sigma=1e5)        # ビールAの標準偏差
    sigma_b = pm.HalfNormal('sigma_b', sigma=1e5)        # ビールBの標準偏差

    ## 尤度
    obs_a = pm.Normal('obs_a', mu=mu_a, sigma=sigma_a, observed=sales_a,
                      dims='id_a')
    obs_b = pm.Normal('obs_b', mu=mu_b, sigma=sigma_b, observed=sales_b,
                      dims='id_b')

    ## 計算値
    diff = pm.Deterministic('diff', mu_b - mu_a)  # ビールAとビールBの売上平均の差

【実行結果】なし

テキスト 6.2 節～6.7 節のように、コードをじっくり眺めていきます。
◆ coords の設定

# coordsの設定
coords = {'id_a': range(len(beer_a)), 'id_b': range(len(beer_b))}

PyMCモデルの変数に「軸」を与えたいときに coords を Python の辞書型などで設定します。
カテゴリ変数のときにその威力を垣間見ることができるでしょう。
$${\mathtt{id\_a}}$$ を例にして深堀りします。
$${\mathtt{id\_a}}$$ は軸名称に相当します。
コロンの後続が $${\mathtt{id\_a}}$$ 軸の各要素に割り振る名称に相当します。
ここでは変数 beer_a の長さの整数連番（０～100）を設定しています。

◆ with pm.Model() as モデル名 の設定

# モデルの定義
with pm.Model(coords=coords) as model:

$${\mathtt{with\ pm.Model()\ as \sim:}}$$ は PyMC モデル設定のおまじないです。
PyMC の Model クラスで「モデル名」のオブジェクトを作ります的な感じです。
$${\mathtt{pm.Model()}}$$ の引数で先ほどの coords を指定します。

◆ pm.Data() の設定

sales_a = pm.Data('sales_a', value=beer_a, dims='id_a')

Stan は data ブロックで観測値などを定義するようです。
PyMC では観測値を改めてモデル内で定義しなくてもよいです。
この $${\mathtt{pm.Data()}}$$ で観測データ（学習データ）を指定しておくことの意義は、未知データによる予測時の操作がしやすくなること（メリット）にあるそうです。

◆ 確率変数（パラメータなど）

mu_a = pm.Normal('mu_a', mu=0, sigma=1e5)

この例はビール A の平均パラメータ mu_a の事前分布に関する設定です。
正規分布クラス $${\mathtt{Normal}}$$ でパラメータ mu と sigma の値を定数で設定しています。
次元・軸は設定していないです。
Stan の場合、事後分布を得たいパラメータは parameters ブロックで定義するようです。

確率分布クラスの設定内容は PyMC 公式サイトで確認しましょう。
こちらは正規分布 Normal ページの引用（ブラウザで翻訳）です。

◆ 確率変数（尤度など）

obs_a = pm.Normal('obs_a', mu=mu_a, sigma=sigma_a, observed=sales_a, dims='id_a')

観測値と紐付く確率変数です。
使われるケースの１つが尤度です。
パラメータを示す確率変数との違いは引数 $${\mathtt{observed}}$$ で観測データを指定します。
観測値で制約される変数なのです。
Stan の model ブロックで定義する変数と似ている感じです。

◆ 決定論的変数

diff = pm.Deterministic('diff', mu_b - mu_a)

四則演算などの計算だけで値が決まる「決定論的変数」は $${\mathtt{Deterministic}}$$ クラスで定義します。
この例は PyMC モデルで定義された mu_b と mu_a の差を設定しています。
Stan の generated quantities ブロックで定義する変数と似ている感じです。

モデルの数式ライク表示とグラフィカルモデル描画をします。

# モデルの表示
model

【実行結果】
決定論的変数 diff は計算式が出てきません（残念…）

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

【実行結果】
ビール A のラインとビール B のラインが分かれて 平均値差 diff でつながる「M 字形」になってます！
尤度の変数 obs_* と観測値 sales_* は軸 id_* の要素 100 個で成り立っています。
diff も複数個の要素をもつ（はず）ですが、未定義のため、軸無しです。

🚀🚀🚀

MCMC の実行

MCMCを実行しましょう。
NUTS サンプラーに nutpie を利用します。

%%time
# p.149 乱数の生成（MCMCの実行）
with model:
    idata = pm.sample(draws=1000, tune=1000, chains=4, random_seed=1,
                      nuts_sampler='nutpie')

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

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

# p.149 結果はこちら（要約統計量の表示）
var_names = ['mu_a', 'sigma_a', 'mu_b', 'sigma_b', 'diff']
az.summary(idata, var_names=var_names, hdi_prob=0.95)

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

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

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

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

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

🚀🚀🚀

平均値の差の評価

変数 diff の 事後分布を見てみます。
要約統計量を再掲します。最下行をご覧ください。

平均値差 diff は ビール B の売り上げ量 からビール A の売り上げ量を引いたものです。
ベイズ流の平均値差 diff は 事後平均 66.7（標準誤差 3.5）、95% HDI [ 59.9, 73.6 ]（０を含まない）です。
ビール B の売り上げ量の平均値はビール A を 66.7 ほど上回っていると言えそうです。

diff の事後分布を描画しましょう。
テキスト 図 2.6.2 に相当します。
arviz の plot_posterior 関数を利用します。

# p.150 図 2.6.2 平均値の差の事後分布
az.plot_posterior(idata, var_names=['diff'], hdi_prob=0.95, round_to=4);

【実行結果】
要約統計量で読んだ事後平均 mean と 95% HDI がビジュアルに確認できました。

平均値の差がある値よりも大きい確率を可視化しましょう。
事後分布の描画時に引数 ref_val でしきい値を設定します。
今回は平均値の差が 55 より大きい確率を見てみます。

# 平均値の差の事後分布 ：販売量55より大きい確率の表示
az.plot_posterior(idata, var_names=['diff'], hdi_prob=0.95, ref_val=55,
                  round_to=4);

【実行結果】
オレンジ垂線がしきい値 55 を示し、図中央のオレンジの％表示が「55より小さい確率」と「55より大きい確率」です。
「55 < 100.0%」の表示ではありますが、55より小さい確率が「0.1% < 55」ですので、平均値差が55より大きい確率は 99.9% 程度でしょう。

以上で 平均値差の評価を終わりにします。
面白かったですね。

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

シリーズの記事

次の記事

前の記事

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

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