pystanでデザイン行列を用いた一般化線形モデルの推定

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

「RとStanで始めるベイズ統計モデリングによるデータ分析入門」の「第3部第4章 デザイン行列を用いた一般化線形モデルの推定」を対象に,公開されているR,Stanのコードをpython,pystanのコードへと書き直した一例です。Stanの代わりにpystanを利用しています。

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

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

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

DataAnalOji.hatena.sample/python_samples/stan/3-4-デザイン行列を用いた一般化線形モデルの推定.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 pandas as pd

%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('ggplot')
plt.rcParams['font.family'] = 'Meiryo'

分析対象のデータ

# csvの読み込み(dataframe型)
file_beer_sales_2 = pd.read_csv('3-2-1-beer-sales-2.csv')

サンプルサイズ

# len()でデータフレームの行数を取得
sample_size = len(file_beer_sales_2)

pandas.dataframe.insert()用いたデザイン行列の作成

デザイン行列の作成

# 応答変数をデータフレームから削除
X = file_beer_sales_2.drop(['sales'], axis=1)
# 先頭列に値が1の(Intercept)を追加
X.insert(0, '(Intercept)', 1)

pandas.dataframe.insert()つかったデザイン行列

# 最初の5行を表示
X.head(n=5)

MCMCの実行

サンプルサイズ

N = len(file_beer_sales_2)

デザイン行列の列数(説明変数の数+1)

K = X.shape[1]

応答変数

Y = file_beer_sales_2['sales']

dictにまとめる

rstanでは入力データをlist形式にまとめたのに対し,pystanでは入力データをdictionary形式にまとめる。

data_list_design = dict(N=N, K=K, Y=Y, X=X)

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)

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

結果の表示

print(mcmc_result_design.stansummary(probs=[0.025, 0.5, 0.975]))
Inference for Stan model: anon_model_40bea9213a37e125c59d5a21d95a3a94.
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]    21.1    0.16   5.98   9.22  21.03  32.89   1426    1.0
b[2]    2.46  7.6e-3   0.29   1.89   2.47   3.04   1438    1.0
sigma  17.07    0.03   1.23  14.85   17.0   19.6   2089    1.0
lp__  -330.1    0.03   1.24 -333.3 -329.8 -328.7   1519    1.0

Samples were drawn using NUTS at Tue Jul 28 23:20:23 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

コメント