実践Data Scienceシリーズ RとStanではじめる ベイズ統計モデリングによるデータ分析入門 (KS情報科学専門書)
「RとStanで始めるベイズ統計モデリングによるデータ分析入門」「実践編第5部第8章 動的一般化線形モデル：二項分布を仮定した例」を対象に，公開されているR，Stanのコードをpython，pystanのコードへと書き直した一例です。Stanの代わりにpystanを利用しています。
この章では，観測方程式に二項分布を用いた動的一般化線形モデル（DGLM)が解説されています。
本ページでは公開されていない書籍の内容については一切触れません。理論や詳しい説明は書籍を参照してください。
なお，こちらで紹介しているコードには誤りが含まれる可能性があります。内容やコードについてお気づきの点等ございましたら，ご指摘いただけると幸いです
（本ページにて紹介しているコードはgithubにて公開しています。）
DataAnalOji.hatena.sample/python_samples/stan/5-8-動的一般化線形モデル：二項分布を仮定した例.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
分析の準備
パッケージの読み込み
plotSSM関数についてはこちらをご参照ください。
import arviz
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
# 自作のplotSSM関数を読み込み
from plotSSM import plotSSMデータの読み込み
こちらで読み込んでいる「boat.csv」はR上で読み込んだ「boat」データセットを，「write.csv(boat, ‘boat.csv’)」によってcsvファイルに出力したものを利用しています。
boat = pd.read_csv('boat.csv',
header=0, # 最初の行をheader指定
index_col=0, # 最小の列をindex指定
names=['victory or defeat'])# 列名変更
display(boat)
二項分布を仮定したDGLMの推定
参考
print(np.logical_not(boat['victory or defeat'].isnull())) # データがあればTrue1 True
2 False
3 False
4 False
5 False
...
179 True
180 True
181 True
182 True
183 True
Name: victory or defeat, Length: 183, dtype: boolind = boat[np.logical_not(boat['victory or defeat'].isnull())].index
print(ind) # データがある時点一覧Int64Index([ 1, 8, 11, 12, 13, 14, 17, 18, 21, 24,
...
174, 175, 176, 177, 178, 179, 180, 181, 182, 183],
dtype='int64', length=155)NAを除く
dropna(axis=0)でNaN値の含まれる行を削除します。
boat_omit_NA = boat.dropna(axis=0)データの準備
data_list = dict(T=len(boat),
len_obs=len(boat_omit_NA),
y=boat_omit_NA['victory or defeat'].astype(int), # stanにデータを渡すためint型へ変換
obs_no=ind
)モデルの推定
# stanコードの記述（5-8-1-dglm-binom.stan）
stan_code = '''
data {
int T; // データ取得期間の長さ
int len_obs; // 観測値が得られた個数
int y[len_obs]; // 観測値
int obs_no[len_obs]; // 観測値が得られた時点
}
parameters {
vector[T] mu; // 状態の推定値
real<lower=0> s_w; // 過程誤差の標準偏差
}
model {
// 弱情報事前分布
s_w ~ student_t(3, 0, 10);
// 状態方程式に従い、状態が遷移する
for(i in 2:T) {
mu[i] ~ normal(mu[i-1], s_w);
}
// 観測方程式に従い、観測値が得られる
// ただし、「観測値が得られた時点」でのみ実行する
for(i in 1:len_obs) {
y[i] ~ bernoulli_logit(mu[obs_no[i]]);
}
}
generated quantities{
vector[T] probs; // 推定された勝率
probs = inv_logit(mu);
}
'''
# モデルのコンパイル
stan_model = pystan.StanModel(model_code=stan_code)
# サンプリング
dglm_binom = stan_model.sampling(data=data_list,
seed=1,
iter=30000,
warmup=10000,
thin=20,
n_jobs=1)推定されたパラメタ
print(dglm_binom.stansummary(pars=["s_w", "lp__"], probs=[0.025, 0.5, 0.975]))Inference for Stan model: anon_model_7ddd0d92330219129ee77bde56cd3307.
4 chains, each with iter=30000; warmup=10000; thin=20;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
s_w 0.74 0.02 0.39 0.21 0.68 1.71 582 1.01
lp__ -96.76 3.76 86.34 -252.7 -104.4 88.46 527 1.01
Samples were drawn using NUTS at Sun Sep 13 15:30:56 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).参考：収束の確認
# 収束確認用のRhatのプロット関数
def mcmc_rhat(dataframe, column='Rhat', figsize=(5, 10)):
plt.figure(figsize=figsize)
plt.hlines(y=dataframe[column].sort_values().index,
xmin=1,
xmax=dataframe[column].sort_values(),
color='b',
alpha=0.5)
plt.vlines(x=1.05, ymin=0, ymax=len(dataframe[column]), linestyles='--')
plt.plot(dataframe[column].sort_values().values,
dataframe[column].sort_values().index,
marker='.',
linestyle='None',
color='b',
alpha=0.5)
plt.yticks(color='None')
plt.tick_params(length=0)
plt.xlabel(column)
plt.show()
# 各推定結果のデータフレームを作成
summary = pd.DataFrame(dglm_binom.summary()['summary'],
columns=dglm_binom.summary()['summary_colnames'],
index=dglm_binom.summary()['summary_rownames'])
# プロット
mcmc_rhat(summary)
print('hmc_diagnostics:\n',
pystan.diagnostics.check_hmc_diagnostics(dglm_binom))hmc_diagnostics of local_level:
{'n_eff': True, 'Rhat': True, 'divergence': True, 'treedepth': True, 'energy': True}参考:トレースプロット
’lp__’（log posterior）のトレースプロットは図示できないため除いています。
arviz.plot_trace(dglm_binom, var_names=["s_w"], legend=True);
参考：推定結果一覧
print(dglm_binom.stansummary(probs=[0.025, 0.5, 0.975]))Inference for Stan model: anon_model_7ddd0d92330219129ee77bde56cd3307.
4 chains, each with iter=30000; warmup=10000; thin=20;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
mu[1] -0.4 0.04 1.95 -4.85 -0.2 2.91 2139 1.0
mu[2] -0.17 0.03 1.91 -4.45 -0.03 3.17 3020 1.0
…省略…
mu[182] -0.59 0.02 1.09 -2.86 -0.55 1.52 3809 1.0
mu[183] -0.81 0.02 1.24 -3.62 -0.74 1.5 3460 1.0
s_w 0.74 0.02 0.39 0.21 0.68 1.71 582 1.01
probs[1] 0.45 6.3e-3 0.29 7.8e-3 0.45 0.95 2150 1.0
probs[2] 0.49 5.2e-3 0.29 0.01 0.49 0.96 3046 1.0
…省略…
probs[182] 0.38 3.3e-3 0.21 0.05 0.37 0.82 3808 1.0
probs[183] 0.35 3.4e-3 0.21 0.03 0.32 0.82 3910 1.0
lp__ -96.76 3.76 86.34 -252.7 -104.4 88.46 527 1.01
Samples were drawn using NUTS at Sun Sep 13 15:30:56 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).推定された状態の図示
時間ラベルの作成
years = pd.date_range('1829-01-01', periods=len(boat), freq='12MS')
print(years[:3])DatetimeIndex(['1829-01-01', '1830-01-01', '1831-01-01'], dtype='datetime64[ns]', freq='12MS')MCMCサンプルの取得
mcmc_sample = dglm_binom.extract()ケンブリッジ大学の勝率の推移のグラフ
fig, ax = plt.subplots(1, 1, figsize=(10, 5))
p = plotSSM(mcmc_sample=mcmc_sample,
time_vec=years,
obs_vec=boat_omit_NA['victory or defeat'].astype(int),
state_name='probs',
graph_title='ケンブリッジ大学の勝率の推移',
y_label='勝率',
axes=ax)
plt.show()
ケンブリッジ大学の平均勝率
boat_omit_NA['victory or defeat'].mean()0.5161290322580645pystanの詳細については公式ページを参照してください。
PyStan — pystan 3.9.1 documentation
pystan.readthedocs.io
コメント