「ベイズ統計モデリングによるデータ分析入門」をPythonとStanで写経 ～ Vol.18 ランダム切片モデル
書籍の著者 馬場真哉 先生
この記事は、書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」第４部第２章「ランダム切片モデル」の Python 写経活動記録です。
第４部は一般化線形混合モデル（通称：GLMM）の階層ベイズモデルです。
今回は「ランダム切片モデル」を２つのツール・手法で取り組みます。
Stanのベイズランダム切片モデル
アディショナルタイム：GPBoostのランダム切片モデル
では書籍を開いてベイズ統計モデリングの旅に出かけましょう🚀
はじめに
このブログシリーズは書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」（講談社、「テキスト」と呼びます）の Python 写経です。
テキストの紹介と引用表記はリンク先の記事に掲載しています。
準備
■ 記事の範囲
この記事はテキスト第４部第２章の以下の節を取り扱います。
2.2 分析の準備
2.5 brms によるランダム切片モデルの推定
2.6 回帰曲線の図示
ベイズ統計モデリングの理論面が気になった場合には、ぜひテキストの第１部および本章の理論面の記載をご覧いただき、本記事との繋がりをご確認下さいませ。
■ コード記述法
Jupyter Notebook 形式でコードを記述します。
■ 利用データ
テキスト・サポートサイトのデータファイルを引用しています。
Jupyter Notebook ファイルと同一フォルダ内の「data」フォルダにデータファイルを格納しています。
■ ライブラリのインポート
この記事で用いるライブラリをインポートします。
# インポート
# 数値計算
import numpy as np
import pandas as pd
# ベイズ統計モデリング
from cmdstanpy import CmdStanModel # stan
import arviz as az # 分析・可視化
# デザイン行列
from patsy import dmatrices
# ユーティリティ
import os
# 可視化
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_theme() # ggplot風のスタイル
plt.rcParams['font.family'] = 'Meiryo'第２章 ランダム切片モデル
前回記事から引き続き、ポアソンGLMMの階層ベイズモデルである「ランダム切片モデル」を実装します。
GLMMおよび階層ベイズモデルの概要はぜひ前回記事をお読み下さい！
★★★前回記事のリンク★★★
🔵🔵🔵
データの読み込みと外観の確認
テキスト 2.2 節に相当します。
テキストの仮想の魚の釣獲尾数データを引用いたします。
csv ファイルを pandas のデータフレーム形式で変数 fish_num_climate_3 に読み込みます。
# p.255 分析対象データの読み込み
# ファイルの読み込み
fish_num_climate_3 = pd.read_csv('./data/4-2-1-fish-num-3.csv')
# 結果の表示
print('fish_num_climate_3.shape: ', fish_num_climate_3.shape)
fish_num_climate_3.head(3)【実行結果】
標本サイズ 100のデータです。
fish_num：釣った魚の数（釣獲尾数）
weather：天気の種類（晴れ：sunny、曇り：cloudy）
temperature：気温
human：釣り人の名前（A ～J までの 10 人）
のちのち使う目的で共通データを作成します。
１つ目は釣り人の名前リストです。
# 釣り人の名前リスト
humans = sorted(fish_num_climate_3['human'].unique())
humans【実行結果】
A ～J までの 10 人の釣り人のデータです。
２つ目は天気ごとの色設定です。
# 天気ごとの色の設定
colors = {'cloudy': 'tomato', 'sunny': 'tab:blue'}
colors【実行結果】
曇りがトマト色（オレンジ寄りの赤）、晴れが青色です。
🔵
天気ごとおよび釣り人ごとのデータ件数をカウントします。
pandas の value_counts メソッドを利用します。
# カテゴリ変数の要素ごとのデータ件数
for col in fish_num_climate_3.columns[fish_num_climate_3.dtypes == 'object']:
display(fish_num_climate_3[col].value_counts().to_frame())【実行結果】
曇り・晴れで 50 件ずつ、釣り人ごとに 10 件ずつのデータです。
データの要約統計量を確認します。
まずは全体（量的変数のみ）です。
# データの要約統計量
fish_num_climate_3.describe().T.round(2)【実行結果】
釣獲尾数は平均 2.5、最小値 0、最大値 15 です。
気温は平均 15.6、最小値 0.4、最大値 29.7 です。
釣獲尾数にポアソン分布を仮定する際には、平均 2.5 と分散 9 （標準偏差の二乗）の乖離、つまり過分散の考慮が必要かもしれません。
次は釣り人別・天気別の釣獲尾数の要約統計量です。
# 釣り人 × 天気 ごとの要約統計量
fish_num_climate_3.groupby(['human', 'weather'])['fish_num'].describe().round(2)【実行結果】
さまざまにばらついている感じです。
数字を追うのは大変💦ですので、後で可視化しましょう。
データを可視化しましょう。
最初に散布図で気温と釣獲尾数の関係を可視化します。
seaborn の scatterplot を利用します。
# 天気別散布図の描画
sns.scatterplot(
data=fish_num_climate_3, x='temperature', y='fish_num',
hue='weather', palette=colors
)
# 修飾
plt.xlabel('気温 [℃]', fontsize=12)
plt.ylabel('釣獲尾数', fontsize=12)
plt.legend(title='天気');【実行結果】
散布図では気温が上昇するにつれて、釣獲尾数の増加傾向とばらつき増加傾向の関係が見られます。
続いて、箱ひげ図で釣り人別・天気別の釣獲尾数のばらつきを確認します。
seaborn の boxplot を利用します。
# 箱ひげ図の描画
# 描画領域の設定
plt.figure(figsize=(10, 4))
# 箱ひげ図の描画
sns.boxplot(data=fish_num_climate_3, x='human', y='fish_num',
hue='weather', palette=colors, fill=False, gap=0.2)
# スウォームプロットの描画
sns.swarmplot(data=fish_num_climate_3, x='human', y='fish_num',
hue='weather', palette=colors, dodge=True)
# 修飾
plt.xlabel('釣り人', fontsize=12)
plt.ylabel('釣獲尾数', fontsize=12)
plt.legend(title='天気');【実行結果】
釣り人によって釣れた魚の数のばらつきが大きく異なっている印象です。
また、曇りの方が釣果がいい人、晴れの方が釣果がいい人に分かれている感じもします。
最後に釣り人・気温と釣獲尾数の散布図を描画します。
seaborn の relplot で釣り人別の小さなチャートを配置します。
# 釣り人別散布図の描画
g = sns.relplot(
data=fish_num_climate_3, x='temperature', y='fish_num', # データ, x軸, y軸
col='human', col_wrap=3, # 列の変数, 列の数
hue='weather', palette=colors, # hueの変数、色
kind='scatter', s=80, height=3, aspect=1.2
)
# 修飾
g.set_titles(col_template='釣り人 {col_name}')
g.set_axis_labels('気温 [℃]', '釣獲尾数')
g._legend.set_title('天気');【実行結果】
A さん、C さんがよく釣れている感じ、D ～ F さんはあまり釣れていない感じ、天気と釣果の関係は…今ひとつ分かりません…。
ベイズモデルに期待しましょう！
なお、この relplot を用いた散布図の描画は、以後の分析でも頻出するので、お楽しみに🍀
🔵🔵🔵
ベイズモデリング by Stan
テキスト 2.5、2.6 節に相当します。
デザイン行列を作成して CmdStanPy でモデリングします。
① モデルの数式
今回は次の数式で示されるランダム切片モデルに取り組みます。
$${r_k}$$ が釣り人 $${k}$$ ごとのランダム効果（ランダム切片）です。
$$
\begin{align*}
r_k &\sim \text{Normal}(0, \sigma_r^2) \\
\log(\lambda_i) &= \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + r_k \\
y_i &\sim \text{Poisson}(\lambda_i)
\end{align*}
$$
💡 ランダム切片の気持ちを探る
ランダム切片 $${r_k}$$ は平均 $${0}$$、標準偏差 $${\sigma_r}$$（分散 $${\sigma_r^2}$$）の正規分布に従うとしています。
つまり、釣り人のランダム切片の効果は１つの正規分布からサンプリングされて決まるという関係になっています。
そして、標準偏差 $${\sigma_r}$$ が小さい場合は各釣り人のランダム効果は平均 $${0}$$ に近く、みなの釣果は似たりよったりになります。
標準偏差 $${\sigma_r}$$ が大きい場合は各釣り人のランダム効果はばらつきが大きく、個性的な釣果が生まれる可能性があります。
一見、ランダム切片は「ダミー変数の係数」に似ているように見えます。
しかし、ダミー変数はカテゴリーごとに独立して推定されるため、共通の確率分布の設計（または縛り）を持っていません。
Gemini はこう言ってます。
ダミー変数のようには、少ないデータで釣り人の実力を決めつけない。この『正規分布という設計図』の範囲内に収まるはずだと考えることで、極端な予測ミスを防ぐ。これこそが、ランダム効果モデルが持つ『賢さ』の正体なのです。
② モデルの概要
デザイン行列、係数ベクトル、ランダム切片を用いて、次のモデルを実装します。
$$
\begin{align*}
\bm Y &\sim \text{Poisson}(\bm \lambda) \\
\log(\bm \lambda) & = \bm {X \beta} + \bm r \\
\bm r &\sim \text{Normal}(0, \sigma_r^2) \\
\end{align*}
$$
目的変数 $${\bm Y}$$（fish_num）はポアソン分布に従うと仮定しています。
ポアソン分布の平均パラメータベクトルの対数 $${\log(\bm \lambda)}$$ は、デザイン行列 $${\bm X}$$ と係数ベクトル $${\bm \beta}$$ の積にランダム切片 $${\bm r}$$ を加えた線形予測子と等しいです。
ランダム切片 $${\bm r}$$ の次元は human の数 10 であり、正規分布に従うと仮定しています。
係数ベクトル $${\bm \beta}$$ と ランダム切片の標準偏差パラメータ $${\sigma_r^2}$$ には事前分布を明示的に設定しません。
③ Stan のモデル設定
Stan ファイル（Stan コード）を作成します。
デザイン行列を用いる書き方を採用しています。
ランダム切片 $${r}$$ の要素「釣り人」を表す標本と同じサイズの変数を human で定義します。
📑ファイル名：4-2-1-glmm-pois-design-matrix.stan
data {
int N; // 標本サイズ
int P; // デザイン行列の列数（説明変数の数＋１）
int K; // 釣り人の数
array[N] int Y; // 目的変数
matrix[N, P] X; // 説明変数
array[N] int human; // 釣り人
}
parameters {
vector[P] b; // 切片を含む係数ベクトル
vector[K] r; // ランダム効果
real<lower=0> sigma_r; // ランダム効果の標準偏差
}
transformed parameters {
vector[N] lam = X * b + r[human];
}
model {
r ~ normal(0, sigma_r);
Y ~ poisson_log(lam);
}【実行結果】なし
③ データセットの作成
patsy ライブラリを用いてデザイン行列等を作成します。
# デザイン行列の作成
# formula構文の設定
formula_pois = 'fish_num ~ weather + temperature'
# 目的変数Y, デザイン行列（説明変数）Xの作成
Y, X = formula_lm = dmatrices(formula_pois, fish_num_climate_3,
return_type='dataframe')
# int型の設定
X = X.astype({'Intercept': int, 'weather[T.sunny]': int})
Y = Y.astype(int)
# デザイン行列の先頭5行の表示
X.head()【実行結果】
デザイン行列は定数項、天気のダミー変数（晴れ）、気温で構成されます。
# 目的変数の先頭5行の表示
Y.head()【実行結果】
目的変数も作ってくれました。
標本サイズ・説明変数の数を算出して、Stan に渡すデータセットを辞書にまとめます。
# データセットの準備
# サンプルサイズ、デザイン行列の列数（説明変数の数＋１）
N, P = X.shape
# 釣り人 ※A,B,C... から 1,2,3...に変換
codes, uniques = pd.factorize(fish_num_climate_3['human'], sort=True)
human = codes + 1 # codesは0はじまりなので+1して1はじまりにする
# 釣り人の数
K = fish_num_climate_3['human'].nunique()
# 辞書にまとめる ※Y:pd.Series, X:pd.DataFrame
data_dict_design = dict(N=N, P=P, K=K, Y=Y['fish_num'], X=X, human=human)【実行結果】なし
④ モデルのコンパイル
モデルのコンパイルを実行します。
%%time
# モデルのコンパイル
# stanプログラムファイルのパス指定
stan_file = '4-2-1-glmm-pois-design-matrix.stan'
current_dir = os.path.abspath(os.getcwd())
stan_path = os.path.join(current_dir, 'stan', stan_file)
# モデルオブジェクトの作成（exeの作成）
model_glmm = CmdStanModel(stan_file=stan_path) # stanのpathを設定【実行結果】
MCMC の準備が整いました！
⑤ MCMC の実行
MCMCを実行しましょう。
%%time
# p.257 MCMCの実行
fit_glmm = model_glmm.sample(
data=data_dict_design, # 対象データ
seed=1, # 乱数の種
sig_figs=18, # 出力CSV等に適用する数値精度
)【実行結果】
⑥ 収束確認
収束の確認をします。
診断メソッド diagnose を利用します。
# 事後分布の診断
print(fit_glmm.diagnose())【実行結果】（１行目のファイルパスは記載省略）
問題は検出されませんでした（no problems detected.）。
CmdStanPy の 出力結果 fit から MCMCサンプルの要約統計量を把握しましょう。
# p.257 結果の表示
fit_glmm.summary(percentiles=[2.5, 50, 97.5]).iloc[[0, 1, 2, 3, 14], :].round(2)【実行結果】
$${\widehat{R}}$$（R_hat）、有効サンプル数（N_Eff）に問題はなさそうです。
b[1] は切片、b[2] は晴れ、b[3] は気温の係数に対応しています。
fit を arviz の idata に変換します。
# arvizのidataに変換
idata_glmm = az.from_cmdstanpy(posterior=fit_glmm, log_likelihood='lp__')
idata_glmm【実行結果】
収束の追加確認をします。
arviz の rhat 関数を用いて、$${\widehat{R} > 1.01}$$ のパラメータがないことを数値で確かめます。
# 収束の確認 r_hat>1.01の確認
# 設定
idata_in = idata_glmm # idata名
threshold = 1.01 # しきい値
# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())【実行結果】
4B の右側の数値が $${\widehat{R} > 1.01}$$ のパラメータの個数です。
すべて $${0}$$ ですので、$${\widehat{R} \leq 1.01}$$ と言えます。
トレースプロットを描画します。
# トレースプロットの描画
az.plot_trace(idata_glmm, var_names=['b', 'sigma_r'], compact=False,
backend_kwargs={'tight_layout': True});【実行結果】
左側のチャートの４本の Chain はほぼ重なっており、かつ１峰です。
右側のチャートがゲジゲジしています。
以上のチェックに基づいて、収束していると考えましょう。
⑥ 推定されたモデルの解釈
線形予測子の数式表現を再掲します。
$$
\log(\lambda_i) = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + r_k \\
$$
両辺の指数 $${\exp}$$ をとると次のようになります。
$$
\begin{align*}
\lambda_i &= \exp(\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + r_k) \\
&= \exp(\beta_0) \times \exp(\beta_1 x_{i1}) \times \exp(\beta_2 x_{i2}) \times \exp(r_k)
\end{align*}
$$
係数（および変数の値）の変動によって、平均パラメータ $${\lambda_i}$$ は指数 $${\exp(\cdot)}$$ 倍変動します。
これを踏まえて、GLMM モデルの事後分布要約統計量を確認しましょう。
前回記事の Gemini 謹製分析テンプレを利用して読み解きを進めます。
1. 気温の影響は「確実」にある
$${\texttt{b[3]}}$$ の 95%信用区間が 0.08 〜 0.11 となっており、0 を含まずにプラス側に振り切っています。
気温が上がると釣果が伸びるという傾向は、この過分散なデータの中でも統計的にしっかりと検出できています。
2. ランダム効果の標準偏差 $${\texttt{sigma\_r}}$$ の存在感
ここが今回の肝です。ランダム切片の標準偏差の平均が 0.65 です。
対数リンク関数を使っているので、この数値は「個体ごとのゆらぎによって、釣果が平均的に（1 標準偏差的に）$${\exp(0.65) \approx 1.9}$$ 倍程度、あるいは $${\exp(-0.65) \approx 0.52}$$ 程度まで変動する」という、ノイズの大きさを物語っています。
3. 天気の効果 $${\texttt{b[2]}}$$
$${\texttt{b[2]}}$$ が -0.52 です。指数をとると $${\exp(-0.52) \approx 0.59}$$。
つまり、「晴れ」は「曇り」に比べて、気温などの条件が同じでも釣果が半分近くまで落ち込む傾向がある、と読み解けます。
⑦ 釣り人の影響力（ランダム切片）の確認
テキスト p.257 の ranef 関数と同等のことを CmdStanPy の fit_glmm に対する summary メソッドで実行します。95%信用区間を表示します。
# p.257 釣り人の影響の大きさ（ランダム切片の要約統計量）
fit_glmm.summary(percentiles=[2.5, 50, 97.5]).iloc[4:14, :].round(3)【実行結果】
A ～ H の釣り人のランダム切片の推定値です。
指数変換前ですので、$${\log(\lambda_i)}$$ への加減算の効果です。
テキストのとおり、Aさんは釣りが得意、D・Fさんは釣りが苦手なのかもしれません。
ちなみに、変換した idata のランダム切片の要約統計量を az.summary 関数で表示しましょう。95%HDI 区間を表示します。
# 参考：arvizのサマリー
az.summary(idata_glmm, var_names=['r'])【実行結果】
ランダム効果を可視化しましょう。
idata のランダム切片をフォレストプロットで表現します。
arviz の plot_forest 関数を利用します。
指数変換を行って、スケールを「倍率」にしてみましょう。
# フォレストプロット（expで「倍率」スケールに変換後）
# 事後分布のMCMCサンプルを指数変換 (exp) する
# これにより、単位が「対数」から「倍率」に変わります
ds_ratio = idata_glmm.posterior[['r']].copy()
ds_ratio['r'] = np.exp(ds_ratio['r'])
# フォレストプロットの描画
az.plot_forest(ds_ratio, var_names=['r'], combined=True, hdi_prob=0.95)
plt.axvline(1, color='tab:red');【実行結果】
値が１倍より大きい釣り人 A（0）、C（2）さんは正の倍率効果がありそうです。
一方で１倍より小さい釣り人 D（3）、F（5）さんは負の倍率効果がありそうです。
その他の釣り人は 95%HDI 区間が１倍を含んでいるので、正・負の判断をしかねます。
⑧ 釣り人別・天気別の釣獲尾数の可視化
テキストの図 4.2.1 に相当します。
釣り人別にチャートを分割して、天気別平均獲釣尾数の平均値と 95%HDI 区間を描画します。
まず、平均獲釣尾数の平均値と 95%HDI 区間を算出します。
Stan で生成したMCMCサンプルを使い、Stan の外で、GLMMの数式に沿って $${\bm \lambda}$$ を計算します。
# 平均釣獲尾数のMCMCサンプルの生成 ※Stanの外で算出
# MCMCサンプルから取り出し
# 係数β0, β1, β2 各shape=(4000,)
intercept, sunny, temper = az.extract(idata_glmm.posterior).b.to_numpy()
# ランダム切片r shape=(10, 4000)
rs = az.extract(idata_glmm.posterior).r.to_numpy()
# x軸の値（気温）の設定
x_val = np.linspace(fish_num_climate_3['temperature'].min(),
fish_num_climate_3['temperature'].max(), 100)
# 平均釣獲尾数のMCMCサンプルの生成
# ※numpy多次元配列：shape=(10, 2, 100, 4000)⇒(釣り人, 晴れ, 気温, MCMCサンプル)
pred_samples = np.exp(
intercept # 切片 shape=(4000)
+ np.outer([0, 1], sunny)[None, :, None, :] # 晴れ shape=(2,4000)
+ np.outer(x_val, temper)[None, None, :, :] # 気温 shape=(100, 4000)
+ rs[:, None, None, :] # ランダム切片 shape=(10, 4000)
)
print('pred_samples.shape:', pred_samples.shape)【実行結果】
４軸の numpy 配列にしています。
第１軸：釣り人 10
第２軸：晴れダミー 2（曇りと晴れ）
第３軸：気温 100 個
第４軸：MCMCサンプル 4000 個
描画します。
seaborn の relplot の axes を取得して、その上に平均値・95%HDI 区間を重ね描きするのがこの記事のクライマックスです！（当社比）
arviz の plot_hdi 関数で 95%HDI 区間を描画します。
# p.259 図4.2.1 釣り人別の回帰曲線
# 描画ヘルパー関数
def plot_helper_function(pred_samples):
# 釣り人別散布図の描画
g = sns.relplot(
data=fish_num_climate_3, x='temperature', y='fish_num', # データ, x軸, y軸
col='human', col_wrap=3, # 列の変数, 列の数
hue='weather', palette=colors, # hueの変数、色
kind='scatter', s=70, height=3, aspect=1.2
)
# 釣り人iごとの平均釣獲尾数の平均値と95% HDI区間の描画
# sns.relplotのaxesを取り出して、釣り人ごとに描画を繰り返し処理
for i, ax in enumerate(g.axes.flat):
# 天気jごとに描画を繰り返し処理
for j, color in enumerate(colors.values()):
# 対象の釣り人iと天気jをMCMCサンプルから取得
pred_sel = pred_samples[i, j, :, :]
# 平均釣獲尾数の平均値の描画
ax.plot(x_val, pred_sel.mean(axis=1), color=color)
# 平均釣獲尾数の95% HDI区間の塗りつぶし描画 ※arvizのplot_hdi利用
az.plot_hdi(x_val, pred_sel.T, hdi_prob=0.95, color=color,
fill_kwargs={'alpha': 0.2}, ax=ax)
# 修飾
ax.set(title=f'釣り人 {humans[i]}', xlabel='', ylabel='')
# チャート全体のx,y軸ラベルの表示
g.figure.supxlabel('気温 [℃]', fontsize=14)
g.figure.supylabel('釣獲尾数', fontsize=14);
# 描画処理
plot_helper_function(pred_samples)【実行結果】
気温が高くなるにつれて釣獲尾数が増える傾向と、曇りの方が晴れよりも釣獲尾数が多い傾向が可視化されています。
先ほど確認した「気温と晴れダミーのパラメータ推定値」を釣り人共通（ランダム効果を除く線形予測子）で適用して、予測値を算出しているので、当然と言えば当然です。
あとは、釣り人の効果としてランダム切片を加味することで、各釣り人の回帰曲線が大きくなったり小さくなったりしています。
🔵🔵🔵
アディショナルタイム：GPBoost 実装例
非ベイズの GLMM をやってみましょう。
前回記事で利用した GPBoost で今回のランダム切片モデルを実装します。
① 追加インポート
# ライブラリの追加インポート
import gpboost as gpb
gpb.__version__【実行結果】
今回利用する GPBoost のバージョンは 1.6.1 です（最新版では無いです）。
② GLMM のモデリング
ランダム切片の変数 human を 引数 group_data に設定します。
ランダム切片と説明変数に pandas の DataFrame、目的変数に pandas の Series を渡すと、分析結果に変数名が表示されるのでおすすめです。
このモデリングでは先ほど作成したデザイン行列を再利用しています。
# ポアソンGLMMモデルによる最尤推定：ラプラス近似を用いた（制限付き）最尤法
# 1. データセットの作成
# ランダム効果（グループ）を作成 ※pandas DataFrame
group_data = fish_num_climate_3[['human']]
# 2. GLMMの実行
# モデルの定義：ランダム効果はグループ（ID）を指定
gp_model = gpb.GPModel(group_data=group_data, likelihood='poisson')
# 学習 ※説明変数には先ほど作成したデザイン行列を利用
gp_model.fit(y=Y_dm.iloc[:, 0], X=X_dm, params={'std_dev': True})
# 結果サマリーの表示 ※ランダム効果のCovarianceは分散（標準偏差でない）
print(gp_model.summary())【実行結果】
③ 予測の実行
GPBoost のモデルで予測を行います。
予測用の説明変数データ（デザイン行列等）を作成し、モデル gp_model に対して predict メソッドを適用して予測し、平均値と95%信頼区間を算出してデータフレーム化します。
# 1. 予測用データの作成
# 予測条件の作成 (10人 × 2天気 × 100気温 = 2000行)
x_val = np.linspace(*np.sort(fish_num_climate_3['temperature'])[[0, -1]], 100)
# 定数項 × 釣り人 × 天気 × 気温 の全組み合わせ（直積）の作成
prod = list(product(range(1,2), humans, range(2), x_val))
# データフレーム化
pred_gp_df = pd.DataFrame(
prod, columns=['Intercept', 'human', 'weather_idx', 'temperature']
)
# 2. GPModel.predict による潜在変数の予測
pred_gp = gp_model.predict(
X_pred=pred_gp_df[['Intercept', 'weather_idx', 'temperature']],
group_data_pred=pred_gp_df['human'],
predict_response=False, # 線形予測子スケール（log(λ)）を取得
predict_var=True # 分散を取得
)
# 'mu' と 'var' を取得
mu_latent = pred_gp['mu']
var_latent = pred_gp['var']
# 3. 期待値と95%信頼区間の算出 (指数変換)
pred_gp_df['mean_fish'] = np.exp(mu_latent)
pred_gp_df['lower'] = np.exp(mu_latent - 1.96 * np.sqrt(var_latent))
pred_gp_df['upper'] = np.exp(mu_latent + 1.96 * np.sqrt(var_latent))
pred_gp_df['weather'] = pred_gp_df['weather_idx'].map({0: 'cloudy', 1: 'sunny'})
pred_gp_df【実行結果】
Intercept ～ temperature までが予測に用いた説明変数、mean_fish が予測平均値、lower と upper が予測平均の 95% 信頼区間です。
④ 回帰曲線の可視化
テキスト 図 4.2.1 のような釣り人別・天気別の回帰曲線を描画します。
# 釣り人別の回帰曲線の描画（図4.2.1 ライク）
# 釣り人別散布図の描画
g = sns.relplot(
data=fish_num_climate_3, x='temperature', y='fish_num', # データ, x軸, y軸
col='human', col_wrap=3, # 列の変数, 列の数
hue='weather', palette=colors, # hueの変数、色
kind='scatter', s=70, height=3, aspect=1.2
)
# 釣り人別の平均釣獲尾数の平均値と95%信頼区間の描画
# sns.relplotのaxesを取り出して、釣り人ごとに描画を繰り返し処理
for human, ax in g.axes_dict.items():
# 天気ごとに描画を繰り返し処理
for weather in colors.keys():
# 対象の釣り人と天気をデータフレームから取得
query = (pred_gp_df['human']==human) & (pred_gp_df['weather']==weather)
subset = pred_gp_df[query]
# 平均釣獲尾数の平均値の描画
ax.plot(subset['temperature'], subset['mean_fish'], color=colors[weather])
# 平均釣獲尾数の95%信頼区間の塗りつぶし描画
ax.fill_between(subset['temperature'], subset['lower'], subset['upper'],
color=colors[weather], alpha=0.2)
# 修飾
ax.set(title=f'釣り人 {human}', xlabel='', ylabel='')
# チャート全体のx,y軸ラベルの表示
g.figure.supxlabel('気温 [℃]', fontsize=14)
g.figure.supylabel('釣獲尾数', fontsize=14)
g._legend.set_title('天気');【実行結果】
ベイズモデルとよく似た描画結果になりました。
🔵🔵🔵
今回の記事は以上です。
楽しかったですね！
シリーズの記事
次の記事
前の記事
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の教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
いいなと思ったら応援しよう！
応援ありがとうございます。これからもがんばって記事を作成します！
