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

書籍の著者　馬場真哉 先生

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

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

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

はじめに

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

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

準備

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

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

■ コード記述法
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'

■ CmdStanPyとArviZのバージョン確認
CmdStanPy はベイズ統計モデリングライブラリです。
ArviZ はベイズ統計モデリングなどの可視化ライブラリです。
記事で用いる両ライブラリのバージョン確認をします。

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

【実行結果】

第４章 Stanの基本

CmdStanPyのインストール

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

CmdStan、CmdStanPy、C++ツールチェーンの整合を保ってインストールする必要があるとのことです。
公式サイトで紹介されている方法「Conda、PyPI、GitHUB のいずれか」でインストールするのが良さそうです。

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

（参考）Python の Stan ライブラリの選択肢
Python の Stan ベースのベイズ統計モデリングライブラリには、CmdStanPy のほかに、PyStan があります。
もしかすると、PyStan の方がメジャーかもしれません。
ただし 最新の PyStan が動作する OS は Linux と macOS です。
本シリーズは Windows 上で（WSLを使わずに）Stan を動かしたくて、CmdStanPy を選択しました。

🔵🔵🔵

データの読み込み

テキスト 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 流の ベイズ統計モデルを記述します。
テキストにならって、Stan ファイルに Stan コードを書きます。

① Stan ファイルを格納するフォルダ
記事では Jupyter Notebook ファイル（.ipynb）と同じフォルダに「Stan」フォルダを作成しています。
Stan ファイルは この Stan フォルダ配下に格納します。

② Stan ファイル（Stan コード）
テキストの Stan ファイル名 および Stan コードを引用いたします。
文字コードは UTF-8 を使います。
📑ファイル名：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 コードの解説はぜひテキストでご確認下さい。

このモデルは、次の尤度の確率分布だけを明示的に書き、パラメータ $${\mu,\ \sigma}$$ の事前分布は省略しています。
（ベイズ統計はパラメータを「確率変数」として扱います！）

🔵🔵🔵

MCMCサンプリングの実行

テキスト 4.11 節と 4.12 節に相当します。

まず Stan に渡す サンプルサイズと売り上げデータを辞書にまとめます。

# p.118 サンプルサイズ
sample_size = len(file_beer_sales_1)
sample_size

【実行結果】

# p.118 リスト、ではなく、辞書にまとめる
data_dict = {
    'N': sample_size,                    # サンプルサイズ
    'sales': file_beer_sales_1['sales']  # 売り上げデータ
}
data_dict

【実行結果】
サンプルサンプル N は整数値、売り上げデータ sales は pandas シリーズ型です。

続いて「モデルのコンパイル」を実行します。
はじめの２行で Stan ファイルのパスを編集しています。
最後の１行で Stan ファイルのパスを受け取って、Stan モデルの実行ファイルを作成するよう指示しています。

%%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 フォルダに２つのファイルが作成されました。

• 2-4-1-calc-mean-variance.exe

• 2-4-1-calc-mean-variance.hpp

.exe ファイルは Stan モデルの実行ファイルです。
.hpp ファイルは実行ファイルを作成するための中間ファイルです。

いよいよ「MCMC サンプリング」の実行です。
モデルコンパイル時につくった「model」に対して sample メソッドを実行し、結果を変数 $${\texttt{fit}}$$ に格納しています。
data 引数には先ほど作った辞書 data_dict を与えています。
その他の引数はテキストに準拠しています。

%%time
# p.118 MCMCを実行する
fit = model.sample(
    data=data_dict,          # 対象データ
    seed=1,                  # 乱数の種
    chains=4,                # チェーン数
    iter_warmup=1000,        # バーンイン期間（１チェーンあたり）
    iter_sampling=1000,      # サンプリング数（１チェーンあたり）
    thin=1,                  # 間引き数（１なら間引き無し）
)

【実行結果】
無事に MCMC サンプリングを実行できたようです。

🔵🔵🔵

推定結果を確認

テキスト 4.13 節に相当します。
結果を格納した $${\mathtt{fit}}$$ に対して summary メソッドを実行し、パラメータ推定値の統計量を表示します。

# p.119 結果の表示 要約統計量
fit.summary(percentiles=(2.5, 50, 97.5)).round(2)

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

【表の見方（概要）】
パラメータ $${\mathtt{mu}}$$ の推定値の平均（Mean）は $${102.18}$$、標準偏差(StdDev)は $${1.82}$$、95%ベイズ信用区間 は $${[98.41, 105.20]}$$ です。
パラメータ $${\mathtt{sigma}}$$ の推定値の平均（Mean）は $${18.20}$$、標準偏差(StdDev)は $${1.29}$$、95%ベイズ信用区間 は $${[15.83, 20.84]}$$ です。

🔵🔵🔵

収束の確認

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

CmdStanPy には診断サマリーを表示する機能があります。
結果を格納した $${\mathtt{fit}}$$ に対して diagnose メソッドを実行し、診断サマリーを表示します。

# p.119 結果の表示 診断
print(fit.diagnose())

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

下の２項目が有効サンプル数（Effective sample size）、$${\widehat{R}}$$（R-hat values）です。

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

① 有効サンプル数
要約表の「N_Eff」が事後分布の有効サンプル数です。
テキスト掲載の目安「$${100}$$ くらい」を満たしています。

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

③ トレースプロット
トレースプロットを描画して事後分布の収束有無を確認します。
まず arvis の from_cmdstanpy 関数を利用して $${\mathtt{fit}}$$ を arviz の idata 形式に変換します。

# arvizのトレースプロット機能を使うために、arvizのidata形式に変換する
idata = az.from_cmdstanpy(posterior=fit, log_likelihood='lp__')
idata

【実行結果】
PyMC 利用者にはおなじみの idata になっています！

トレースプロットを描画します。
arviz ライブラリの plot_trace 関数を利用します。

# p.120 図2.4.1 トレースプロットの描画　※バーンイン期間なし
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);

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

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

（注）
CmdStamPy の summary() にはベイズ信用区間が表示されますが、arviz の summary() や 事後分布プロットなどには HDI（最高密度区間）が表示されます。

モデリング、MCMCサンプリング、推定結果確認、収束確認の一連の手続きを概観できました。
CmdStanPy の使い方はテキストの Stan と少々異なるようです。
Python・CmdStanPy コード作成に慣れるまでの間は、今回の基本コードをコピペ利用するのも「アリ」です。

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

シリーズの記事

次の記事

前の記事

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

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