「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
コメント