「ベイズ統計モデリングによるデータ分析入門」をPythonとBambiで写経 ～ Vol.11 Bambiの使い方
書籍の著者 馬場真哉 先生
この記事は、書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」第３部第５章「brmsの使い方」の Python 写経活動記録です。
書籍の第３部は「一般化線形モデル」のベイズ統計モデリングです。
引き続き「単回帰モデル」を堪能します。
テキストの「R 言語と Stan ⇒ brms」の関係のように、Python と PyMC で簡単にベイズ統計モデリングを行えるライブラリは「Bambi」です。
今回記事は Bambi に取り組みます🦌
では書籍を開いてベイズ統計モデリングの旅に出かけましょう🚀
はじめに
このブログシリーズは書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」（講談社、「テキスト」と呼びます）の Python 写経です。
テキストの紹介と引用表記はリンク先の記事に掲載しています。
準備
■ 記事の範囲
この記事はテキスト第３部第５章の以下の節を取り扱います。
5.4 分析の準備
5.5 brmsによる単回帰モデルの推定
5.6 brmsの基本的な使い方
5.7 事前分布の変更
5.12 brmsによる事後分布の可視化
5.13 brmsによる予測
5.14 補足：predict関数を使わない予測の実装
5.15 回帰直線の図示
ベイズ統計モデリングの理論面が気になった場合には、ぜひテキストの第１部、第３部第１章をご覧いただき、本記事との繋がりをご確認下さいませ。
■ コード記述法
Jupyter Notebook 形式でコードを記述します。
■ 利用データ
テキスト・サポートサイトのデータファイルを引用しています。
Jupyter Notebook ファイルと同一フォルダ内の「data」フォルダにデータファイルを格納しています。
■ ライブラリのインポート
この記事で用いるライブラリをインポートします。
# インポート
# 数値計算
import numpy as np
import pandas as pd
# ベイズ統計モデリング
import bambi as bmb # formula文によるモデル定義
import arviz as az # 分析・可視化
# 統計処理
import scipy.stats as stats # 正規分布乱数の生成
# 可視化
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_theme() # ggplot風のスタイル
plt.rcParams['font.family'] = 'Meiryo'第５章 Bambiの使い方
Bambiの概要
Bambi は PyMC と連携したベイズ統計モデル構築ライブラリです。
Bambi 公式サイトの紹介文章をご覧ください。
Bambiは、Pythonで書かれた高水準のベイズモデル構築インターフェースです。確率プログラミングフレームワーク PyMC と連携し、生物学、社会科学、その他の分野で一般的に用いられるベイズ混合効果モデルの適合を極めて容易にするように設計されています。
公式サイトは次のリンクで！
公式サイト（部分、ブラウザで翻訳）の一部をご紹介します。
インストールは pip または conda を用いて実行するようです。
この記事で Bambi の使い方を確認していきましょう。
🦌🦌🦌
データの読み込み
テキスト 5.2 節に相当します。
テキストの仮想のビールの売り上げデータを引用いたします。
前々々回記事から利用しているデータです。
csv ファイルを pandas のデータフレーム形式で変数 file_beer_sales_2 に読み込みます。
# p.186 分析対象データの読み込み
# ファイルの読み込み
file_beer_sales_2 = pd.read_csv('./data/3-2-1-beer-sales-2.csv')
# 結果の表示
print('file_beer_sales_2.shape: ', file_beer_sales_2.shape)
file_beer_sales_2.head(3)【実行結果】
標本サイズ 100、変数 sales は 売り上げ（単位：万円）、temperature は気温（おそらく摂氏℃）です。
気温とビール売り上げの関係を単回帰モデルで分析します。
データの要約統計量を確認します。
Bambi モデリング時に参照しましょう。
# データの要約統計量
file_beer_sales_2.describe().T.round(3)【実行結果】
🦌🦌🦌
Bambiによる単回帰モデルの推定
テキスト 5.5 節に相当します。
早速ですが、次の formula 構文を用いて、Bambi で目的変数 sales、説明変数 temperature の単回帰モデルを構築しましょう。
$$
\mathtt{sales} \sim \mathtt{temperature}
$$
① ベイズモデリング
Bambi は短縮名 bmb で import しています。
$${\mathtt{bmb.Model()}}$$ でモデルの定義は完了です！
# モデリング
simple_lm_bmb = bmb.Model(
formula='sales ~ temperature', # フォーミュラ式
data=file_beer_sales_2, # データ
# family='gaussian', # 確率分布（デフォルトは正規分布gaussian）
# link='identity', # リンク（省略可）
)【実行結果】なし
【コードの補足】
Python の統計ライブラリ statsmodels の formula API のコードに似ています。
引数に formula と data（=データフレーム名）を設定するだけ！
目的変数の確率分布に正規分布を用いる場合、確率分布の設定は不要です。
また確率分布のデフォルトリンク関数を用いる場合、リンク関数の設定は不要です。
【ちなみに情報】
Bambi は formula の処理で formulae ライブラリを利用しています。
formulae サイトには Bambi のために作成されたと記載されています。
② モデルの表示
Bambi モデル名 simple_lm_bmb を実行するとモデルの内容を表示します。
# モデルの表示
simple_lm_bmb【実行結果】
モデル定義（インプット）に比べて、モデル内容（アウトプット）は情報が多いですね…
読み方を Gemini に訊いてみましょう。
Bambiが自動生成した各項目の意味と、背後にあるロジックを詳細に解説します。
モデルの基本定義
Formula：$${\mathtt{sales} \sim \mathtt{temperature}}$$
temperature（気温）を説明変数として、sales（売上）を予測する線形モデルです。
Family：$${\mathtt{gaussian}}$$ / Link：$${\mathtt{mu = identity}}$$
正規分布（ガウス分布）を誤差構造に採用しています。
リンク関数が identity（恒等関数）であるため、予測値 $${\mu}$$ は $${\mathtt{Intercept + temperature * beta}}$$ という線形結合の形をそのまま取ります。
Common-level effects（共通効果）の事前分布
これらは、データ全体に対して推定される固定的な係数です。Bambiはデータのスケールを見て、以下の正規分布（Normal）を自動設定しています。Intercept（切片）：$${\mathtt{Normal(mu: 69.938, sigma: 193.1861)}}$$
ロジック：切片の平均（mu）は、通常、目的変数 sales の平均値付近に設定されます。
解釈：データの中心が 69.938 付近にあるため、そこを基準に探索を開始します。標準偏差（sigma）が 193.1861 と非常に大きく設定されているのは、特定の値を決め打ちせず、広い範囲を許容する「弱情報事前分布」にするためです。
temperature（傾き）：$${\mathtt{Normal(mu: 0.0, sigma: 9.3318)}}$$
ロジック：係数の平均（mu）はデフォルトで 0.0 に設定されます。これは「気温が売上に影響を与えない（無相関）」という帰無仮説に近い状態をスタート地点とするためです。
解釈：標準偏差 9.3318 は、説明変数と目的変数の比率に基づいて計算されています。これにより、極端に大きすぎる傾き（気温が1度上がると売上が100万増えるなど）を抑制しつつ、データから最適な傾きを推定させます。
Auxiliary parameters（補助パラメータ）
モデルの線形予測子（$${\mu}$$）以外のパラメータです。
sigma（残差の標準偏差）：$${\mathtt{HalfStudentT(nu: 4.0, sigma: 22.177)}}$$役割：実際のデータが予測線（$${\mu}$$）からどれだけバラついているかを表す「誤差の大きさ」です。
分布の選択：標準偏差は必ず正の値になるため、0 以上の値のみを取る HalfStudentT（半t分布）が採用されています。
ロジック：
nu: 4.0（自由度）：裾の長い分布にすることで、外れ値に対して比較的頑健（ロバスト）な推論を可能にします。
sigma: 22.1772：目的変数 sales の標準偏差に基づいて計算されています。誤差がこの程度の範囲に収まるだろうという緩やかな制約を与えています。
◆ まとめ：Bambiが何をしたか
Bambiはデータ（100サンプル）をスキャンし、「売上の平均は約70で、標準偏差は22くらいだ」という情報を読み取りました。
その情報を元に、
切片は売上の平均（約70）付近。
傾きはとりあえず0（影響なし）から開始し、幅を持たせる。
誤差（sigma）は売上のバラつき（約22）を上限の目安にする。
という、「データに対して矛盾しないが、特定の偏見も持たない」絶妙な初期設定を自動で行ったことになります。これにより、ユーザーは事前分布を細かく指定せずとも、精度の高いサンプリング（MCMC）を実行できるようになっています。
モデルの定義で全く触れていないことを、Bambi が丁寧に設定していることが分かりました。
特にパラメータの事前分布は「弱情報事前分布」となっています。
弱情報事前分布は「推論結果をデータの情報に委ねつつ、現実離れした極端な値（異常な推定値）を排除するために、あえて広い裾野（緩やかな制約）を持たせた事前分布」」（Gemini談）です。
🦌
③ グラフィカルモデルの描画
PyMC と同様にグラフィカルモデル形式でモデルを可視化できます。
モデルに対してまず bulid メソッドで PyMC モデルオブジェクトを生成してから、graph メソッドで可視化します。
# モデルの可視化
simple_lm_bmb.build()
simple_lm_bmb.graph()【実行結果】
綺麗な単回帰モデルが描かれています！
④ 事前分布の可視化
「事前分布」サンプルを生成して、パラメータの事前分布を可視化します。
モデルに対して plot_priors メソッドを適用します。
# 事前分布プロットの描画
simple_lm_bmb.plot_priors();【実行結果】
事前分布の形状はいかがでしょう？（違和感無し？）
⑤ MCMC の実行
MCMCを実行しましょう。
モデルに対して fit メソッドを適用します。
NUTS サンプラーに nutpie を利用します。
引数は PyMC と同等です。
%%time
# p.187 MCMCの実行 ※nutpieのvar_namesワーニングは無視する...
idata_bmb1 = simple_lm_bmb.fit(
draws=1000, tune=1000, chains=4, random_seed=123, nuts_sampler='nutpie')【実行結果】（nutpie ワーニングは省略）
Divergences（ダイバージェンス）は０件です。
⑥ 要約統計量の表示
PyMC と同様に arviz を利用します。
# p.187 要約統計量の表示
var_names = ['Intercept', 'temperature', 'sigma']
az.summary(idata_bmb1, var_names=var_names, hdi_prob=0.95)【実行結果】
$${\widehat R}$$（r_hat）、有効サンプル数（ess_bulk、ess_tail）に問題はなさそうです。
前回記事の PyMC の単回帰モデルの要約統計量と比べると、Bambi モデルの方が有効サンプル数が圧倒的に大きい結果になっています。
（参考：前回記事のデザイン行列を用いた単回帰モデルの要約統計量）
⑦ トレースプロットの描画
PyMC と同様に arviz を利用します。
# p.188 トレースプロットの描画
az.plot_trace(idata_bmb1, compact=False, var_names=var_names,
backend_kwargs={'tight_layout': True});【実行結果】
右側のチャートがゲジゲジしています。
以上のチェックに基づいて、収束していると考えましょう。
Bambi によるベイズ統計モデリングの基本ステップを確認できました！
🦌🦌🦌
FamilyとLinkとFormula
テキスト 5.6 節に相当します。
一般化線形モデルを構成する主要３要素は、目的変数の確率分布、リンク関数、線形予測子です。
① モデル定義と３要素
３要素は Bambi のモデル定義で指定できます。
確率分布は family、リンク関数は link、線形予測子は formula の各引数で指定します。
② Bambi の family と link
Bambi にプリセットされている確率分布 と 確率分布のデフォルトリンク関数の一部を公式サイトから引用いたします。
multinomial の下にも一覧は続きますし、サンプルコードのリンクが「Example notebook」にありますので、ぜひ、次のリンクで公式サイトを訪ねてみましょう。
③ ３要素のクラス
複雑な線形予測子の記述やカスタム要件を実装する際には、$${\mathtt{Family}}$$ クラス、$${\mathtt{Link}}$$ クラス、$${\mathtt{Formula}}$$ クラスを利用します。
公式サイトには $${\mathtt{Formula}}$$ クラスを利用して formula を２つ設定する例が掲載されています。
構築したいモデルの内容によって、複雑な formula が描写されるようです。
詳しくは公式サイトをご確認下さい。
④ Formula クラスを使ってみる
$${\mathtt{Formula}}$$ クラスを用いて、テキストの bf 関数を真似てみましょう。
さきほどのモデリングで用いた formula 構文を使って、単回帰モデルを再度実装します。
# 複雑なformula文はbmb.Formula関数で作成
simple_lm_formula = bmb.Formula('sales ~ temperature')
simple_lm_formula【実行結果】
続いてどんどんモデリングを進めます。
結果は先ほどの単回帰モデルと全く同じになります。
・モデリング
# モデリング
simple_lm_bmb2 = bmb.Model(
formula=simple_lm_formula, # bmb.Formula形式のフォーミュラ式
data=file_beer_sales_2, # データ
)【実行結果】なし
・モデルの表示
# モデルの表示
simple_lm_bmb2【実行結果】
・モデルの可視化
# モデルの可視化
simple_lm_bmb2.build()
simple_lm_bmb2.graph()【実行結果】
・MCMC の実行
%%time
# MCMCの実行 ※nutpieのvar_namesワーニングは無視する...
idata_bmb2 = simple_lm_bmb2.fit(
draws=1000, tune=1000, chains=4, random_seed=123, nuts_sampler='nutpie')【実行結果】（nutpie ワーニングは省略）
・要約統計量の表示
# 要約統計量の表示
var_names = ['Intercept', 'temperature', 'sigma']
az.summary(idata_bmb2, var_names=var_names, hdi_prob=0.95)【実行結果】
・トレースプロットの描画
# トレースプロットの描画
az.plot_trace(idata_bmb2, compact=False, var_names=var_names,
backend_kwargs={'tight_layout': True});【実行結果】
🦌🦌🦌
Bambi の事前分布の変更
テキスト 5.7 節に相当します。
Bambi が自動設定する弱情報事前分布を使わずに、個別に事前分布を指定したい場合には、
$${\mathtt{Prior}}$$ クラスで事前分布を定義する
モデル定義時に prior 引数で定義した事前分布を指定する
をコード化します。
以下の幅の広い正規分布を事前分布に設定してモデリングしてみます。
$$
\begin{align*}
Intercept &\sim \text{Normal}\ (0,\ (1e5)^2) \\
temperature &\sim \text{Normal}\ (0,\ (1e5)^2) \\
\sigma &\sim \text{HalfNormal}\ ((1e5)^2) \\
\end{align*}
$$
・モデリング
# モデリング
# 事前分布の設定
uninformed_prior = bmb.Prior('Normal', mu=0, sigma=1e5) # 係数の無情報事前分布
sigma_prior = bmb.Prior('HalfNormal', sigma=1e5) # 標準偏差の無情報事前分布
priors = {'Intercept': uninformed_prior, # パラメータごとの事前分布を
'temperature': uninformed_prior, # 辞書にまとめる
'sigma': sigma_prior}
# モデルの定義
simple_lm_bmb3 = bmb.Model(
formula='sales ~ temperature', # フォーミュラ式
data=file_beer_sales_2, # データ
priors=priors, # 事前分布（辞書）
)【実行結果】なし
【コードの補足】
事前分布の定義は $${\mathtt{bmb.Prior}}$$ クラスを用いて、引数に確率分布の名称、当該確率分布のパラメータを設定します。
uninformed_prior = bmb.Prior('Normal', mu=0, sigma=1e5) # 係数の無情報事前分布モデルに与える準備として、パラメータと事前分布の対応を Python の辞書形式でまとめます。
priors = {'Intercept': uninformed_prior,
'temperature': uninformed_prior,
'sigma': sigma_prior}モデルの定義 $${\mathtt{bmb.Model}}$$ では priors 引数に辞書形式でまとめた事前分布の設定を与えます。
simple_lm_bmb3 = bmb.Model(
...(途中省略)...
priors=priors,
)以降は淡々と Bambi のコードを動かしていきます。
・モデルの表示
# モデルの表示
simple_lm_bmb3【実行結果】
$${\mathtt{Intercept,\ temperarure,\ sigma}}$$ の事前分布には先ほど設定した内容が表示されています。
・モデルの可視化
# モデルの可視化
simple_lm_bmb3.build()
simple_lm_bmb3.graph()【実行結果】
・MCMC の実行
%%time
# MCMCの実行 ※nutpieのvar_namesワーニングは無視する...
idata_bmb3 = simple_lm_bmb3.fit(
draws=1000, tune=1000, chains=4, random_seed=123, nuts_sampler='nutpie')【実行結果】（nutpie ワーニングは省略）
・要約統計量の表示
# 要約統計量の表示
var_names = ['Intercept', 'temperature', 'sigma']
az.summary(idata_bmb3, var_names=var_names, hdi_prob=0.95)【実行結果】
・トレースプロットの描画
# トレースプロットの描画
az.plot_trace(idata_bmb3, compact=False, var_names=var_names,
backend_kwargs={'tight_layout': True});【実行結果】
🦌🦌🦌
事後分布の可視化
テキスト 5.12 節に相当します。
テキストは係数の「信用区間」をフォレストプロットの形式で描画しています。
Bambi 特有のフォレストプロットは無さそうです。
ここではPyMC モデリングと同様に、arviz の plot_forest を利用して 係数の「95% HDI 区間」を描画します。
# p.196 図3.5.2 係数の95%HDI区間
az.plot_forest(idata_bmb1, var_names=['Intercept', 'temperature'],
hdi_prob=0.95, combined=True);【実行結果】
🦌🦌🦌
Bambi による予測
テキスト 5.13 節に相当します。
Bambi モデルの predict メソッドを利用して、以下の予測を行います。
パラメータの事後分布サンプルから算出する mu の予測値
事後予測サンプル（ばらつきを含む）に基づく sales の予測値
① 予測用の説明変数の作成
予測のための説明変数を作成します。
テキストと同じ 20 度を pandas のデータフレームに設定します。
# p.196 予測のための説明変数
new_data = pd.DataFrame({'temperature': [20]})
new_data【実行結果】
② mu の事後分布からのサンプリング
predict メソッドの引数 kind に「response_params」を設定します。
# p.196 muの事後分布からのサンプル（信用区間付きの予測値）
# デフォルトの動き：結果をidataのposteriorグループに追加する
ipred_post = simple_lm_bmb.predict(
idata_bmb1, # idata
data=new_data, # 予測用のデータ
inplace=False, # 結果をidata内に上書き保存するか
kind='response_params', # response_params：事後分布からのサンプル
)【実行結果】なし
ipred_post のグループ posterior 配下に mu の予測値が格納されました。
ipred_post.posterior.mu【実行結果】
③ mu の事後分布サンプルの平均値と 95% HDI
arviz の summary 関数を用いて、mu の事後分布サンプルの統計量を表示します。
# p.196 統計量の表示（回帰直線のHDI区間付きの予測値）
az.summary(ipred_post.posterior.mu, hdi_prob=0.95, kind='stats')【実行結果】
④ sales の事後予測サンプリング
predict メソッドの引数 kind に「response」を設定します。
# p.197 salesの予測分布からのサンプル（予測区間付きの予測値）
# デフォルトの動き：結果をidataのposterior_preditiveグループに追加する
ipred_pred = simple_lm_bmb.predict(
idata_bmb1, # idata
data=new_data, # 予測用のデータ
inplace=False, # 結果をidata内に上書き保存するか
kind='response', # response：事後予測分布からのサンプル
)【実行結果】なし
ipred_pred のグループ posterior_predictive 配下に sales の予測値が格納されました。
ipred_pred.posterior_predictive.sales【実行結果】
⑤ sales の事後予測サンプルの平均値と 95% HDI
arviz の summary 関数を用いて、sales の事後予測サンプルの統計量を表示します。
# p.197 統計量の表示（事後予測 HDI区間付きの予測値）
az.summary(ipred_pred.posterior_predictive.sales, hdi_prob=0.95, kind='stats')【実行結果】
mu よりも 95% HDI 区間が広いことが分かります。
🦌🦌🦌
Bambi や PyMC の予測関数を使わない予測の実装
テキスト 5.14 節に相当します。
MCMC サンプルを活用して、Bambi や PyMC の外で mu の事後分布サンプルや sales の事後予測サンプルを作成します。
① MCMC サンプルの取り出し
最初のモデルから MCMC サンプルを取り出して numpy 配列に変換します。
まず、arviz の extract 関数を利用して、idata のMCMCサンプルの次元を「フラット」にします（まだ xarray データセットのままです）。
# p.197 MCMCサンプルを取り出す
# idataのchainとdrawを平坦にする
idata_extract = az.extract(idata_bmb1.posterior)
idata_extract【実行結果】
パラメータ Intercepr、tempereture、sigma の形状は sample=4000 の平坦な状態になりました。
３つのパラメータをそれぞれ numpy 配列に変換します。
# p.197 推定されたパラメータ別に保存しておく
mcmc_b_intercept = idata_extract.Intercept.data
mcmc_b_temperature = idata_extract.temperature.data
mcmc_sigma = idata_extract.sigma.data【実行結果】なし
② mu の事後分布サンプルの算出
予測に用いる説明変数：気温を 20 度とすると、線形予測子 mu は
$$
\mathtt{mu} = \mathtt{Intercept} + 20 \times \mathtt{temprerature}
$$
で計算できます。計算しましょう。
# p.197 気温20度のときの売上の予測値
saigen_fitted = mcmc_b_intercept + 20 * mcmc_b_temperature
saigen_fitted【実行結果】
③ mu の再現（fitted の再現）
テキスト p.198 の fitted の再現に掲載の「mu の事後平均」「95% 信用区間」に加えて、arviz の summary 関数による要約統計量（HDI区間使用）を計算・表示します。
# p.198 fittedの再現
print('mean :', saigen_fitted.mean())
print('95%信用区間:', np.quantile(saigen_fitted, q=[0.025, 0.975]))
display(az.summary(saigen_fitted, hdi_prob=0.95, kind='stats'))【実行結果】
要約統計量は Bambi の predict 関数の要約統計量と一致しています！
（参考：Bambi の predict による mu の事後分布サンプルの統計量）
④ sales の事後予測サンプルの取得
sales の事後分布サンプルは、sales が従う正規分布から乱数を取得して用意します。
mu の事後分布サンプルの再現値 saigen_fitted を平均パラメータ、sigma のMCMC サンプルを標準偏差パラメータにして正規分布乱数を生成します。
# p.198 予測分布のMCMCサンプルを得る
saigen_predict = stats.norm.rvs(loc=saigen_fitted.mean(), scale=mcmc_sigma,
random_state=1)
saigen_predict【実行結果】
⑤ sales の再現（predict の再現）
テキスト p.198 の predict の再現に掲載の「sales の事後予測平均」「95% 予測区間」に加えて、arviz の summary 関数による要約統計量（HDI区間使用）を計算・表示します。
# p.198 predictの再現 ※乱数が異なるので再現できません
print('mean :', saigen_predict.mean())
print('95%予測区間:', np.quantile(saigen_predict, q=[0.025, 0.975]))
display(az.summary(saigen_predict, hdi_prob=0.95, kind='stats'))【実行結果】
要約統計量は Bambi の predict 関数の要約統計量と一致して⋯
いません（泣）
（参考：Bambi の predict による sales の事後予測サンプルの統計量）
【独白！再現性の裏にはあの乱数シードが⋯】
Bambi の predict 時に実は⋯
乱数シードを固定していませんでした。
自環境の Bambi バージョン 0.15.0 は predict の際に乱数シードを設定できないのです⋯
しかし朗報です！
Bambi バージョン 0.16.0 以降は、predict に引数 random_seed が追加されました！
きっと乱数シードを固定でき、再現可能になると思います（希望）！
🦌🦌🦌
Bambi による回帰直線の描画
テキスト 5.15 節に相当します。
さてさて気を取り直して、楽しい可視化タイムに進みましょう！
Bambi の可視化関数でスッキリしましょう！
最初に、mu の事後分布の「平均値＝回帰直線」と 「95% 信用区間の塗りつぶし」の描画です。
Bambi の interpret.plot_predictions 関数を利用します。
引数 pps=False で事後平均値、use_hdi=False で信用区間を指定します。
散布図部分は matplotlib の scatter 関数を利用します。
# p.199 図3.5.3 回帰直線：95%信用区間区間付き
# 描画領域の設定
fig, ax = plt.subplots(figsize=(8, 4))
# 観測値の散布図の描画
ax.scatter(file_beer_sales_2['temperature'], file_beer_sales_2['sales'],
color='navy', s=10, zorder=10)
# ★mu事後分布サンプルの平均値（回帰直線）と95%信用区間の塗りつぶし
bmb.interpret.plot_predictions(
simple_lm_bmb, # モデル
idata_bmb1, # idata
conditional='temperature', # 条件付けする共変量（変化させる変数）
pps=False, # 事後予測サンプルのプロットかどうか
use_hdi=False, # True: HDI、False: 分位数（信用区間）
prob=0.95, # 区間の確率
ax=ax
)
# 修飾
ax.set_xlabel('温度 [℃]', fontsize=12)
ax.set_ylabel('売上', fontsize=12)
ax.set_xticks(range(10, 31, 2));【実行結果】
はい！サクッと描画できました！
続いて、sales の事後予測の平均値曲線と 「95% 予測区間の塗りつぶし」の描画です。
mu のケースと同様に、Bambi の interpret.plot_predictions 関数を利用します。
事後予測の場合は、引数 pps=True を指定します。
散布図部分は matplotlib の scatter 関数を利用します。
# p.200 図3.5.4 回帰直線95%予測区間付き
# 描画領域の設定
fig, ax = plt.subplots(figsize=(8, 4))
# 観測値の散布図の描画
ax.scatter(file_beer_sales_2['temperature'], file_beer_sales_2['sales'],
color='navy', s=10, zorder=10)
# ★sales事後予測サンプルの平均値と95%予測区間の塗りつぶし
bmb.interpret.plot_predictions(
simple_lm_bmb, # モデル
idata_bmb1, # idata
conditional='temperature', # 条件付けする共変量（変化させる変数）
pps=True, # 事後予測サンプルのプロットかどうか
use_hdi=False, # True: HDI、False: 分位数（信用区間）
prob=0.95, # 区間の確率
ax=ax
)
# 修飾
ax.set_xlabel('温度 [℃]', fontsize=12)
ax.set_ylabel('売上', fontsize=12)
ax.set_xticks(range(10, 31, 2));【実行結果】
両方ともに、明示的にサンプル算出のコードを書いておらず、単純に描画関数の設定だけで可視化できました。
モデルのチェック時に手軽に利用できますね！
以上で Bambi による単回帰モデルの最初の学びを終わりにします。
PyMC のモデリングよりもシンプルなコードで、PyMC と同等の分析ができるのはとても有意義です！
今回の記事は以上です。
楽しかったですね！
シリーズの記事
次の記事（PyMC版）
前の記事（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の教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
いいなと思ったら応援しよう！
応援ありがとうございます。これからもがんばって記事を作成します！