見出し画像

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

ネイピア DS

書籍の著者　馬場真哉 先生

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

テキストのこの章は Stan コードを詳しく解説しています。
Stan コードの解説はテキストをお読みいただくとして…
この記事では、Python・CmdStanPy で「平均値の差」のベイズ流推論を実践します！

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


はじめに

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

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

準備

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

6.8 平均値の差の評価と generated quantities ブロック

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

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

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

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

# インポート

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

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

以下のコードは基本的にテキストの 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}$$ 検定へ！
ではなくて、ベイズ流の平均値差の推論へ進みます。

🔵🔵🔵

モデリング

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

$$
\begin{align*}
\text{sales}_A &\sim \text{Normal}\ (\mu_A,\ \sigma_A^2) \\
\text{sales}_B &\sim \text{Normal}\ (\mu_B,\ \sigma_B^2) \\
\text{diff} &= \mu_B - \mu_A
\end{align*}
$$

$${\text{sales}_A,\ \text{sales}_B}$$ はそれぞれビール A、ビール B の売り上げ量に関する尤度であり、正規分布に従うと仮定しています。
こちらは model ブロックで記述します。

２つの正規分布の平均パラメータ $${\mu_A,\ \mu_B}$$ および標準偏差パラメータ $${\sigma_A,\ \sigma_B}$$ には事前分布を明示的に設定しません。
こちらは parameters ブロックに変数と型を記述します。

最後の変数 $${\text{diff}}$$ はビール B の売り上げ量の期待値とビール A の売り上げ量の期待値の差であり、「平均値差」を表します。
こちらは generated quantities ブロックで記述します。

② Stan のモデル設定
テキストと同様に銘柄別にデータを分割し、Stan に渡すデータセットを辞書にまとめます。

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

# ビールの種類別にデータを分ける
sales_a = file_beer_sales_ab[file_beer_sales_ab['beer_name']=='A']['sales']
sales_b = file_beer_sales_ab[file_beer_sales_ab['beer_name']=='B']['sales']

# 標本サイズ（A,B同じ）
N = len(sales_a)

# 辞書にまとめる
data_dict_ab = dict(N=N, sales_a=sales_a, sales_b=sales_b)

【実行結果】なし

Stan ファイル（Stan コード）を作成します。
テキストの Stan ファイル名 および Stan コードを引用いたします。
📑ファイル名：2-6-5-difference-mean.stan

data {
    int N;                  // サンプルサイズ
    vector[N] sales_a;      // ビールAの売上データ
    vector[N] sales_b;      // ビールBの売上データ
}

parameters {
    real mu_a;              // ビールAの平均
    real<lower=0> sigma_a;  // ビールAの標準偏差
    real mu_b;              // ビールBの平均
    real<lower=0> sigma_b;  // ビールBの標準偏差
}

model {
    // 平均mu, 標準偏差sigma の正規分布に従ってデータが得られたと仮定
    sales_a ~ normal(mu_a, sigma_a);
    sales_b ~ normal(mu_b, sigma_b);
}

generated quantities {
    real diff;              // ビールAとBの売上平均の差
    diff = mu_b - mu_a;
}

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

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

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

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

【実行結果】（右側のファイルパスは記載省略）

🔵🔵🔵

MCMC の実行

MCMCを実行しましょう。

%%time
# p.149 乱数の生成（MCMCの実行）
fit = model.sample(
    data=data_dict_ab,       # 対象データ
    seed=1,                  # 乱数の種
    sig_figs=18,             # 出力CSV等に適用する数値精度
)

【実行結果】

収束の確認をします。
診断メソッド diagnose を利用します。

# 事後分布の診断
print(fit.diagnose())

【実行結果】（１行目のファイルパスは記載省略）
問題は検出されませんでした（no problems detected.）。

  1. treedepth（木の深さ）：
    satisfactory（満たしている、OK）

  2. divergences（ダイバージェンス）：
    No Divergent（ダイバージェンスなし、OK）

  3. E-BFMI（期待される情報不足のベイズ的割合）：
    satisfactory（満たしている、OK）

  4. Effective sample size（有効サンプル数）：
    satisfactory（満たしている、OK）

  5. R-hat values（$${\widehat{R}}$$ 値）：
    satisfactory（満たしている、OK）

Stan の MCMC サンプルを arviz の idata 形式に変換して、具体的な値を見ていきましょう。

# arvizのidataに変換
idata = az.from_cmdstanpy(posterior=fit, log_likelihood='lp__')

【実行結果】なし

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, compact=False, var_names=var_names,
              backend_kwargs={'tight_layout': True});

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

以上のチェックに基づいて、収束している！と噛み締めましょう。

🔵🔵🔵

平均値の差の評価

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

平均値差 diff は ビール B の売り上げ量 からビール A の売り上げ量を引いたものです。
ベイズ流の平均値差 diff は 事後平均 66.7（標準誤差 3.5）、95% HDI [ 60.0, 73.5 ]（０を含まない）です。
ビール 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より大きい確率は 100% です。

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

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

シリーズの記事

次の記事

前の記事

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

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

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

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

この記事が参加している募集

#

最近の学び

31万件