「RとStanで始めるベイズ統計モデリングによるデータ分析入門」「実践編第8章 ポアソン回帰モデル」を対象に,公開されているR,Stanのコードをpython,pystanのコードへと書き直した一例です。Stanの代わりにpystanを利用しています。
本ページでは公開されていない書籍の内容については一切触れません。理論や詳しい説明は書籍を参照してください。
なお,こちらで紹介しているコードには誤りが含まれる可能性があります。内容やコードについてお気づきの点等ございましたら,ご指摘いただけると幸いです。
(本ページにて紹介しているコードはgithubにて公開しています。)
DataAnalOji.hatena.sample/python_samples/stan/3-8-ポアソン回帰モデル.ipynb at master · Data-Anal-Ojisan/DataAnalOji.hatena.sample
samples for my own blog. Contribute to Data-Anal-Ojisan/DataAnalOji.hatena.sample development by creating an account on GitHub.
分析の準備
パッケージの読み込み
import pystan
import numpy as np
import pandas as pd
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('ggplot')
plt.rcParams['font.family'] = 'Meiryo'
import seaborn as sns
データの読み込みと図示
分析対象のデータ
fish_num_climate = pd.read_csv('3-8-1-fish-num-1.csv')
fish_num_climate.head(n=3)
データの要約
fish_num_climate.describe(include='all')
図示
# グラフ描画領域の作成
plt.figure(figsize=(10, 5))
# 散布図の描画
sns.scatterplot(x='temperature', y='fish_num', data=fish_num_climate, hue='weather')
# タイトルの変更
plt.title('釣獲尾数と気温・天気の関係')
# グラフの描画
plt.show()
brmsによる分散分析モデルの推定
pythonにbrmsの代用可能なパッケージが知っている範囲内で存在しないため,brmsを利用したコードは省略しています。
補足:brmsを用いない実装の方法
参考:デザイン行列の作成
# 応答変数の削除とダミー変数化処理
design_mat = pd.get_dummies(fish_num_climate.drop(['fish_num'], axis=1),
drop_first=True)
# 列名を書籍準拠に
design_mat = design_mat[['weather_sunny', 'temperature']]
# (Intercept)列追加
design_mat.insert(0, '(Intercept)', 1)
display(design_mat)
参考:データの作成
data_list_1 = dict(N=len(fish_num_climate),
fish_num=fish_num_climate['fish_num'],
temp=fish_num_climate['temperature'],
sunny=design_mat['weather_sunny'])
display(data_list_1)
{'N': 100,
'fish_num': 0 0
1 2
2 5
3 1
4 3
..
95 0
96 0
97 0
98 0
99 0
Name: fish_num, Length: 100, dtype: int64,
'temp': 0 5.5
1 21.1
2 17.2
3 5.0
4 28.3
...
95 5.9
96 12.8
97 2.8
98 3.5
99 13.2
Name: temperature, Length: 100, dtype: float64,
'sunny': 0 0
1 0
2 0
3 0
4 0
..
95 1
96 1
97 1
98 1
99 1
Name: weather_sunny, Length: 100, dtype: uint8}
参考:自分で変換処理を入れる
# stanコードの記述
stan_code = '''
data {
int N; // サンプルサイズ
int fish_num[N]; // 釣獲尾数
vector[N] temp; // 気温データ
vector[N] sunny; // 晴れダミー変数
}
parameters {
real Intercept; // 切片
real b_temp; // 係数(気温)
real b_sunny; // 係数(晴れの影響)
}
model {
vector[N] lambda = exp(Intercept + b_temp*temp + b_sunny*sunny);
fish_num ~ poisson(lambda);
}
'''
# モデルのコンパイル
stan_model = pystan.StanModel(model_code=stan_code)
# サンプリング
glm_pois_stan_exp = stan_model.sampling(data=data_list_1, seed=2, n_jobs=1)
結果の表示
print(glm_pois_stan_exp.stansummary(probs=[0.025, 0.5, 0.975]))
Inference for Stan model: anon_model_10b93b39bc827a53de74eb5431c131b8.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
Intercept -0.79 6.0e-3 0.24 -1.3 -0.79 -0.35 1628 1.0
b_temp 0.08 2.4e-4 0.01 0.06 0.08 0.1 1682 1.0
b_sunny -0.59 3.7e-3 0.17 -0.93 -0.59 -0.26 2089 1.0
lp__ -37.73 0.03 1.22 -40.95 -37.41 -36.32 1412 1.0
Samples were drawn using NUTS at Tue Sep 1 19:30:40 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
参考:poisson_log関数を使用
# stanコードの記述
stan_code = '''
data {
int N; // サンプルサイズ
int fish_num[N]; // 釣獲尾数
vector[N] temp; // 気温データ
vector[N] sunny; // 晴れダミー変数
}
parameters {
real Intercept; // 切片
real b_temp; // 係数(気温)
real b_sunny; // 係数(晴れの影響)
}
model {
vector[N] lambda = Intercept + b_temp*temp + b_sunny*sunny;
fish_num ~ poisson_log(lambda);
}
generated quantities {
vector[N] lambda_sunny;
vector[N] lambda_cloudy;
vector[N] fish_num_sunny;
vector[N] fish_num_cloudy;
for(i in 1:N){
lambda_sunny[i] = Intercept + b_temp*temp[i] + b_sunny*1;
lambda_cloudy[i] = Intercept + b_temp*temp[i] + b_sunny*0;
fish_num_sunny[i] = poisson_log_rng(lambda_sunny[i]);
fish_num_cloudy[i] = poisson_log_rng(lambda_cloudy[i]);
}
}
'''
# モデルのコンパイル
stan_model = pystan.StanModel(model_code=stan_code)
# サンプリング
glm_pois_stan = stan_model.sampling(data=data_list_1, n_jobs=1)
参考:結果の表示
Inference for Stan model: anon_model_6d78a2071bf51cbe1e42c5c01664ea40.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
Intercept -0.79 6.1e-3 0.24 -1.28 -0.79 -0.32 1553 1.0
b_temp 0.08 2.4e-4 0.01 0.06 0.08 0.1 1755 1.0
b_sunny -0.59 3.9e-3 0.17 -0.93 -0.59 -0.26 1939 1.0
lambda_sunny[1] -0.93 4.4e-3 0.21 -1.35 -0.92 -0.54 2248 1.0
・・・省略・・・
lambda_sunny[100] -0.29 3.1e-3 0.16 -0.61 -0.29 0.01 2596 1.0
lambda_cloudy[1] -0.34 4.9e-3 0.19 -0.72 -0.34 0.03 1557 1.0
・・・省略・・・
lambda_cloudy[100] 0.3 3.2e-3 0.13 0.04 0.3 0.55 1690 1.0
fish_num_sunny[1] 0.41 0.01 0.65 0.0 0.0 2.0 3641 1.0
・・・省略・・・
fish_num_sunny[100] 0.74 0.01 0.86 0.0 1.0 3.0 3788 1.0
fish_num_cloudy[1] 0.73 0.01 0.86 0.0 1.0 3.0 4072 1.0
・・・省略・・・
fish_num_cloudy[100] 1.35 0.02 1.17 0.0 1.0 4.0 3936 1.0
lp__ -37.78 0.03 1.25 -40.99 -37.46 -36.31 1352 1.0
Samples were drawn using NUTS at Tue Sep 1 19:31:19 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
ポアソン回帰の回帰曲線
95%ベイズ信頼区間付きのグラフ
MCMCサンプルの抽出
mcmc_sample = glm_pois_stan.extract()
事後予測分布を天気別にデータフレームへまとめる
result_lambda_cloudy = pd.DataFrame(np.zeros([len(fish_num_climate), 3]))
result_lambda_sunny = pd.DataFrame(np.zeros([len(fish_num_climate), 3]))
# 2.5・50・97.5パーセンタイルの値をデータフレームに代入
for i in range(len(fish_num_climate)):
result_lambda_cloudy.iloc[i,:] = np.percentile(mcmc_sample['lambda_cloudy'][:,i], q=[2.5, 50, 97.5])
result_lambda_sunny.iloc[i,:] = np.percentile(mcmc_sample['lambda_sunny'][:,i], q=[2.5, 50, 97.5])
# データフレームの列名を変更
result_lambda_cloudy.columns = ["2.5%", "50%", "97.5%"]
result_lambda_sunny.columns = ["2.5%", "50%", "97.5%"]
# 温度情報をデータフレームに追加
result_lambda_cloudy['temperature'] = fish_num_climate['temperature']
result_lambda_sunny['temperature'] = fish_num_climate['temperature']
# 温度情報で並び順をソート
result_lambda_cloudy = result_lambda_cloudy.sort_values('temperature')
result_lambda_sunny = result_lambda_sunny.sort_values('temperature')
回帰直線の描画
# グラフ描画領域の作成
plt.figure(figsize=(10, 5))
# 散布図の描画
sns.scatterplot(x='temperature',
y='fish_num',
data=fish_num_climate,
hue='weather',
hue_order=['cloudy', 'sunny'],
palette=['red', 'blue'])
# 曇りの日のビールの売り上げの平均値を描画
plt.plot(result_lambda_cloudy['temperature'],
np.exp(result_lambda_cloudy["50%"]),
color='red')
# 95%ベイズ信頼区間を描画
plt.fill_between(x=result_lambda_cloudy['temperature'],
y1=np.exp(result_lambda_cloudy["2.5%"]),
y2=np.exp(result_lambda_cloudy["97.5%"]),
color='red',
alpha=0.3)
# 晴れの日のビールの売り上げの平均値を描画
plt.plot(result_lambda_sunny['temperature'],
np.exp(result_lambda_sunny["50%"]),
color='blue')
# 95%ベイズ信頼区間を描画
plt.fill_between(x=result_lambda_sunny['temperature'],
y1=np.exp(result_lambda_sunny["2.5%"]),
y2=np.exp(result_lambda_sunny["97.5%"]),
color='blue',
alpha=0.3)
# グラフタイトル追加
plt.title('ポアソン回帰曲線:信用区間')
# グラフの描画
plt.show()
99%ベイズ予測区間付きのグラフ
事後予測分布を天気別にデータフレームへまとめる
result_fish_num_cloudy = pd.DataFrame(np.zeros([len(fish_num_climate), 3]))
result_fish_num_sunny = pd.DataFrame(np.zeros([len(fish_num_climate), 3]))
# 0.5・50・99.5パーセンタイルの値をデータフレームに代入
for i in range(len(fish_num_climate)):
result_fish_num_cloudy.iloc[i,:] = np.percentile(mcmc_sample['fish_num_cloudy'][:,i], q=[0.5, 50, 99.5])
result_fish_num_sunny.iloc[i,:] = np.percentile(mcmc_sample['fish_num_sunny'][:,i], q=[0.5, 50, 99.5])
# データフレームの列名を変更
result_fish_num_cloudy.columns = ["0.5%", "50%", "99.5%"]
result_fish_num_sunny.columns = ["0.5%", "50%", "99.5%"]
# 温度情報をデータフレームに追加
result_fish_num_cloudy['temperature'] = fish_num_climate['temperature']
result_fish_num_sunny['temperature'] = fish_num_climate['temperature']
# 温度情報で並び順をソート
result_fish_num_cloudy = result_fish_num_cloudy.sort_values('temperature')
result_fish_num_sunny = result_fish_num_sunny.sort_values('temperature')
回帰直線の描画
# グラフ描画領域の作成
plt.figure(figsize=(10, 5))
# 散布図の描画
sns.scatterplot(x='temperature',
y='fish_num',
data=fish_num_climate,
hue='weather',
hue_order=['cloudy', 'sunny'],
palette=['red', 'blue'])
# 曇りの日のビールの売り上げの平均値を描画
plt.plot(result_fish_num_cloudy['temperature'],
result_fish_num_cloudy["50%"],
color='red')
# 95%ベイズ信頼区間を描画
plt.fill_between(x=result_fish_num_cloudy['temperature'],
y1=result_fish_num_cloudy["0.5%"],
y2=result_fish_num_cloudy["99.5%"],
color='red',
alpha=0.3)
# 晴れの日のビールの売り上げの平均値を描画
plt.plot(result_fish_num_sunny['temperature'],
result_fish_num_sunny["50%"],
color='blue')
# 95%ベイズ信頼区間を描画
plt.fill_between(x=result_fish_num_sunny['temperature'],
y1=result_fish_num_sunny["0.5%"],
y2=result_fish_num_sunny["99.5%"],
color='blue',
alpha=0.3)
# グラフタイトル追加
plt.title('ポアソン回帰曲線:予測区間')
# グラフの描画
plt.show()
補足:デザイン行列を使ったモデルの推定
参考:Stanに渡すデータ
data_list_2 = dict(N=len(fish_num_climate),
K=design_mat.shape[1],
Y=fish_num_climate['fish_num'],
X=design_mat)
display(data_list_2)
{'N': 100,
'K': 3,
'Y': 0 0
1 2
2 5
3 1
4 3
..
95 0
96 0
97 0
98 0
99 0
Name: fish_num, Length: 100, dtype: int64,
'X': (Intercept) weather_sunny temperature
0 1 0 5.5
1 1 0 21.1
2 1 0 17.2
3 1 0 5.0
4 1 0 28.3
.. ... ... ...
95 1 1 5.9
96 1 1 12.8
97 1 1 2.8
98 1 1 3.5
99 1 1 13.2
[100 rows x 3 columns]}
参考:MCMCの実行
# stanコードの記述
stan_code = '''
data {
int N; // サンプルサイズ
int K; // デザイン行列の列数(説明変数の数+1)
int Y[N]; // 応答変数(整数型)
matrix[N, K] X; // デザイン行列
}
parameters {
vector[K] b; // 切片を含む係数ベクトル
}
model {
vector[N] lambda = X * b;
Y ~ poisson_log(lambda);
}
'''
# モデルのコンパイル
stan_model = pystan.StanModel(model_code=stan_code)
# サンプリング
glm_pois_stan_design_mat = stan_model.sampling(data=data_list_2, seed=1, n_jobs=1)
参考:結果の表示
Inference for Stan model: anon_model_8fb9bb9e4fdb448c421e95a5e990a395.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
b[1] -0.81 6.4e-3 0.24 -1.3 -0.8 -0.36 1426 1.0
b[2] -0.59 3.5e-3 0.17 -0.93 -0.59 -0.27 2284 1.0
b[3] 0.08 2.7e-4 0.01 0.06 0.08 0.1 1433 1.0
lp__ -37.69 0.03 1.22 -40.94 -37.38 -36.31 1291 1.0
Samples were drawn using NUTS at Tue Sep 1 19:32:01 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
pystanの詳細については公式ページを参照してください。
PyStan — pystan 3.10.0 documentation
コメント