pystanで正規線形モデル

実践Data Scienceシリーズ RとStanではじめる ベイズ統計モデリングによるデータ分析入門 (KS情報科学専門書)

「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.

分析の準備

パッケージの読み込み

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)

f:id:RYHR:20200728232735p:plain

データの要約

sales_climate.describe(include='all')
f:id:RYHR:20200728232811p:plain

図示

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()

f:id:RYHR:20200728232839p:plain

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)
f:id:RYHR:20200728232917p:plain

補足: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;                // デザイン行列の列数(説明変数の数+1)
  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;                // デザイン行列の列数(説明変数の数+1)
  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()
f:id:RYHR:20200728233328p:plain

pystanの詳細については公式ページを参照してください。

PyStan — pystan 3.9.1 documentation

コメント