「ベイズ統計モデリングによるデータ分析入門」をPythonとPyMCで写経 ～ Vol.14 ポアソン回帰モデル

書籍の著者　馬場真哉 先生

この記事は、書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」第３部第８章「ポアソン回帰モデル」の Python 写経活動記録です。

書籍の第３部は「一般化線形モデル」のベイズ統計モデリングです。
今回は「ポアソン回帰モデル」を３つのツール・手法で取り組みます。

• GLMによるポアソン回帰モデル

• Bambiのベイズポアソン回帰モデル

• PyMCのベイズポアソン回帰モデル

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

はじめに

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

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

準備

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

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

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

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

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

# インポート

# 数値計算
import numpy as np
import pandas as pd

# 統計処理
import scipy.stats as stats            # ポアソン分布乱数の生成

# ベイズ統計モデリング
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

# 可視化
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_theme()                        # ggplot風のスタイル
plt.rcParams['font.family'] = 'Meiryo'

第８章 ポアソン回帰モデル

ポアソン回帰モデルとは？

ポアソン回帰モデルは一般化線形モデルの文脈において、

• 線形予測子：複数の説明変数を用いる

• リンク関数：対数関数

• 確率分布：ポアソン分布

で構成される統計モデルです。
目的変数が０以上の整数の場合に用いられます。

🚀🚀🚀

データの読み込みと外観の確認

テキスト 8.4 節に相当します。
テキストの仮想の魚の釣獲尾数データを引用いたします。
csv ファイルを pandas のデータフレーム形式で変数 fish_num_climate に読み込みます。

# p.213 分析対象データの読み込み

# ファイルの読み込み
fish_num_climate = pd.read_csv('./data/3-8-1-fish-num-1.csv')

# 結果の表示
print('fish_num_climate.shape: ', fish_num_climate.shape)
fish_num_climate.head(3)

【実行結果】
標本サイズ 100、変数 fish_num は 釣った魚の数（釣獲尾数）、weather は天気の種類（晴れ、曇り）、temperature は気温です。
天気・気温と釣獲尾数の関係をポアソン回帰モデルで分析します。

天気ごとのデータ件数をカウントします。
pandas の value_counts メソッドを利用します。

# 天気ごとのデータ件数
fish_num_climate['weather'].value_counts().to_frame()

【実行結果】
曇り、晴れそれぞれ 50 件のデータがあります。

データの要約統計量を確認します。
まずは全体（量的変数のみ）です。

# データの要約統計量
fish_num_climate.describe().T.round(2)

【実行結果】
釣獲尾数は平均 1.6、最小値 0、最大値 8 です。
気温は平均 14.8、最小値 0.2、最大値 29.7 です。

次は天気別釣獲尾数の要約統計量です。

# 天気別の釣獲尾数の要約統計量
fish_num_climate.groupby(['weather'])['fish_num'].describe().round(2)

【実行結果】
曇りの方が大きい感じがします。

次は天気別気温の要約統計量も念のため見ておきましょう。

# 天気別の気温要約統計量
fish_num_climate.groupby(['weather'])['temperature'].describe().round(2)

【実行結果】
曇りの方が気温が高い感じがします。

データを可視化しましょう。
最初にジョイントプロットで、要約統計量の感覚値を可視化してみます。
seaborn の jointplot を利用します。

# ジョイントプロットの描画
colors = {'cloudy': 'tomato', 'sunny': 'tab:blue'}
sns.jointplot(
    data=fish_num_climate, x='temperature', y='fish_num',
    hue='weather', palette=colors
);

【実行結果】
散布図では気温の上昇と釣獲尾数の増加＆ばらつき増加の関係が見られます。
気温のKDEプロット（上）では、天気別気温の分布は似ています。
釣獲尾数のKDEプロット（右）では、晴れのばらつきが小さく、曇りのばらつきが大きいことが分かります。

続いて、テキスト 図 3.8.1 に相当する散布図です。
seaborn の scatterplot を利用します。

# p.214 図3.8.1 釣獲尾数と天気・気温の散布図

