pystanでダミー変数と分散分析モデル

実践Ｄａｔａ　Ｓｃｉｅｎｃｅシリーズ　ＲとＳｔａｎではじめる　ベイズ統計モデリングによるデータ分析入門 (ＫＳ情報科学専門書)

「RとStanで始めるベイズ統計モデリングによるデータ分析入門実践編第6章ダミー変数と分散分析モデル」を対象に，公開されているR，Stanのコードをpython，pystanのコードへと書き直した一例です。Stanの代わりにpystanを利用しています。

本ページでは公開されていない書籍の内容については一切触れません。理論や詳しい説明は書籍を参照してください。

なお，こちらで紹介しているコードには誤りが含まれる可能性があります。内容やコードについてお気づきの点等ございましたら，ご指摘いただけると幸いです。

（本ページにて紹介しているコードはgithubにて公開しています。）

DataAnalOji.hatena.sample/python_samples/stan/3-6-ダミー変数と分散分析モデル.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_weather = pd.read_csv('3-6-1-beer-sales-3.csv')
sales_weather.head(n=3)

データの要約

display(sales_weather.describe(include='all'))

図示

plt.figure()
sns.violinplot(x='weather',
               y='sales',
               data=sales_weather,
               inner='point',
               order=['cloudy', 'rainy', 'sunny'],
               color='white')
sns.stripplot(x='weather',
              y='sales',
              data=sales_weather,
              order=['cloudy', 'rainy', 'sunny'],
              palette=['r', 'g', 'b'],
              jitter=0)
plt.title('ビールの売り上げと天気の関係')
plt.show()

brmsによる分散分析モデルの推定

pythonにbrmsの代用可能なパッケージが知っている範囲内で存在しないため，brmsを利用したコードは省略しています。

分散分析モデルのデザイン行列

デザイン行列の作成

# 天気をダミー変数化
design_mat = pd.get_dummies(sales_weather['weather'], drop_first=True)
# 先頭列に(Intercept)を追加
design_mat.insert(0, '(Intercept)', 1)
design_mat.head(n=3)

stanに渡すdictの作成

data_list = dict(N=len(sales_weather),
                 K=design_mat.shape[1],
                 Y=sales_weather['sales'],
                 X=design_mat)

stanに渡すデータの表示

data_list

{'N': 150,
 'K': 3,
 'Y': 0      48.5
 1      64.8
 2      85.8
 3      45.0
 4      60.8
        ... 
 145    50.2
 146    69.9
 147    68.2
 148    47.9
 149    45.0
 Name: sales, Length: 150, dtype: float64,
 'X':      (Intercept)  rainy  sunny
 0              1      0      0
 1              1      0      0
 2              1      0      0
 3              1      0      0
 4              1      0      0
 ..           ...    ...    ...
 145            1      1      0
 146            1      1      0
 147            1      1      0
 148            1      1      0
 149            1      1      0
 
 [150 rows x 3 columns]}

補足：brmsを使わない分散分析モデルの推定

pystanで分散分析モデルを実行

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

# サンプリング
anova_stan = stan_model.sampling(data=data_list, seed=1, n_jobs=1)

結果の確認

print(anova_stan.stansummary(probs=[0.025, 0.5, 0.975]))

Inference for Stan model: anon_model_ff3814f23c22278ffe4dee1df56226cc.
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]   63.11    0.06   2.45  58.36  63.06  68.09   1880    1.0
b[2]   -0.45    0.08   3.44  -7.39  -0.39   6.11   2076    1.0
b[3]   19.91    0.07   3.38  13.26  19.85  26.57   2066    1.0
sigma  16.93    0.02    1.0  15.13  16.87  19.13   3313    1.0
lp__  -495.7    0.03   1.44 -499.3 -495.3 -493.9   1892    1.0

Samples were drawn using NUTS at Mon Jul 20 21:46:33 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).

補足：推定された天気別の平均売り上げのグラフ

推定用のデザイン行列の作成

# 1列目は(Intercept)に対応
design_cloudy = [[1, 0, 0]]
design_rainy = [[1, 1, 0]]
design_sunny = [[1, 0, 1]]
pred_mat = np.concatenate((design_cloudy, 
                           design_rainy, 
                           design_sunny))
print(pred_mat)

[[1 0 0]
 [1 1 0]
 [1 0 1]]

stanに渡すdictの作成

# stanに渡すdictの作成
data_list = dict(N=len(sales_weather),
                 K=design_mat.shape[1],
                 Y=sales_weather['sales'],
                 X=design_mat,
                 X_pred=pred_mat)

pystanで分散分析モデルを実行し予測事後分布を推定

# stanコードの記述
stan_code = '''
data { 
    int N;                // サンプルサイズ
    int K;                // デザイン行列の列数(説明変数の数＋１)
    vector[N] Y;          // 応答変数 
    matrix[N, K] X;       // デザイン行列 
    matrix[K, K] X_pred;    // 予測デザイン行列
} 

parameters { 
    vector[K] b;          // 切片を含む係数ベクトル
    real<lower=0> sigma;  // データのばらつきを表す標準偏差
} 

model { 
    vector[N] mu = X * b;
    Y ~ normal(mu, sigma);
} 

generated quantities {
    vector[K] mu_pred = X_pred * b;
}

'''

# モデルのコンパイル
stan_model = pystan.StanModel(model_code=stan_code)

# サンプリング
anova_stan = stan_model.sampling(data=data_list, seed=1, n_jobs=1)

生成された乱数を格納

mcmc_sample= anova_stan.extract()

推定された天気別の平均売り上げのグラフ

plt.figure()
plt.boxplot(mcmc_sample['mu_pred'],
            labels=['cloudy', 'rainy', 'sunny'],
            showfliers=False)
plt.show()

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

PyStan — pystan 3.10.0 documentation