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

ネイピア DS

2026年1月23日 11:49

書籍の著者　馬場真哉先生

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

書籍のR・StanコードをPython・PyMC化して、ベイズ統計モデリングをスタートします！

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

はじめに

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

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

準備

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

4.2 Stanのインストール
4.7 Stan: Stanファイルの実装例
4.10 R: CSVファイルから分析対象となるデータを読み込む
4.12 R: Stanと連携してMCMCを実行する
4.13 R: 推定結果を確認する
4.14 R: 収束の確認

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

■ コード記述法
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とArviZのバージョン確認
PyMC はベイズ統計モデリングライブラリです。
ArviZ はベイズ統計モデリングなどの可視化ライブラリです。
記事で用いる両ライブラリのバージョン確認をします。

# バージョン情報
print('PyMC', pm.__version__)
print('ArviZ', az.__version__)

【実行結果】

$$
\mathtt{PyMC \quad 5.23.0} \\
\mathtt{ArviZ\ \ 0.21.0}
$$

第４章 PyMCの基本

PyMCのインストール

テキスト 4.2 節に相当します。
PyMC 公式サイトの情報を参考にして、PyMCをインストールします。

公式サイトは Anaconda または Miniforge を使用して環境構築することを勧めているようです。

公式サイトのインストールページを引用します。ご参考まで。
（ブラウザの機能で日本語翻訳表示しています）

PyMC 公式サイトより引用（https://www.pymc.io/projects/docs/en/latest/installation.html）

🚀🚀🚀

データの読み込み

テキスト 4.10 節に相当します。
テキストの「架空のビール売り上げデータ」を引用いたします。
csv ファイルを pandas のデータフレーム形式で変数 file_beer_sales_1 に読み込みます。
read_csv 関数を利用します。

# p.117 分析対象のデータ

# ファイルの読み込み
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」のみのデータです。

データを可視化しましょう。
１変数ですので、ひとまずヒストグラムを描画します。
pandas の plot メソッドを利用します。

# ヒストグラムの描画 by pandas plot
file_beer_sales_1.plot.hist();

【実行結果】
平均 100 前後、60 ～ 140 の範囲で正規分布に似た形状をしています。

要約統計量を見ておきます。
pandas の describe メソッドを利用します。

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

【実行結果】
平均 102.18、55.71 ～ 148.03 の範囲の売り上げです。

🚀🚀🚀

いきなりモデリング

テキスト 4.7 節に相当します。
今回のモデリングで留意したい Stan と PyMC の相違点を記載します。

ファイルを別個に作成する必要がない
Stanの場合、Stanコードを記述したファイルを R や Python の実行コードファイルとは別に作成しますが、PyMCは別ファイル化する必要はありません。
Jupyter Notebook のコードセルで PyMC のモデルを記述します。
ブロックに分けて記述する必要がない
Stanの場合、data ブロック、parametars ブロック、model ブロックなどのブロックに分けてモデルを記述する必要があります。
しかし、PyMCではブロック分けする必要がありません。
モデルを読みやすくするには、コードやコメントの書き方で工夫します。
パラメータの事前分布の記述を省略しない
Stanの場合、パラメータの事前分布の記述を省略できます。
しかし、PyMCは明示する必要があります。
今回のモデルでは十分に範囲の広い正規分布を設定します。

以上のことを念頭に置いて、次のモデルを実装しましょう。

$$
\begin{array}{lll}
Y \sim \text{Normal}(\mu,\ \sigma^2) &\text{尤度（テキストでは sales）} \\
\mu \sim \text{Normal}(0,\ (1e5)^2) &\text{平均パラメータの事前分布} \\
\sigma \sim \text{HalfNormal}((1e5)^2) &\text{標準偏差パラメータの事前分布} \\
\end{array}
$$

モデルを噛み砕いてみます。

尤度
売り上げデータ $${Y}$$（テキストの $${\mathtt{sales}}$$）は平均パラメータ $${\mu}$$、標準偏差パラメータ $${\sigma}$$ の正規分布に従います。
MCMCサンプリングによって、パラメータ $${\mu,\ \sigma}$$ の事後分布を推定します。
パラメータの事前分布
ベイズ統計はパラメータを「確率変数」として扱います！
平均パラメータ $${\mu}$$ は平均 $${0}$$、標準偏差 $${100000}$$ の正規分布に従います。
標準偏差パラメータ $${\sigma}$$ は標準偏差 $${100000}$$ の半正規分布に従います。

