「ベイズ統計モデリングによるデータ分析入門」をPythonとPyMCで写経 ～ Vol.15 ロジスティック回帰モデル
書籍の著者 馬場真哉 先生
この記事は、書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」第３部第９章「ロジスティック回帰モデル」の Python 写経活動記録です。
書籍の第３部は「一般化線形モデル」のベイズ統計モデリングです。
今回は「ロジスティック回帰モデル」を３つのツール・手法で取り組みます。
GLMによるロジスティック回帰モデル
Bambiのベイズロジスティック回帰モデル
PyMCのベイズロジスティック回帰モデル
では書籍を開いてベイズ統計モデリングの旅に出かけましょう🚀
はじめに
このブログシリーズは書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」（講談社、「テキスト」と呼びます）の Python 写経です。
テキストの紹介と引用表記はリンク先の記事に掲載しています。
準備
■ 記事の範囲
この記事はテキスト第３部第９章の以下の節を取り扱います。
9.4 データの読み込みと可視化
9.5 brmsによるロジスティック回帰モデルの推定
9.6 推定されたモデルの解釈
9.7 回帰曲線の図示
9.8 補足：ロジスティック回帰モデルのためのStanファイルの実装
ベイズ統計モデリングの理論面が気になった場合には、ぜひテキストの第１部、第３部第１章をご覧いただき、本記事との繋がりをご確認下さいませ。
■ コード記述法
Jupyter Notebook 形式でコードを記述します。
■ 利用データ
テキスト・サポートサイトのデータファイルを引用しています。
Jupyter Notebook ファイルと同一フォルダ内の「data」フォルダにデータファイルを格納しています。
■ ライブラリのインポート
この記事で用いるライブラリをインポートします。
# インポート
# 数値計算
import numpy as np
import pandas as pd
# 統計処理
from scipy.special import expit # ロジスティック関数
# ベイズ統計モデリング
import pymc as pm # pymc
import bambi as bmb # bambi
import arviz as az # 分析・可視化
# 統計モデリング
import statsmodels.api as sm
import statsmodels.formula.api as smf
# デザイン行列
from patsy import dmatrices, dmatrix
# 可視化
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_theme() # ggplot風のスタイル
plt.rcParams['font.family'] = 'Meiryo'第９章 ロジスティック回帰モデル
ロジスティック回帰モデルとは？
ロジスティック回帰モデルは一般化線形モデルの文脈において、
線形予測子：複数の説明変数を用いる
リンク関数：ロジット関数
確率分布：二項分布
で構成される統計モデルです。
目的変数が０以上の上限のある整数の場合に用いられます。
二項ロジスティック回帰モデルともよばれます。
今回のロジスティック回帰モデルは、目的変数が二値の分類問題に適用されるモデルでは無いことにご留意ください。
🚀🚀🚀
データの読み込みと外観の確認
テキスト 9.4 節に相当します。
テキストの仮想の植物の種子発芽数データを引用いたします。
csv ファイルを pandas のデータフレーム形式で変数 germination_dat に読み込みます。
# p.221 分析対象データの読み込み
# ファイルの読み込み
germination_dat = pd.read_csv('./data/3-9-1-germination.csv')
# 結果の表示
# germination: 発芽数、size: 植木鉢にまいた種子の数, solar: 日照の有無,
# nutrition: 栄養素の量
print('germination_dat.shape: ', germination_dat.shape)
germination_dat.head(3)【実行結果】
標本サイズ 100のデータです。
germination：種子の発芽数（成功数）
size：植木鉢に蒔いた種子の数（試行回数、すべて 10 ）
solar：日照の有無（あり：sunshine、なし：shade）
nutrition：栄養素の量
日照の有無・栄養素量と発芽数の関係をロジスティック回帰モデルで分析します。
日照の有無ごとのデータ件数をカウントします。
pandas の value_counts メソッドを利用します。
# 日照の有無ごとのデータ件数
germination_dat['solar'].value_counts().to_frame()【実行結果】
日照なし、ありで、それぞれ 50 件のデータがあります。
データの要約統計量を確認します。
まずは全体（量的変数のみ）です。
# データの要約統計量
germination_dat.describe().T.round(2)【実行結果】
発芽数は平均 2.8、最小値 0、最大値 10 です。
栄養素量は平均 5.5、最小値 1、最大値 10 です。
次は日照の有無別の発芽数の要約統計量です。
# 日照有無別の発芽数の要約統計量
germination_dat.groupby(['solar'])['germination'].describe().round(2)【実行結果】
日照あり（sunshine）の発芽数の方が大きいです！
次は日照の有無別の栄養素量の要約統計量も念のため見ておきましょう。
# 日照の有無別の栄養素量の要約統計量
germination_dat.groupby(['solar'])['nutrition'].describe().round(2)【実行結果】
違いは見られません。
データを可視化しましょう。
最初にジョイントプロットで、要約統計量の感覚値を可視化してみます。
seaborn の jointplot を利用します。
# ジョイントプロットの描画
colors = {'shade': 'tomato', 'sunshine': 'tab:blue'}
sns.jointplot(
data=germination_dat, x='nutrition', y='germination',
hue='solar', palette=colors
);【実行結果】
散布図では栄養素量の増加と発芽数の増加＆ばらつき増加の関係が見られ、しかも、日照ありの方が発芽数が大きい傾向が見られます。
栄養素量のKDEプロット（上）では、日照の有無で同じ分布です。
発芽数のKDEプロット（右）では、日照なしのばらつきが小さく、日照ありのばらつきが大きいことが分かります。
続いて、テキスト 図 3.9.1 に相当する散布図です。
seaborn の scatterplot を利用します。
# p.222 図3.9.1 種子の発芽数と日照・栄養素の散布図
# 描画領域の設定
plt.figure(figsize=(8, 4))
# 散布図
colors = {'shade': 'tomato', 'sunshine': 'tab:blue'}
sns.scatterplot(data=germination_dat, x='nutrition', y='germination',
hue='solar', palette=colors)
# 修飾
plt.title('種子の発芽数と日照の有無・栄養素の量の関係', loc='left');【実行結果】
ジョイントプロットで見てしまいましたが…栄養素量の増加と発芽数の増加＆ばらつきの増加の関係が見られます。
この散布図に天気ごとの回帰曲線を重ねてみましょう。
seaborn の lmplot を利用します。
lowess 引数を ON にして局所的な傾きの違いを表現します。
# 種子の発芽数と日照・栄養素の散布図（回帰曲線付き）
# 回帰曲線付き散布図
colors = {'shade': 'tomato', 'sunshine': 'tab:blue'}
sns.lmplot(data=germination_dat, x='nutrition', y='germination', lowess=True,
hue='solar', palette=colors, height=4, aspect=1.8)
# 修飾
plt.title('種子の発芽数と日照の有無・栄養素の量の関係（回帰曲線付き）', loc='left');【実行結果】
日照有無にかかわらず、栄養素の増加と種子数の増加の関係が見られます。
また日照ありの種子数は相対的に大きいです。
🚀🚀🚀
GLMによるロジスティック回帰モデル
① モデルの概要
まずGLMでロジスティック回帰モデルを確認しておき、あとでベイズ統計モデルの結果と比べてみましょう。
Python の統計ライブラリ statsmodels を利用して、次のモデルを実装します！
$$
\begin{align*}
\eta_i &= \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} \\
\text{logit}(p_i) &= \eta_i \\
p_i &= \text{logistic}(\eta_i) \\
y_i &\sim \text{Binomial}(10, p_i)
\end{align*}
$$
$${\eta}$$ は線形予測子、$${\beta}$$ は係数、$${x_{i1}}$$ は日照ありダミー、$${x_{i2}}$$ は栄養素量です。
$${\text{logit}(p)}$$ はロジットリンク関数、$${\text{logistic}(\eta)}$$ は逆リンク関数（ロジスティック関数、シグモイド関数）です。
$${p}$$ は二項分布の成功率パラメータ、$${y}$$ は目的変数となる発芽数です。
このロジスティック回帰モデルの線形予測子を statsmodels の一般化線形モデル glm にあたえる formula 構文に変換します。
$${\sim}$$ を挟んで、左辺は目的変数、右辺は説明変数です。
formula の目的変数に「成功数 + 失敗数」を与えるのがポイントです。
$$
\texttt{germination + I(size - germination)} \sim \texttt{solar} + \texttt{nutrition}
$$
② ロジスティック回帰の実行
ではロジスティック回帰を実行します！
また確率分布 family を sm.families.Binomial() で指定します。
二項分布のデフォルトリンク関数はロジットリンク関数です（指定不要）。
# ロジスティック回帰モデル by statsmodels
# formulaの設定：目的変数は「成功数 + 失敗数」
formula = 'germination + I(size - germination) ~ solar + nutrition'
# 確率分布：二項分布とリンク関数：デフォルト（ロジット関数）の設定
family=sm.families.Binomial()
# ロジスティック回帰の実行
res_sm = smf.glm(formula=formula, data=germination_dat, family=family).fit()
res_sm.summary()【実行結果】
こちらは GLM のサマリーです。
③ 推定値の確認
係数の推定値にフォーカスしてみます。
# 要約表から係数の推定値の部分を取り出し
res_sm.summary().tables[1]【実行結果】
【ちょっと深堀り】
ロジスティック回帰の線形予測子の係数の解釈は少々ややこしいです。
生成 AI にやさしい解説をお願いしました🤖
キーワード「オッズ」は $${\cfrac{p}{1-p}}$$ です。
テキストでは「オッズとは失敗するよりも何倍成功しやすいかを表した指標」と説明しています。
では Gemini の解説をどうぞ！
【解説】ロジスティック回帰の係数：正体は「オッズの掛け算」
ロジスティック回帰の結果（係数 $${\beta}$$）を見て、「1増えると確率が〇％上がる」と解釈していませんか？実はその解釈、厳密には少し違います。
ロジスティック回帰の係数を正しく読み解くポイントは、「足し算ではなく掛け算」、そして「確率ではなくオッズ」で考えることです。
1. 係数を「指数（exp）」に変換する
統計ソフトが出力する係数 $${\beta}$$ は、そのままでは使いにくい数値です。まずこれを $${e^\beta}$$（指数） に変換しましょう。
この変換後の数値が、「オッズ比（オッズが何倍になるか）」を表します。
2. オッズ比（$${e^\beta}$$）の読み方
$${e^\beta = 2.0}$$ の場合：
変数が 1 増えると、起こるオッズが 2倍 になる。$${e^\beta = 0.5}$$ の場合：
変数が 1 増えると、起こるオッズが 0.5倍（半分） になる。
3. 「確率」ではなく「オッズ」である理由
「オッズが2倍」といっても、「確率が2倍」になるとは限りません。
もともとの確率が 1% なら、オッズ2倍で確率は 約 2%（ほぼ2倍）。
もともとの確率が 50% なら、オッズ2倍で確率は 約 67%（1.3倍程度）。
確率は $${100\%}$$ を超えられないため、係数が増えるほど「確率の伸び」は緩やかになるという性質があります（これがロジスティック回帰特有の $${S}$$ 字カーブの正体です）。
4. まとめ
ロジスティック回帰の係数を見たら、まずは np.exp(beta) を計算！
出てきた数字は「発生しやすさ（オッズ）が何倍ブーストされるか」の指標だと考えましょう。
「オッズが $${\exp(係数の推定値)}$$ 倍」と覚えておきましょう。
日照ダミーの係数の場合、日照なしと比べて日照ありはオッズが $${\exp(4.0113) \approx 55.2}$$ 倍（大きすぎ !?）。
栄養素量の場合、栄養素の１単位の増加でオッズが $${\exp(0.7142) \approx 2.0}$$ 倍。
④ 日照の有無別・栄養素量別の平均発芽数の可視化
最後に日照の有無別・栄養素量別の平均発芽数の95%信頼区間を可視化しましょう。
最初に平均釣獲尾数の予測値を算出します。
ロジスティック回帰の結果 res_sm に対して get_prediction メソッドを適用するなどして、予測値を得ます。
## 平均発芽数の予測値の算出
# 予測用データの作成
n_samples = 100
test_data = pd.DataFrame({
'solar': np.repeat(['shade', 'sunshine'], n_samples), # 日照の有無
'nutrition': np.tile(np.linspace(0, 10, n_samples), 2) # 栄養素の量（連続値）
})
# 予測値の取得
prediction = res_sm.get_prediction(test_data)
pred_df = prediction.summary_frame(alpha=0.05) # alpha=0.05で95%信頼区間
# 予測値データフレームの作成
pred_df = pd.concat([test_data, pred_df], axis=1)
pred_df【実行結果】
それでは描画します。
## 平均発芽数の95%信頼区間の描画
# 色の設定
colors = {'shade': 'tomato', 'sunshine': 'tab:blue'}
# 描画領域の設定
plt.figure(figsize=(8, 4))
# 観測値の散布図の描画
sns.scatterplot(data=germination_dat, x='nutrition', y='germination',
hue='solar', palette=colors)
# 日照の有無ごとに平均発芽数の予測値と95%信頼区間の描画を繰り返し処理
for sol, color in colors.items():
# 描画データの取得
temp, mean_germ, lower, upper = (
pred_df[pred_df['solar']==sol].iloc[:, [1, 2, 4, 5]].to_numpy().T
)
# 平均発芽数の予測値（点推定）の描画
plt.plot(temp, mean_germ*10, color=color)
# 平均発芽数の95%信頼区間の塗りつぶし描画
plt.fill_between(temp, lower*10, upper*10, color=color, alpha=0.2)
# 修飾
plt.title('平均発芽数の95%信頼区間')
plt.xlabel('栄養素の量', fontsize=12)
plt.ylabel('発芽数', fontsize=12);【実行結果】
栄養素の量が大きくなるにつれて発芽数は S 字形で大きくなることが分かります。
日照ありの方が明らかに発芽数が大きいです。
ベイズ流のロジスティック回帰モデルへ進みます。
🚀🚀🚀
ベイズモデリング by Bambi
テキスト 9.5、9.7 節に相当します。
brms の代わりに Bambi を利用します。
① モデルの概要
確率分布を二項分布、リンク関数をロジットリンク関数（二項分布のデフォルト）、線形予測子を次の formula で実装します。
$$
\texttt{p(germination, size)} \sim \texttt{solar} + \texttt{nutrition}
$$
目的変数には $${\texttt{p(成功数, 試行回数)}}$$ を設定します。
公式サイトの実装例によると「成功数を試行回数で割った結果得られる比率をモデル化する」とのこと。
テキストの無情報事前分布に合わせるべく、Bambi に以下の事前分布情報を与えます。
$$
\begin{align*}
intercept &\sim \text{Normal}\ (0, (1e5)^2) \\
solar &\sim \text{Normal}\ (0, (1e5)^2) \\
nutrition &\sim \text{Normal}\ (0, (1e5)^2) \\
\end{align*}
$$
② モデル定義
上式のモデルを Bambi で記述します。
bmb.prior() で事前分布を定義します。複数ある場合は辞書でまとめます。
bmb.Model() において、prior 引数で定義した事前分布を与えます。
# モデリング
# 無情報事前分布（想定）の設定
uninformed_prior = bmb.Prior('Normal', mu=0, sigma=1e5)
# 事前分布を辞書にとりまとめ
priors = {'Intercept': uninformed_prior, 'solar': uninformed_prior,
'nutrition': uninformed_prior}
# モデルの定義
# formulaの左辺は目的変数と試行回数の比率を設定
# 書き方は３通り：proportion(y, n), prop(y, n), p(y, n)
model_bmb = bmb.Model(
formula='p(germination, size) ~ solar + nutrition', # フォーミュラ式
data=germination_dat, # データ
family='binomial', # 確率分布：二項分布
# link='logit' # リンク関数：ロジット
priors=priors, # 事前分布（辞書）
)【実行結果】なし
モデルの内容を表示します。
# モデルの表示
model_bmb【実行結果】
確率分布 Family に 二項分布 binomial、リンク関数 Link にロジットリンク関数 p= logit が設定されました。
モデルをグラフィカルモデルで描画します。
# モデルの可視化
model_bmb.build()
model_bmb.graph()【実行結果】
solar 係数の次元は１です（solar_dim (1) の 1 ）。
ダミー変数化するときに要素の１つを除外しています。
③ MCMC の実行
MCMCを実行しましょう。
NUTS サンプラーに nutpie を利用します。
%%time
# MCMCの実行
idata_bmb = model_bmb.fit(
draws=1000, tune=1000, chains=4, random_seed=123, nuts_sampler='nutpie'
)【実行結果】
Divergences（ダイバージェンス）は０件です。
④ 収束確認
収束の確認をします。
MCMC サンプルの要約表を表示します。
# p.223 要約統計量の表示
var_names = ['Intercept', 'solar', 'nutrition']
az.summary(idata_bmb, var_names=var_names, hdi_prob=0.95)【実行結果】
$${\widehat{R}}$$（r_hat）、有効サンプル数（ess_bulk、ess_tail）に問題はなさそうです。
ちなみに $${\exp(4.059) \approx 57.92}$$ 倍、$${\exp(0.722) \approx 2.06}$$ 倍です。
GLMの推定値の表を並べましょう。
（参考：statsmodels の推定値）
手法間で係数推定値は近い値になっています。
トレースプロットを描画します。
# トレースプロットの描画
az.plot_trace(idata_bmb, compact=False, var_names=var_names,
backend_kwargs={'tight_layout': True});【実行結果】
左側のチャートの４本の Chain はほぼ重なっており、かつ１峰です。
右側のチャートがゲジゲジしています。
以上のチェックに基づいて、収束していると考えましょう。
⑤ 推定されたモデルの解釈
テキスト 9.6 節に相当します。
テキストに準拠して進めます。つまり、223 ページ最終行からの「回帰係数とオッズ比の関係を」Python「で確認します」！
a. 説明変数を作る
発芽率（成功率） $${p}$$ の予測値を算出するために、予測用の説明変数を作成します。
# p.224 説明変数を作る
# solarとnutritionのデータフレームを作成
newdata_1 = pd.DataFrame({
'solar': ['shade', 'sunshine', 'sunshine'],
'nutrition': [2, 2, 3],
})
print('【説明変数】')
display(newdata_1)【実行結果】
３行のデータを作成しました。
b. 発芽率 $${p}$$ の予測値の計算
model_bmb に対して predict メソッドを適用して $${p}$$ の予測値を計算します。
# p.224 発芽率の予測値（点推定値）の計算
# 利用するMCMCサンプルの番号
num_samples = 2
# 成功確率の予測値の計算（ロジスティック関数適用済み） predictメソッド利用
fit = (
model_bmb.predict(idata_bmb, data=newdata_1, inplace=False)
.posterior['p']
.stack(sample=('chain', 'draw'))
.to_numpy()[:, num_samples]
)
print('【発芽率（成功確率）の予測値】')
print(fit)【実行結果】
３つのデータの発芽率です。
c. オッズの計算
発芽率にロジスティック関数（シグモイド関数）を適用してオッズを計算します。
# p.224 オッズの計算
oddses = fit / (1 - fit)
print('【オッズ】')
print(oddses)【実行結果】
３つのデータのオッズです。
d. モデルの係数を取得
MCMCサンプルから係数の推定値を取得します。
# p.225 モデルの係数の取得
coef = (
idata_bmb.posterior[['Intercept', 'solar', 'nutrition']]
.sel({'chain': 0, 'draw': num_samples})
.to_array()
.to_numpy()
.squeeze()
)
print('【モデルの係数】')
print(coef)【実行結果】
左から順に切片、solar、nutrition の係数推定値です。
e. オッズ比と exp(係数) の一致を確認
データ No.0 から No.1 へ変化（solar が日照なしから日照ありに変化）するときのオッズ比が、solar の係数推定値の指数 exp と一致することを確認します。
# p.225
# solarがshadeからsunshineに変わったときのオッズ比は
# solarの係数にexpを取ったものと等しくなる
print('【solar：shade → sunshineに変わったときのオッズ比】')
print(f'odds ratio : {oddses[1] / oddses[0]:.10f}')
print(f'exp(coef) : {np.exp(coef[1]):.10f}')【実行結果】
続いて、データ No.1 から No.2 へ変化（nutrition が 2 から 3 に変化）するときのオッズ比が、nutrition の係数推定値の指数 exp と一致することを確認します。
# p.225
# nutritionが2から3に変わったときのオッズ比は
# nutritionの係数にexpを取ったものと等しくなる
print('【nutrition：2 → 3に変わったときのオッズ比】')
print(f'odds ratio : {oddses[2] / oddses[1]:.10f}')
print(f'exp(coef) : {np.exp(coef[2]):.10f}')【実行結果】
以上でオッズ比の計算例を終了します。
⑥ 日照の有無ごとの平均発芽数データの作成
Bambi モデル model_bmb に対して predict メソッドを適用することで平均発芽数のMCMCサンプルを取得します。
# 日照の有無ごとの平均発芽数のMCMCサンプルの算出
## 設定
# 予測データの個数
n_pred = 100
# 試行回数
n_trials = 10
# 栄養素量の値
x_val = np.linspace(germination_dat['nutrition'].min(),
germination_dat['nutrition'].max(), n_pred)
# 日照の有無の要素
solars = ['shade', 'sunshine']
## 予測の実行
results_n_pred = []
# 日照の有無ごとに予測を繰り返し処理
for solar in solars:
# 予測データの作成
new_data = pd.DataFrame({
'size': [n_trials]*n_pred, 'solar': [solar]*n_pred, 'nutrition': x_val}
)
# 平均確率の事後予測サンプルの取得
result_p_pred = model_bmb.predict(
idata=idata_bmb, data=new_data, kind='response_params', inplace=False
)
# 平均発芽数に変換して結果をリストに格納
results_n_pred.append(az.extract(result_p_pred.posterior).p.data.T * n_trials)【実行結果】なし
⑦ 日照の有無ごとの平均発芽数の可視化
テキストの図 3.9.2 に相当します。
平均発芽数の平均値と 95%HDI 区間を描画します。
# p.226 図3.9.2 ロジスティック回帰モデルの回帰曲線：HDI区間付き
## 平均発芽数データと描画色を日照の有無ごとの辞書に格納
germins = {'shade': results_n_pred[0], 'sunshine': results_n_pred[1]}
colors = {'shade': 'tomato', 'sunshine': 'tab:blue'}
## 描画
# 描画領域の設定
plt.figure(figsize=(8, 4))
# 観測値の散布図の描画
sns.scatterplot(data=germination_dat, x='nutrition', y='germination',
hue='solar', palette=colors)
# 日照の有無ごとに平均発芽数と95%HDI区間の描画を繰り返し処理
for solar in germins.keys():
# 平均発芽数の平均値の描画
plt.plot(x_val, germins[solar].mean(axis=0), color=colors[solar])
# 平均発芽数の95%HDI区間の塗りつぶし描画
az.plot_hdi(x_val, germins[solar].reshape(4, 1000, -1), hdi_prob=0.95,
color=colors[solar], fill_kwargs={'alpha': 0.2})
# 修飾
plt.title('平均発芽数の95%HDI区間')
plt.xlabel('栄養素の量', fontsize=12)
plt.ylabel('発芽数', fontsize=12);【実行結果】
栄養素の量が大きくなるにつれて発芽数は S 字形で大きくなることが分かります。
日照ありの方が明らかに発芽数が大きいです。
statsmodels の 95% 信頼区間とも似ている感じです。
（参考：statsmodels の95%信頼区間）
🚀🚀🚀
ベイズモデリング by PyMC
テキスト 9.7、9.8 節に相当します。
デザイン行列を作成して PyMC でモデリングします。
① モデルの概要
次のモデルを実装します。
$$
\begin{align*}
\bm Y &\sim \text{Binomial} (10, \bm p) \\
\text{logit}(\bm p) &= \bm{X \beta} \\
\bm \beta &\sim \text{Normal} (0,\ (1e5)^2,\ \text{dims}=3) \\
\end{align*}
$$
② データセットの作成
patsy ライブラリを用いてデザイン行列等を作成します。
# デザイン行列の作成
# formula構文の設定
formula_binom = 'germination ~ solar + nutrition'
# 目的変数Y, デザイン行列（説明変数）Xの作成
Y_dm, X_dm = dmatrices(formula_binom, germination_dat, return_type='dataframe')
# デザイン行列の先頭5行の表示
X_dm.head()【実行結果】
デザイン行列は定数項、日照有無のダミー変数（日照あり）、栄養素の量で構成されます。
# 目的変数の先頭5行の表示
Y_dm.head()【実行結果】
目的変数も作ってくれました。
③ モデル定義
上式のモデルを PyMC で記述します。
# モデリング
# coordsの設定
coords = {'id': germination_dat.index.values, # 観測データの識別子
'coefs': ['intercept', 'sunshine', 'nutrition']} # 係数ベクトルbの識別子
# モデルの定義
with pm.Model(coords=coords) as model_pm:
## dataの設定
# 目的変数: 発芽数データ
Y = pm.Data('Y', value=Y_dm.values.flatten(), dims='id')
# 試行回数
n = pm.Data('n', value=germination_dat['size'].values, dims='id')
# 説明変数: デザイン行列
X = pm.Data('X', value=X_dm.values, dims=('id', 'coefs'))
## 事前分布: 無情報事前分布的な分布
# 係数ベクトル
b = pm.Normal('b', mu=0, sigma=1e5, dims='coefs')
## 尤度関数: 二項分布を仮定
p = pm.Deterministic('p', pm.math.sigmoid(X @ b), dims='id')
obs = pm.Binomial('obs', n=n, p=p, observed=Y, dims='id')
## （参考）尤度関数: logit_pを使う方法があります
# logit_p = pm.Deterministic('logit_p', X @ b, dims='id')
# obs = pm.Binomial('obs', n=n, logit_p=logit_p, observed=Y, dims='id')【実行結果】なし
モデルを数式ライクに表示します。
# モデルの表示
model_pm【実行結果】
モデルをグラフィカルモデルで描画します。
# モデルの可視化
pm.model_to_graphviz(model_pm)【実行結果】
④ MCMC の実行
MCMCを実行しましょう。
NUTS サンプラーに nutpie を利用します。
%%time
# MCMCの実行
with model_pm:
idata_pm = pm.sample(draws=1000, tune=1000, chains=4, random_seed=1,
nuts_sampler='nutpie')【実行結果】
Divergences（ダイバージェンス）は０件です。
⑤ 収束確認
収束の確認をします。
MCMC サンプルの要約表を表示します。
# 要約統計量の表示
var_names = ['b']
az.summary(idata_pm, var_names=var_names, hdi_prob=0.95)【実行結果】
$${\widehat{R}}$$（r_hat）、有効サンプル数（ess_bulk、ess_tail）に問題はなさそうです。
ちなみに $${\exp(4.042) \approx 56.94}$$ 倍、$${\exp(0.719) \approx 2.05}$$ 倍です。
（参考：Bambiモデルの要約統計量）
パラメータの事後分布はBambiと似た値になっています。
ただ…有効サンプル数（ess_bulk、ess_tail）がBambiよりも小さくなっている点は気になりますね…
トレースプロットを描画します。
# トレースプロットの描画
az.plot_trace(idata_pm, compact=False, var_names=var_names,
backend_kwargs={'tight_layout': True});【実行結果】
左側のチャートの４本の Chain はほぼ重なっており、かつ１峰です。
右側のチャートがゲジゲジしています。
以上のチェックに基づいて、収束していると考えましょう。
⑥ 日照の有無ごとの平均発芽数の可視化
テキストの図 3.9.3 に相当します。
平均発芽数の平均値と 95%HDI 区間を描画します。
パラメータのMCMCサンプルを使って、PyMCの外で平均発芽数の予測値を計算します。
# p.226 図3.9.2 ロジスティック回帰モデルの回帰曲線：HDI区間付き
## 平均発芽数データの作成
# MCMCサンプルからβ0, β1, β2を取り出し
intercept, w_sun, nutri = az.extract(idata_pm.posterior).b.values
# x軸の値
x_val = np.linspace(X_dm['nutrition'].min(), X_dm['nutrition'].max(), 1001)
# MCMCサンプルから日照の有無ごとの平均発芽数を算出
num_trials = 10
shade_nums = num_trials * expit(intercept + w_sun*0 + np.outer(x_val, nutri)).T
suns_nums = num_trials * expit(intercept + w_sun*1 + np.outer(x_val, nutri)).T
# 平均発芽数データと描画色を日照の有無ごとの辞書に格納
germins = {'shade': shade_nums, 'sunshine': suns_nums}
colors = {'shade': 'tomato', 'sunshine': 'tab:blue'}
## 描画
# 描画領域の設定
plt.figure(figsize=(8, 4))
# 観測値の散布図の描画
sns.scatterplot(data=germination_dat, x='nutrition', y='germination',
hue='solar', palette=colors)
# 日照の有無ごとに平均発芽数の平均値と95%HDI区間の描画を繰り返し処理
for solar in germins.keys():
# 平均発芽数の平均値の描画
plt.plot(x_val, germins[solar].mean(axis=0), color=colors[solar])
# 平均発芽数の95%HDI区間の塗りつぶし描画
az.plot_hdi(x_val, germins[solar].reshape(4, 1000, -1), hdi_prob=0.95,
color=colors[solar], fill_kwargs={'alpha': 0.2})
# 修飾
plt.title('平均発芽数の95%HDI区間')
plt.xlabel('栄養素の量', fontsize=12)
plt.ylabel('発芽数', fontsize=12);【実行結果】
Bambi モデルの 95% HDI区間とほぼほぼ似た形状になりました。
（参考：Bambiモデルの 95% HDI区間）
🚀🚀🚀
今回の記事は以上です。
楽しかったですね！
シリーズの記事
次の記事
前の記事
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の教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
いいなと思ったら応援しよう！
応援ありがとうございます。これからもがんばって記事を作成します！