# 描画領域の設定
plt.figure(figsize=(8, 4))
# 散布図
colors = {'cloudy': 'tomato', 'sunny': 'tab:blue'}
sns.scatterplot(data=fish_num_climate, x='temperature', y='fish_num',
                hue='weather', palette=colors)
# 修飾
plt.title('釣獲尾数と気温・天気の関係', loc='left');

【実行結果】
ジョイントプロットで見てしまいましたが…気温の上昇と釣獲尾数の増加＆ばらつきの増加の関係が見られます。

この散布図に天気ごとの回帰曲線を重ねてみましょう。
seaborn の lmplot を利用します。
lowess 引数を ON にして局所的な傾きの違いを表現します。

# 釣獲尾数と天気・気温の散布図（回帰曲線付き）

# 回帰曲線付き散布図
colors = {'cloudy': 'tomato', 'sunny': 'tab:blue'}
sns.lmplot(data=fish_num_climate, x='temperature', y='fish_num', lowess=True,
           hue='weather', palette=colors, height=4, aspect=1.8)
# 修飾
plt.title('釣獲尾数と気温・天気の関係（回帰曲線付き）', loc='left');

【実行結果】
すべての天気で気温上昇と釣獲尾数増加の関係が見られます。
また曇りの釣獲尾数は相対的に大きいです。

🚀🚀🚀

GLMによるポアソン回帰モデル

① モデルの概要
まずGLMでポアソン回帰モデルを確認しておき、あとでベイズ統計モデルの結果と比べてみましょう。
Python の統計ライブラリ statsmodels を利用して、次のモデルを実装します！

$${\eta}$$ は線形予測子、$${\beta}$$ は係数、$${x_{i1}}$$ は晴れダミー、$${x_{i2}}$$ は気温です。
$${\log(\lambda)}$$ は対数リンク関数、$${\exp(\eta)}$$ は逆リンク関数（指数関数）です。
$${\lambda}$$ はポアソン分布の平均パラメータ、$${y}$$ は目的変数となる釣獲尾数です。

このポアソン回帰モデルの線形予測子を statsmodels の一般化線形モデル glm にあたえる formula 構文に変換します。
$${\sim}$$ を挟んで、左辺は目的変数、右辺は説明変数です。

② ポアソン回帰の実行
ではポアソン回帰を実行します！
確率分布 family で sm.families.Poisson() を指定します。
ポアソン分布のデフォルトリンク関数は対数リンク関数です（指定不要）。

# ポアソン回帰モデル by statsmodels
formula = 'fish_num ~ weather + temperature'   # formulaの設定
family=sm.families.Poisson()                   # 確率分布とリンク関数の設定
res_sm = smf.glm(formula=formula, data=fish_num_climate, family=family).fit()
res_sm.summary()

【実行結果】
こちらは GLM のサマリーです。

③ 推定値の確認
係数の推定値にフォーカスしてみます。

# 要約表から係数の推定値の部分を取り出し
res_sm.summary().tables[1]

【実行結果】

ポアソン回帰の線形予測子の係数の解釈は少々ややこしいです。
ポアソン分布の平均パラメータ $${\lambda}$$ と線形予測子の関係を解きほぐします。

平均パラメータ $${\lambda}$$ は $${\exp(\beta)}$$ の掛け算になります。

曇りと比べて晴れの場合は $${\exp(\beta_{1, 晴れ})=\exp(-0.5869) \approx 0.56}$$ 倍です。曇りの半分近くになるようです。

気温が１度上昇する場合には：

となり、$${\exp(\beta_2)=\exp(0.0827) \approx 1.09}$$ 倍です。
気温１度上昇で $${1.09}$$ 倍（$${9\%}$$）釣獲尾数が増えるようです。

「$${\exp(係数の推定値)}$$ 倍」と覚えておきましょう。

④ 天気別気温別の平均釣獲尾数の可視化
最後に天気別気温別の平均釣獲尾数の95%信頼区間を可視化しましょう。