それでは 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')

【実行結果】なし

【コード補足】

モデル構文
この記事では、PyMCのモデル記述に Python の with 構文を用います。
$${\mathtt{with\ pm.Model()\ as\ モデル名:}}$$ の下にモデルを記述します。
ちなみに今回のモデル名は $${\mathtt{model}}$$ です。
data の設定
Stan の data ブロックに寄せて書いてみました。
PyMC の Data() コンテナを用いて、売り上げデータを変数名 $${Y}$$ に設定しています。
なお、観測データを Data で定義しなくてもモデル記述できます。
事前分布の設定
Stan の parameters ブロックに寄せて書いてみました。
$${\mathtt{変数名=pm.確率分布名(PyMC内の変数名, 確率分布のパラメータ値)}}$$
みたいな感じで書きます。
このモデルでは、平均パラメータ $${\mu}$$ に平均 $${0}$$、標準偏差 $${100000}$$ の正規分布 $${\mathtt{Normal}}$$、標準偏差パラメータ $${\sigma}$$ に標準偏差 $${100000}$$ の半正規分布 $${\mathtt{HalfNormal}}$$ を設定します。
尤度（観測値と連動）の設定
Stan の model ブロックに寄せて書いてみました。
事前分布と似た構文ですが、観測データと関連付けるパラメータ $${\mathtt{observed}}$$ が特徴的です。
特に断りしない場合、このシリーズ記事では変数名 $${\mathtt{obs}}$$（観測値：observed data）を使っていきます。
coords の設定
冒頭の coords はデータの座標（または次元・軸）の名称や値です。
このモデルでは $${\mathtt{id}}$$ に観測データのインデックス（0から始まる連番）を設定して、各変数の引数 $${\mathtt{dim}}$$ で用います。
なお、coords を設定しなくてもモデル記述できます。

🚀🚀🚀

モデルの内容を表示してみましょう。
「モデル名」を実行するとモデル内容（数式ライク）を見られます。

# モデルの表示
model

【実行結果】
変数名と確率分布の内容が表示されました。

続いてグラフィカルモデルの形式でモデルを描画します。
テキスト 4.16 節のモデル表記法と比べてみましょう。

この機能は、Graphviz ライブラリがインストールされている場合に利用できます。
Graphviz のインストールは、例えばこちらのサイトが参考になります。
ありがとうございます！

$${\mathtt{pm.model\_to\_graphviz(モデル名)}}$$ を指定して実行します。

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

【実行結果】
図が現れました！
グラフィカルモデルの方がモデルの直感的理解に役立ちそうです。

【図の見方】
角丸四角の枠の右下の $${\mathtt{id}}$$ は coords で指定した座標です。
変数 $${Y}$$ と $${\text{obs}}$$ はデータ数 $${100}$$ の $${\mathtt{id}}$$ 座標（軸と呼ぶほうがいいかも）で包まれています。
$${\text{obs}}$$ に向かってパラメータを示す変数 $${\text{mu},\ \text{sigma}}$$ がいます。この変数たちは座標を設定しておらず、それぞれデータ数は $${1}$$ です。

図形の形状や色の意味合いは以下の表を参考にして下さい。

🚀🚀🚀

MCMCサンプリングの実行

テキスト 4.12 節に相当します。
モデリングと同様に with 構文を使います。
sample 関数を使って MCMC サンプリングを行います。

%%time
# MCMCの実行 ※高速化したいときは nuts_sampler='nutpie'
with model:
    idata = pm.sample(draws=1000, tune=1000, chains=4, random_seed=1)

【実行結果】
自環境がいまいちなので処理時間に１分以上かかっていますが、ちゃんとした環境であればおそらく秒殺で処理が完了します。

【コード補足】
$${\texttt{pm.sample()}}$$ の引数はテキストのStan引数と合わせています。

draw：バーンイン期間後のMCMCサンプルの数
tune：バーンイン期間（Stan の warmup に相当）
chains：チェーン数（Stan の chains に相当）
random_seed：乱数の種（Stan の seed に相当）

この処理によって、MCMCサンプルは $${\mathtt{chains}:4 \times \mathtt{draw}:1000 = 4000}$$ 個作成されます。
MCMC サンプルは「パラメータの事後分布」から取得した乱数に相当します。

