実践Ｄａｔａ Ｓｃｉｅｎｃｅシリーズ ＲとＳｔａｎではじめる ベイズ統計モデリングによるデータ分析入門 (ＫＳ情報科学専門書)
「RとStanで始めるベイズ統計モデリングによるデータ分析入門 実践編第7章 正規線形モデル」を対象に，公開されているR，Stanのコードをpython，pystanのコードへと書き直した一例です。Stanの代わりにpystanを利用しています。
本ページでは公開されていない書籍の内容については一切触れません。理論や詳しい説明は書籍を参照してください。
なお，こちらで紹介しているコードには誤りが含まれる可能性があります。内容やコードについてお気づきの点等ございましたら，ご指摘いただけると幸いです。
（本ページにて紹介しているコードはgithubにて公開しています。）
DataAnalOji.hatena.sample/python_samples/stan/3-7-正規線形モデル.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.
github.com
分析の準備
パッケージの読み込み
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データの読み込みと図示
分析対象のデータ
sales_climate = pd.read_csv('3-7-1-beer-sales-4.csv')
sales_climate.head(n=3)
データの要約
sales_climate.describe(include='all')
図示
plt.figure(figsize=(10, 5))
sns.scatterplot(x='temperature',
y='sales',
data=sales_climate,
hue='weather',
hue_order=['cloudy', 'rainy', 'sunny'],
palette=['red', 'green', 'blue'])
plt.title('ビールの売り上げと気温・天気の関係')
plt.show()
brmsによる正規線形モデルの推定
pythonにbrmsの代用可能なパッケージが知っている範囲内で存在しないため，brmsを利用したコードは省略しています。
補足：正規線形モデルのデザイン行列
デザイン行列の作成
design_mat = pd.get_dummies(sales_climate.drop('sales',axis=1), drop_first=True) # ダミー変数処理
design_mat.insert(0, '(Intercept)', 1) # (Intercept)列追加
design_mat = design_mat[['(Intercept)', 'weather_rainy', 'weather_sunny', 'temperature']] # 列名並び替え
design_mat.head(n=3)
補足：brmsを使わない正規線形モデルの推定
dictにまとめる
data_list_design = dict(N=len(design_mat),
K=design_mat.shape[1],
Y=sales_climate['sales'],
X=design_mat)MCMCの実行
# stanコードの記述
stan_code = '''
data {
int N; // サンプルサイズ
int K; // デザイン行列の列数(説明変数の数＋１)
vector[N] Y; // 応答変数
matrix[N, K] X; // デザイン行列
}
parameters {
vector[K] b; // 切片を含む係数ベクトル
real<lower=0> sigma; // データのばらつきを表す標準偏差
}
model {
vector[N] mu = X * b;
Y ~ normal(mu, sigma);
}
'''
# モデルのコンパイル
stan_model = pystan.StanModel(model_code=stan_code)
# サンプリング
lm = stan_model.sampling(data=data_list_design, seed=1, n_jobs=1)MCMCの結果の確認
print(lm)Inference for Stan model: anon_model_97e723c47136639268ab4817ab64884e.
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% 25% 50% 75% 97.5% n_eff Rhat
b[1] 20.32 0.12 5.1 10.46 16.97 20.2 23.7 30.69 1788 1.0
b[2] -3.55 0.07 3.21 -9.85 -5.74 -3.49 -1.42 2.58 2258 1.0
b[3] 29.39 0.07 3.2 23.28 27.25 29.34 31.49 35.88 2257 1.0
b[4] 2.54 4.9e-3 0.22 2.09 2.4 2.54 2.69 2.98 2024 1.0
sigma 16.04 0.02 0.93 14.32 15.39 16.01 16.62 17.97 3069 1.0
lp__ -487.9 0.04 1.58 -491.8 -488.7 -487.5 -486.7 -485.8 1957 1.0
Samples were drawn using NUTS at Tue Jul 28 23:20:47 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).補足：回帰直線
dictにまとめる
10~30℃の温度に対応する天気別のビールの売り上げに関する事後予測分布を得て，回帰直線を描画するために10~30℃の温度データを入力データに加える。
data_list_design = dict(N=len(design_mat),
K=design_mat.shape[1],
Y=sales_climate['sales'],
X=design_mat,
N_temperature = len(np.arange(10,31,1)), # 10-30℃までのデータ長
temperature = np.arange(10,31,1)) # 10-30℃までの温度MCMCの実行
# stanコードの記述
stan_code = '''
data {
int N; // サンプルサイズ
int K; // デザイン行列の列数(説明変数の数＋１)
vector[N] Y; // 応答変数
matrix[N, K] X; // デザイン行列
int N_temperature; // 温度データの長さ
vector[N_temperature] temperature; // 温度データ
}
parameters {
vector[K] b; // 切片を含む係数ベクトル
real<lower=0> sigma; // データのばらつきを表す標準偏差
}
model {
vector[N] mu = X * b;
Y ~ normal(mu, sigma);
}
generated quantities {
vector[N_temperature] mu_pred_cloudy; // 曇りの日の平均売り上げ
vector[N_temperature] mu_pred_rainy; // 雨の日の平均売り上げ
vector[N_temperature] mu_pred_sunny; // 晴れの日の平均売り上げ
for(i in 1:N_temperature){
mu_pred_cloudy[i] = (b[1] + b[2]*0 + b[3]*0 + b[4]*temperature[i]);
mu_pred_rainy[i] = (b[1] + b[2]*1 + b[3]*0 + b[4]*temperature[i]);
mu_pred_sunny[i] = (b[1] + b[2]*0 + b[3]*1 + b[4]*temperature[i]);
}
}
'''
# モデルのコンパイル
stan_model = pystan.StanModel(model_code=stan_code)
# サンプリング
lm = stan_model.sampling(data=data_list_design, seed=1, n_jobs=1)MCMCの結果の確認
print(lm)Inference for Stan model: anon_model_28f6803b7c29c57caa237a3515c03927.
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% 25% 50% 75% 97.5% n_eff Rhat
b[1] 20.32 0.12 5.1 10.46 16.97 20.2 23.7 30.69 1788 1.0
b[2] -3.55 0.07 3.21 -9.85 -5.74 -3.49 -1.42 2.58 2258 1.0
b[3] 29.39 0.07 3.2 23.28 27.25 29.34 31.49 35.88 2257 1.0
b[4] 2.54 4.9e-3 0.22 2.09 2.4 2.54 2.69 2.98 2024 1.0
sigma 16.04 0.02 0.93 14.32 15.39 16.01 16.62 17.97 3069 1.0
mu_pred_cloudy[1] 45.76 0.08 3.26 39.5 43.56 45.76 47.95 52.23 1802 1.0
省略・・・
mu_pred_cloudy[21] 96.65 0.06 3.07 90.79 94.58 96.65 98.73 102.7 3000 1.0
mu_pred_rainy[1] 42.22 0.06 3.23 35.85 40.02 42.21 44.36 48.48 2661 1.0
省略・・・
mu_pred_rainy[21] 93.1 0.06 3.13 86.76 91.05 93.11 95.17 99.27 2565 1.0
mu_pred_sunny[1] 75.16 0.06 3.1 69.11 73.1 75.06 77.22 81.34 2983 1.0
省略・・・
mu_pred_sunny[21] 126.04 0.07 3.21 119.75 123.88 126.07 128.2 132.47 2299 1.0
lp__ -487.9 0.04 1.58 -491.8 -488.7 -487.5 -486.7 -485.8 1957 1.0
Samples were drawn using NUTS at Tue Jul 28 23:21:26 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).
MCMCサンプルの抽出
mcmc_sample = lm.extract()天気別ビール売り上げ平均値の事後予測分布をデータフレームにまとめる
result_cloudy = pd.DataFrame(np.zeros([len(np.arange(10,31,1)),3]))
result_rainy = pd.DataFrame(np.zeros([len(np.arange(10,31,1)),3]))
result_sunny = pd.DataFrame(np.zeros([len(np.arange(10,31,1)),3]))
# 2.5・50・97.5パーセンタイルの値をデータフレームに代入
for i in range(21):
result_cloudy.iloc[i,:] = np.percentile(mcmc_sample['mu_pred_cloudy'][:,i], q=[2.5, 50, 97.5])
result_rainy.iloc[i,:] = np.percentile(mcmc_sample['mu_pred_rainy'][:,i], q=[2.5, 50, 97.5])
result_sunny.iloc[i,:] = np.percentile(mcmc_sample['mu_pred_sunny'][:,i], q=[2.5, 50, 97.5])
# データフレームの列名を変更
result_cloudy.columns = ["2.5%", "50%", "97.5%"]
result_rainy.columns = ["2.5%", "50%", "97.5%"]
result_sunny.columns = ["2.5%", "50%", "97.5%"]
# 温度情報をデータフレームに追加
result_cloudy['temperature'] = np.arange(10,31,1)
result_rainy['temperature'] = np.arange(10,31,1)
result_sunny['temperature'] = np.arange(10,31,1)回帰直線の描画
# グラフ描画領域の作成
plt.figure(figsize=(10, 5))
# 散布図の描画
sns.scatterplot(x='temperature',
y='sales',
data=sales_climate,
hue='weather',
hue_order=['cloudy', 'rainy', 'sunny'],
palette=['red', 'green', 'blue'])
# 曇りの日のビールの売り上げの平均値を描画
plt.plot(result_cloudy['temperature'], result_cloudy["50%"], color='red')
# 95%ベイズ信頼区間を描画
plt.fill_between(x=result_cloudy['temperature'],
y1=result_cloudy["2.5%"],
y2=result_cloudy["97.5%"],
color='red',alpha=0.3)
# 雨の日のビールの売り上げの平均値を描画
plt.plot(result_rainy['temperature'], result_rainy["50%"], color='green')
# 95%ベイズ信頼区間を描画
plt.fill_between(x=result_rainy['temperature'],
y1=result_rainy["2.5%"],
y2=result_rainy["97.5%"],
color='green',alpha=0.3)
# 晴れの日のビールの売り上げの平均値を描画
plt.plot(result_sunny['temperature'], result_sunny["50%"], color='blue')
# 95%ベイズ信頼区間を描画
plt.fill_between(x=result_sunny['temperature'],
y1=result_sunny["2.5%"],
y2=result_sunny["97.5%"],
color='blue',alpha=0.3)
# グラフタイトル追加
plt.title('ビールの売り上げと気温・天気の関係')
# グラフの描画
plt.show()
pystanの詳細については公式ページを参照してください。
PyStan — pystan 3.10.0 documentation
pystan.readthedocs.io
コメント