最初に平均釣獲尾数の予測値を算出します。
ポアソン回帰の結果 res_sm に対して get_prediction メソッドを適用するなどして、予測値を得ます。

## 平均釣獲尾数予測値の算出
# 予測用データの作成
n_samples = 100
test_data = pd.DataFrame({
    'weather': np.repeat(['cloudy', 'sunny'], n_samples),     # 天気
    'temperature': np.tile(np.linspace(0, 30, 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 = {'cloudy': 'tomato', 'sunny': 'tab:blue'}
# 描画領域の設定
plt.figure(figsize=(8, 4))
# 観測値の散布図の描画
sns.scatterplot(data=fish_num_climate, x='temperature', y='fish_num',
                hue='weather', palette=colors)
# 天気ごとに平均釣獲尾数の予測値と95%信頼区間の描画を繰り返し処理
for weather, color in colors.items():
    # 描画データの取得
    temp, mean_sales, lower, upper = (
        pred_df[pred_df['weather']==weather].iloc[:, [1, 2, 4, 5]].to_numpy().T
    )
    # 平均釣獲尾数の予測値（点推定）の描画
    plt.plot(temp, mean_sales, color=color)
    # 平均釣獲尾数の95%信頼区間の塗りつぶし描画
    plt.fill_between(temp, lower, upper, color=color, alpha=0.2)
# 修飾
plt.title('釣獲尾数の95%信頼区間')
plt.xlabel('気温', fontsize=12)
plt.ylabel('釣獲尾数', fontsize=12);

【実行結果】
気温が大きくなるにつれて釣獲尾数の増加数が大きくなる＝指数関数的に増えていることが分かります。
曇りの方が釣獲尾数が大きい感じがします。
また気温が大きくなるにつれて 95% 信頼区間の幅が広くなり、ばらつきが大きくなっていることを示しています。

ベイズ流のポアソン回帰モデルへ進みます。

🚀🚀🚀

ベイズモデリング by Bambi

テキスト 8.5、8.7 節に相当します。
brms の代わりに Bambi を利用します。

① モデルの概要
確率分布をポアソン分布、リンク関数を対数リンク関数（ポアソン分布のデフォルト）、線形予測子を次の formula で実装します。

テキストの無情報事前分布に合わせるべく、Bambi に以下の事前分布情報を与えます。

② モデル定義
上式のモデルを Bambi で記述します。
bmb.prior() で事前分布を定義します。複数ある場合は辞書でまとめます。
bmb.Model() において、prior 引数で定義した事前分布を与えます。

# モデリング

# 無情報事前分布（想定）の設定
uninformed_prior = bmb.Prior('Normal', mu=0, sigma=1e5)
# 事前分布を辞書にとりまとめ
priors = {'Intercept': uninformed_prior, 'weather': uninformed_prior,
          'temperature': uninformed_prior}

# モデルの定義
model_bmb = bmb.Model(
    formula='fish_num ~ weather + temperature',   # フォーミュラ式
    data=fish_num_climate,                        # データ
    family='poisson',                             # 確率分布：ポアソン分布
    # link='log'                                    # リンク関数：対数（デフォルト）
    priors=priors,                                # 事前分布（辞書）
)

【実行結果】なし

モデルの内容を表示します。

# モデルの表示
model_bmb

【実行結果】
確率分布 Family に ポアソン分布 poisson、リンク関数 Link に対数リンク関数 mu = log が設定されました。

モデルをグラフィカルモデルで描画します。

# モデルの可視化
model_bmb.build()
model_bmb.graph()

【実行結果】
weather 係数の次元は１です（weather_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.215 要約統計量の表示
var_names = ['Intercept', 'weather', 'temperature']
az.summary(idata_bmb, var_names=var_names, hdi_prob=0.95)

【実行結果】
$${\widehat{R}}$$（r_hat）、有効サンプル数（ess_bulk、ess_tail）に問題はなさそうです。

ちなみに $${\exp(-0.591) \approx 0.55}$$ 倍、$${\exp(0.083) \approx 1.09}$$ 倍です。

GLMの推定値の表を並べましょう。
（参考：statsmodels の推定値）

手法間で係数推定値は近い値になっています。

トレースプロットを描画します。

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

【実行結果】
左側のチャートの４本の Chain はほぼ重なっており、かつ１峰です。
右側のチャートがゲジゲジしています。

以上のチェックに基づいて、収束していると考えましょう。

⑤ 気温別天気別の平均釣獲尾数・予測釣獲尾数データの作成
Bambi モデル model_bmb に対して predict メソッドを適用することで平均釣獲尾数（mu の事後分布に相当）、予測釣獲尾数（y の事後予測に相当）のMCMCサンプルを取得します。

# 曇り、晴れの平均釣獲尾数、予測釣獲尾数の算出

## 設定
# 予測データの個数
n_pred = 100
# 気温の値
x_val = np.linspace(fish_num_climate['temperature'].min(),
                    fish_num_climate['temperature'].max(), n_pred)
# 天気の要素
weathers = ['cloudy', 'sunny']

## 予測の実行
results_mu_pred, results_y_pred = [], []
# 天気ごとに予測を繰り返し処理
for weather in weathers:
    # 予測データの作成
    new_data = pd.DataFrame({'weather': [weather]*n_pred, 'temperature': x_val})
    # muの事後予測サンプル（平均釣獲尾数）の取得
    result_mu_pred = model_bmb.predict(
        idata=idata_bmb, data=new_data, kind='response_params', inplace=False)
    # fish_numの事後予測サンプル（予測釣獲尾数）の取得
    result_y_pred = model_bmb.predict(
        idata=idata_bmb, data=new_data, kind='response', inplace=False)
    # 結果をリストに格納
    results_mu_pred.append(result_mu_pred)
    results_y_pred.append(result_y_pred)

【実行結果】なし

⑥ 気温別天気別の平均釣獲尾数の可視化
テキストの図 3.8.2 に相当します。
平均獲釣尾数の中央値と 95% ベイズ信用区間を描画します。

# p.216 図3.8.2 ポアソン回帰曲線：信用区間

## 平均釣獲尾数と描画色を天気別の辞書に格納
fishes = {weather: az.extract(res.posterior).mu.data.T
          for weather, res in zip(weathers, results_mu_pred)}
colors = {'cloudy': 'tomato', 'sunny': 'tab:blue'}

## 描画
# 描画領域の設定
plt.figure(figsize=(8, 4))
# 観測値の散布図の描画
sns.scatterplot(data=fish_num_climate, x='temperature', y='fish_num',
                hue='weather', palette=colors)
# 天気ごとに平均釣獲尾数の中央値と95%信用区間の描画を繰り返し処理
for weather in fishes.keys():
    # 平均釣獲尾数の中央値の描画
    plt.plot(x_val, np.median(fishes[weather], axis=0), color=colors[weather])
    # 平均釣獲尾数の95%信用区間の塗りつぶし描画
    ci95 = np.quantile(fishes[weather], q=[0.025, 0.975], axis=0).T
    plt.fill_between(x_val, ci95[:, 0], ci95[:, 1], color=colors[weather],
                     alpha=0.2)

# 修飾
plt.title('平均釣獲尾数の95%信用区間')
plt.xlabel('気温', fontsize=12)
plt.ylabel('釣獲尾数', fontsize=12);

【実行結果】
曇りの平均釣獲尾数は晴れと比べて大きいです。
また気温が大きくなるにつれて 95% 信用区間の幅が広くなり、ばらつきが大きくなっていることを示しています。
statsmodels の 95% 信頼区間とも似ている感じです。

（参考：statsmodels の95%信頼区間）

⑦ 気温別天気別の予測釣獲尾数の可視化
テキストの図 3.8.3 に相当します。
予測獲釣尾数の中央値と 99% 予測区間を描画します。
観測値のデータ点が 99% 予測区間に収まっていることを確認します。

# p.217 図3.8.3 ポアソン回帰曲線：予測区間

## 予測釣獲尾数と描画色を天気別の辞書に格納
fishes = {weather: az.extract(res.posterior_predictive).fish_num.data.T
          for weather, res in zip(weathers, results_y_pred)}
colors = {'cloudy': 'tomato', 'sunny': 'tab:blue'}

## 描画
# 描画領域の設定
plt.figure(figsize=(8, 4))
# 観測値の散布図の描画
sns.scatterplot(data=fish_num_climate, x='temperature', y='fish_num',
                hue='weather', palette=colors)
# 天気ごとに釣獲尾数の中央値と99%予測区間の描画を繰り返し処理
for weather in fishes.keys():
    # 釣獲尾数の中央値の描画
    plt.plot(x_val, np.median(fishes[weather], axis=0), color=colors[weather])
    # 釣獲尾数の99%予測区間の塗りつぶし描画
    ci99 = np.quantile(fishes[weather], q=[0.005, 0.995], axis=0).T
    plt.fill_between(x_val, ci99[:, 0], ci99[:, 1], color=colors[weather],
                     alpha=0.2)
# 修飾
plt.title('釣獲尾数の99%予測区間')
plt.xlabel('気温', fontsize=12)
plt.ylabel('釣獲尾数', fontsize=12)
plt.legend(loc='upper left', title='weather');

【実行結果】
観測値のデータ点が 99% 予測区間に収まっています。
ところで、平均釣獲尾数と比べると、とてもガタガタした形状です。
テキストによると「整数値をとるポアソン分布に従う乱数を使っている」からガタガタしているとのこと。

🚀🚀🚀

ベイズモデリング by PyMC

テキスト 8.7、8.8、8.9 節に相当します。
デザイン行列を作成して PyMC でモデリングします。

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

② データセットの作成
patsy ライブラリを用いてデザイン行列等を作成します。

# p.219 デザイン行列の作成

# formula構文の設定
formula_pois = 'fish_num ~ weather + temperature'

# 目的変数Y, デザイン行列（説明変数）Xの作成
Y_dm, X_dm = dmatrices(formula_pois, fish_num_climate, return_type='dataframe')

# デザイン行列の先頭5行の表示
X_dm.head()

【実行結果】
デザイン行列は定数項、天気のダミー変数（晴れ）、気温で構成されます。

# 目的変数の先頭5行の表示
Y_dm.head()

【実行結果】
目的変数も作ってくれました。

③ モデル定義
上式のモデルを PyMC で記述します。

# モデリング

# coordsの設定
coords = {'id': fish_num_climate.index.values,             # 観測データの識別子
          'coefs': ['intercept', 'sunny', 'temperature']}  # 係数ベクトルbの識別子

# モデルの定義
with pm.Model(coords=coords) as model_pm:
    
    ## dataの設定
    # 目的変数: 釣獲尾数データ
    Y = pm.Data('Y', value=Y_dm.values.flatten(), dims='id')
    # 説明変数: デザイン行列
    X = pm.Data('X', value=X_dm.values, dims=('id', 'coefs'))

    ## 事前分布: 無情報事前分布的な分布
    # 係数ベクトル
    b = pm.Normal('b', mu=0, sigma=1e5, dims='coefs')

    ## 尤度関数: ポアソン分布を仮定
    lam = pm.Deterministic('lam', pm.math.exp(X @ b), dims='id')
    obs = pm.Poisson('obs', mu=lam, 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(-0.587) \approx 0.56}$$ 倍、$${\exp(0.084) \approx 1.09}$$ 倍です。

（参考：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.8.2 に相当します。
平均獲釣尾数の中央値と 95% ベイズ信用区間を描画します。
パラメータのMCMCサンプルを使って、PyMCの外で平均釣獲尾数の予測値を計算します。

# p.216 図3.8.2 ポアソン回帰曲線：信用区間

## 準備
# MCMCサンプルからβ0, β1, β2を取り出し
intercept, w_sunny, temper = az.extract(idata_pm.posterior).b.values
# x軸の値
x_val = np.linspace(X_dm['temperature'].min(), X_dm['temperature'].max(), 100)

## 平均釣獲尾数データの作成
# MCMCサンプルから曇り、晴れの平均釣獲尾数λを算出
cloudy_lam = np.exp(intercept + w_sunny*0 + np.outer(x_val, temper)).T
sunny_lam = np.exp(intercept + w_sunny*1 + np.outer(x_val, temper)).T
# 平均釣獲尾数データと描画色を天気別の辞書に格納
fishes = {'cloudy': cloudy_lam, 'sunny': sunny_lam}
colors = {'cloudy': 'tomato', 'sunny': 'tab:blue'}

## 描画
# 描画領域の設定
plt.figure(figsize=(8, 4))
# 観測値の散布図の描画
sns.scatterplot(data=fish_num_climate, x='temperature', y='fish_num',
                hue='weather', palette=colors)
# 天気ごとに平均釣獲尾数の中央値と95%信用区間の描画を繰り返し処理
for weather in fishes.keys():
    # 平均釣獲尾数の中央値の描画
    plt.plot(x_val, np.median(fishes[weather], axis=0), color=colors[weather])
    # 平均釣獲尾数の95%信用区間の塗りつぶし描画
    ci95 = np.quantile(fishes[weather], q=[0.025, 0.975], axis=0).T
    plt.fill_between(x_val, ci95[:, 0], ci95[:, 1], color=colors[weather],
                     alpha=0.2)

# 修飾
plt.title('釣獲尾数の95%信用区間')
plt.xlabel('気温', fontsize=12)
plt.ylabel('釣獲尾数', fontsize=12);

【実行結果】
Bambi モデルの 95% 信用区間とよく似た形状になりました。

（参考：Bambiモデルの 95% 信用区間）

⑦ 気温別天気別の予測釣獲尾数の可視化
テキストの図 3.8.3 に相当します。
予測獲釣尾数の中央値と 99% 予測区間を描画します。
観測値のデータ点が 99% 予測区間に収まっていることを確認します。
パラメータのMCMCサンプルを使って、PyMCの外でポアソン分布乱数を生成して予測釣獲尾数を得ます。

# p.217 図3.8.3 ポアソン回帰曲線：予測区間

## 予測釣獲尾数データの作成
# 曇り、晴れの平均釣獲尾数λから予測値を算出：ポアソン分布乱数を生成
cloudy_pred = stats.poisson.rvs(mu=cloudy_lam, random_state=123)
sunny_pred = stats.poisson.rvs(mu=sunny_lam, random_state=123)
# 予測釣獲尾数データと描画色を天気別の辞書に格納
fishes = {'cloudy': cloudy_pred, 'sunny': sunny_pred}
colors = {'cloudy': 'tomato', 'sunny': 'tab:blue'}

## 描画
# 描画領域の設定
plt.figure(figsize=(8, 4))
# 観測値の散布図の描画
sns.scatterplot(data=fish_num_climate, x='temperature', y='fish_num',
                hue='weather', palette=colors)
# 天気ごとに釣獲尾数の中央値と99%予測区間の描画を繰り返し処理
for weather in fishes.keys():
    # 釣獲尾数の中央値の描画
    plt.plot(x_val, np.median(fishes[weather], axis=0), color=colors[weather])
    # 釣獲尾数の99%予測区間の塗りつぶし描画
    ci99 = np.quantile(fishes[weather], q=[0.005, 0.995], axis=0).T
    plt.fill_between(x_val, ci99[:, 0], ci99[:, 1], color=colors[weather],
                     alpha=0.2)
# 修飾
plt.title('釣獲尾数の99%予測区間')
plt.xlabel('気温', fontsize=12)
plt.ylabel('釣獲尾数', fontsize=12)
plt.legend(loc='upper left', title='weather');

【実行結果】
Bambi モデルの 99% 予測区間とそこそこ似た形状になりました。

（参考：Bambiモデルの 99% 予測区間）

🚀🚀🚀

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

シリーズの記事

次の記事

前の記事

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

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