MCMC サンプルや MCMC にまつわる情報は $${\texttt{pm.sample()}}$$ を代入する変数に格納されます。
今回は変数 $${\mathtt{idata}}$$ に格納しています。
idata は「Inference Data」（推論データ）の略称です。
arviz ライブラリの InferenceData 形式で格納されます。

【出力結果の確認点】
Divergences（ダイバージェンス）が $${0}$$ になっていること
（MCMCサンプリングが不安定のとき、$${1}$$ 以上になります）

🚀🚀🚀

推定結果を確認

テキスト 4.13 節に相当します。
パラメータ推定値の統計量を arviz ライブラリの summary 関数で表示します。
$${\mathtt{idata}}$$ の中身を要約します。

# 結果の表示
az.summary(idata, hdi_prob=0.95)

【実行結果】
テキスト p.119 青領域の中段の表に似ています。
4000 個の MCMC サンプルの統計量が表示されています。

【表の見方（概要）】
パラメータ $${\mathtt{mu}}$$ の推定値の平均（mean）は $${102.166}$$、標準偏差(sd)は $${1.813}$$、95%HDI は $${[98.438, 105.483]}$$ です。
パラメータ $${\mathtt{sigma}}$$ の推定値の平均（mean）は $${18.200}$$、標準偏差(sd)は $${1.303}$$、95%HDI は $${[15.818, 20.798]}$$ です。

固いことを抜きにすると…
ざっくり売り上げの平均値は 102 (95%HDI [98, 105]) です。
ざっくり売り上げの標準偏差（ばらつき）は 18.2 (95%HDI [15.8, 20.8]) です。

🚀🚀🚀

【HDIとベイズ信用区間の違い】
テキストの 95% 区間は「95% ベイズ信用区間」です。
一方、この記事の区間は「95% HDI（最高密度区間）」です。。
PyMC・ArviZ 界隈は HDI を用いるようですので、ここで補足いたします。
下の図のイメージをご覧ください。

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

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

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

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

🚀🚀🚀

収束の確認

テキスト 4.14 節に相当します。
テキストは３つの方法で MCMC サンプル（事後分布）の収束状況をチェックしています。

最初の２つの指標はさきほどの要約表に掲載されています。

① 有効サンプル数
要約表の「ess_bulk」が事後分布の中央部分の有効サンプル数、「ess_tail」は事後分布の裾部分（極端な値）の有効サンプル数です。
テキスト掲載の目安「$${100}$$ くらい」を満たしています。

② $${\widehat{R}}$$（アールハット）
要約表の「r_hat」が $${\widehat{R}}$$ です。
MCMC サンプル（を取得した事後分布）が収束しているかどうかの判断に用いられる指標の１つであり、テキスト掲載の目安「$${1.1}$$ 未満」を満たしています。

③ トレースプロット
トレースプロットを描画して事後分布の収束有無を確認します。
arviz ライブラリの plot_trace 関数を利用します。

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

【実行結果】
右のゲジゲジチャートがトレースプロットです見ましょう。
４本の chains が色違いで描画されています。
テキストの「４本のチェーンがまじりあっていれば良いです」の状況になっている（ゲジゲジしている）ので満たしています。
なお、左の分布は MCMC サンプルのKDE プロットです。
４本のチェーンがほぼ同じ分布になっていることが分かります。

ちなみにパラメータの事後分布のKDEプロットは次のコードでも描画できます。
arviz ライブラリの plot_posterior 関数を利用します。

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

【実行結果】
推定値（数字の目盛り）が表示されていて、パラメータの事後分布の形状や値が見やすくなっています。

さきほどの要約表で得た売り上げの「ざっくり感」を事後分布プロットからも感じ取れますね！

固いことを抜きにすると…
ざっくり売り上げの平均値は 102 (95%HDI [98, 105]) です。
ざっくり売り上げの標準偏差（ばらつき）は 18.2 (95%HDI [15.8, 20.8]) です。

モデリング、MCMCサンプリング、推定結果確認、収束確認の一連の手続きを概観できました。
Python・PyMC コード作成に慣れるまでの間は、今回の基本コードをコピペ利用するのも「アリ」です。

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

シリーズの記事

次の記事

前の記事

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

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

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

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