pystanでトレンドの構造

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

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


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


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


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

DataAnalOji.hatena.sample/python_samples/stan/5-5-トレンドの構造.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 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'

データの読み込みと図示

データの読み込み

sales_df_3 = pd.read_csv('5-5-1-sales-ts-3.csv')
sales_df_3['date'] = pd.to_datetime(sales_df_3['date'])
sales_df_3.head(n=3)

図示

plt.figure(figsize=(15,5))
plt.plot(sales_df_3['sales'].values, color='black')
plt.show()

ローカルレベルモデルの推定

データの準備

data_list = dict(y=sales_df_3['sales'], T=len(sales_df_3))

ローカルレベルモデルの推定

# stanコードの記述
stan_code = '''
data {
  int T;         // データ取得期間の長さ
  vector[T] y;   // 観測値
}

parameters {
  vector[T] mu;       // 状態の推定値(水準成分)
  real<lower=0> s_w;  // 過程誤差の標準偏差
  real<lower=0> s_v;  // 観測誤差の標準偏差
}

model {
  // 状態方程式に従い、状態が遷移する
  for(i in 2:T) {
    mu[i] ~ normal(mu[i-1], s_w);
  }
  
  // 観測方程式に従い、観測値が得られる
  for(i in 1:T) {
    y[i] ~ normal(mu[i], s_v);
  }
}

'''

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

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

ローカルレベルモデルの推定結果

print(local_level.stansummary(pars=['s_w', 's_v', 'lp__'],
                              probs=[0.025, 0.5, 0.975]))
Inference for Stan model: anon_model_5912d60dc72bb8ffe174964dc2cf39dd.
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
s_w    4.27    0.04   0.83   2.94   4.18   6.12    477    1.0
s_v    7.58    0.01   0.72   6.21   7.56   9.09   2710    1.0
lp__ -438.7    0.72  14.66 -467.8 -438.6 -411.6    415    1.0

Samples were drawn using NUTS at Thu Sep  3 21:54:42 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).

平滑化トレンドモデルの推定

平滑化トレンドモデルの推定

# stanコードの記述
stan_code = '''
data {
  int T;        // データ取得期間の長さ
  vector[T] y;  // 観測値
}

parameters {
  vector[T] mu;       // 水準+ドリフト成分の推定値
  real<lower=0> s_z;  // ドリフト成分の変動の大きさを表す標準偏差
  real<lower=0> s_v;  // 観測誤差の標準偏差
}

model {
  // 状態方程式に従い、状態が遷移する
  for(i in 3:T) {
    mu[i] ~ normal(2 * mu[i-1] - mu[i-2], s_z);
  }
  
  // 観測方程式に従い、観測値が得られる
  for(i in 1:T) {
    y[i] ~ normal(mu[i], s_v);
  }

}

'''

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

# サンプリング
smooth_trend = stan_model_st.sampling(data=data_list,
                                      seed=1,
                                      n_jobs=1,
                                      iter=8000,
                                      warmup=2000,
                                      thin=6,
                                      control={
                                          'adapt_delta': 0.9,
                                          'max_treedepth': 15
                                      })

平滑化トレンドモデルの推定結果

print(smooth_trend.stansummary(pars=['s_z', 's_v', 'lp__'],
                               probs=[0.025, 0.5, 0.975]))
Inference for Stan model: anon_model_7270b55e94b176c72908fca84c5355ab.
4 chains, each with iter=8000; warmup=2000; thin=6; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

       mean se_mean     sd   2.5%    50%  97.5%  n_eff   Rhat
s_z    0.26  6.4e-3   0.12   0.13   0.23   0.58    349   1.01
s_v    8.43    0.01   0.64   7.27   8.39    9.8   3891    1.0
lp__ -171.1    1.95  36.92 -256.2 -166.8 -110.2    357   1.01

Samples were drawn using NUTS at Thu Sep  3 22:05:50 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).

ローカル線形トレンドモデルの推定

ローカル線形トレンドモデルの推定

# stanコードの記述
stan_code = '''
data {
  int T;        // データ取得期間の長さ
  vector[T] y;  // 観測値
}

parameters {
  vector[T] mu;       // 水準+ドリフト成分の推定値
  vector[T] delta;    // ドリフト成分の推定値
  real<lower=0> s_w;  // 水準成分の変動の大きさを表す標準偏差
  real<lower=0> s_z;  // ドリフト成分の変動の大きさを表す標準偏差
  real<lower=0> s_v;  // 観測誤差の標準偏差
}

model {
  // 弱情報事前分布
  s_w ~ normal(2, 2);
  s_z ~ normal(0.5, 0.5);
  s_v ~ normal(10, 5);
  
  // 状態方程式に従い、状態が遷移する
  for(i in 2:T) {
    mu[i] ~ normal(mu[i-1] + delta[i-1], s_w);
    delta[i] ~ normal(delta[i-1], s_z);
  }
  
  // 観測方程式に従い、観測値が得られる
  for(i in 1:T) {
    y[i] ~ normal(mu[i], s_v);
  }

}

'''

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

# サンプリング
local_linear_trend = stan_model_llt.sampling(data=data_list,
                                            seed=1,
                                            n_jobs=1,
                                            iter=8000,
                                            warmup=2000,
                                            thin=6)

ローカル線形トレンドモデルの推定結果

print(local_linear_trend.stansummary(pars=['s_w', 's_z', 's_v', 'lp__'],
                                     probs=[0.025, 0.5, 0.975]))
Inference for Stan model: anon_model_ddb0a2063558b071e59ccdc656993435.
4 chains, each with iter=8000; warmup=2000; thin=6; 
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    1.46    0.07   0.94   0.12    1.3   3.68    186   1.02
s_z    0.29  7.8e-3   0.13   0.13   0.26   0.63    293   1.01
s_v    8.29    0.02   0.66    7.1   8.26   9.63   1915    1.0
lp__ -235.7   11.45  89.63 -387.5 -246.1 -19.92     61   1.04

Samples were drawn using NUTS at Thu Sep  3 22:20:58 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_local = pd.DataFrame(local_level.summary()['summary'],
                             columns=local_level.summary()['summary_colnames'],
                             index=local_level.summary()['summary_rownames'])
summary_smooth = pd.DataFrame(smooth_trend.summary()['summary'],
                              columns=smooth_trend.summary()['summary_colnames'],
                              index=smooth_trend.summary()['summary_rownames'])
summary_linear = pd.DataFrame(local_linear_trend.summary()['summary'],
                              columns=local_linear_trend.summary()['summary_colnames'],
                              index=local_linear_trend.summary()['summary_rownames'])

# プロット
mcmc_rhat(summary_local)
mcmc_rhat(summary_smooth)
mcmc_rhat(summary_linear)
print('hmc_diagnostics of local_level:\n',pystan.diagnostics.check_hmc_diagnostics(local_level))
print('hmc_diagnostics of smooth_trend:\n',pystan.diagnostics.check_hmc_diagnostics(smooth_trend))
print('hmc_diagnostics of local_linear_trend:\n',pystan.diagnostics.check_hmc_diagnostics(local_linear_trend))
hmc_diagnostics of local_level:
 {'n_eff': True, 'Rhat': True, 'divergence': True, 'treedepth': True, 'energy': True}

hmc_diagnostics of smooth_trend:
 {'n_eff': True, 'Rhat': True, 'divergence': True, 'treedepth': True, 'energy': True}

WARNING:pystan:5 of 4000 iterations ended with a divergence (0.125 %).
WARNING:pystan:Try running with adapt_delta larger than 0.8 to remove the divergences.
WARNING:pystan:699 of 4000 iterations saturated the maximum tree depth of 10 (17.5 %)
WARNING:pystan:Run again with max_treedepth larger than 10 to avoid saturation
WARNING:pystan:Chain 1: E-BFMI = 0.0958
WARNING:pystan:Chain 2: E-BFMI = 0.191
WARNING:pystan:Chain 3: E-BFMI = 0.176
WARNING:pystan:E-BFMI below 0.2 indicates you may need to reparameterize your model

hmc_diagnostics of local_linear_trend:
 {'n_eff': True, 'Rhat': True, 'divergence': False, 'treedepth': False, 'energy': False}

参考:推定結果一覧

print(local_level.stansummary(probs=[0.025, 0.5, 0.975]))
print(smooth_trend.stansummary(probs=[0.025, 0.5, 0.975]))
print(local_linear_trend.stansummary(probs=[0.025, 0.5, 0.975]))
Inference for Stan model: anon_model_5912d60dc72bb8ffe174964dc2cf39dd.
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
mu[1]    91.29    0.07   4.94  81.33  91.48 100.68   4999    1.0
mu[2]    90.39    0.08   4.45  81.42  90.47  98.81   3494    1.0
…省略…

参考:トレースプロット

arviz.plot_trace(local_level,
                 var_names=['s_w', 's_v'],
                 legend=True);
arviz.plot_trace(smooth_trend,
                 var_names=['s_z', 's_v'],
                 legend=True);
arviz.plot_trace(local_linear_trend,
                 var_names=['s_w', 's_z', 's_v'],
                 legend=True);

推定された状態の図示

MCMCサンプルの取得

mcmc_sample_ll = local_level.extract()
mcmc_sample_st = smooth_trend.extract()
mcmc_sample_llt = local_linear_trend.extract()

補足:全時刻の95%ベイズ信用区間と中央値

# パラメータを指定
state_name = 'mu'

# 2.5・50・97.5パーセンタイルの値を
# 入れるためのデータフレームを作成
result_ll = pd.DataFrame(np.zeros([len(sales_df_3), 3]))
result_st = pd.DataFrame(np.zeros([len(sales_df_3), 3]))
result_llt = pd.DataFrame(np.zeros([len(sales_df_3), 3]))

# 2.5・50・97.5パーセンタイルの値をデータフレームに格納
for i in range(len(sales_df_3)):
    result_ll.iloc[i, :] = np.percentile(mcmc_sample_ll[state_name][:, i],
                                         q=[2.5, 50, 97.5])
    result_st.iloc[i, :] = np.percentile(mcmc_sample_st[state_name][:, i],
                                         q=[2.5, 50, 97.5])
    result_llt.iloc[i, :] = np.percentile(mcmc_sample_llt[state_name][:, i],
                                          q=[2.5, 50, 97.5])

# 列名の変更
result_ll.columns = ["lwr", "fit", "upr"]
result_st.columns = ["lwr", "fit", "upr"]
result_llt.columns = ["lwr", "fit", "upr"]

# 時間情報の追加
result_ll['time'] = sales_df_3['date']
result_st['time'] = sales_df_3['date']
result_llt['time'] = sales_df_3['date']

# 実測値を追加
result_ll['obs'] = sales_df_3['sales']
result_st['obs'] = sales_df_3['sales']
result_llt['obs'] = sales_df_3['sales']

図示

# 描画領域の作成
fig = plt.figure(figsize=(15, 15))

# ローカルレベルモデル
ax_ll = fig.add_subplot(3, 1, 1)
# 実測値をプロット
ax_ll.plot(result_ll['time'],
           result_ll['obs'],
           marker='.',
           linewidth=0,
           color='black')
# 50パーセンタイルの値をプロット
ax_ll.plot(result_ll['time'], result_ll['fit'], color='black')
# 2.5-97.5%区間の間を埋める
ax_ll.fill_between(x=result_ll['time'],
                   y1=result_ll['upr'],
                   y2=result_ll['lwr'],
                   color='gray',
                   alpha=0.5)
# x軸ラベルの設定
ax_ll.set_xlabel('time')
# y軸ラベルの設定
ax_ll.set_ylabel('sales')
# タイトルの設定
ax_ll.set_title('ローカルレベルモデル')

# 平滑化トレンドモデル
ax_st = fig.add_subplot(3, 1, 2)
# 実測値をプロット
ax_st.plot(result_st['time'],
           result_st['obs'],
           marker='.',
           linewidth=0,
           color='black')
# 50パーセンタイルの値をプロット
ax_st.plot(result_st['time'], result_st['fit'], color='black')
# 2.5-97.5%区間の間を埋める
ax_st.fill_between(x=result_st['time'],
                   y1=result_st['upr'],
                   y2=result_st['lwr'],
                   color='gray',
                   alpha=0.5)
# x軸ラベルの設定
ax_st.set_xlabel('time')
# y軸ラベルの設定
ax_st.set_ylabel('sales')
# タイトルの設定
ax_st.set_title('平滑化トレンドモデル')

# ローカル線系トレンドモデル
ax_llt = fig.add_subplot(3, 1, 3)
# 実測値をプロット
ax_llt.plot(result_llt['time'],
            result_llt['obs'],
            marker='.',
            linewidth=0,
            color='black')
# 50パーセンタイルの値をプロット
ax_llt.plot(result_llt['time'], result_llt['fit'], color='black')
# 2.5-97.5%区間の間を埋める
ax_llt.fill_between(x=result_llt['time'],
                    y1=result_llt['upr'],
                    y2=result_llt['lwr'],
                    color='gray',
                    alpha=0.5)
# x軸ラベルの設定
ax_llt.set_xlabel('time')
# y軸ラベルの設定
ax_llt.set_ylabel('sales')
# タイトルの設定
ax_llt.set_title('ローカル線形トレンドモデル')

# グラフの描画
fig.tight_layout()
fig.show()

ドリフト成分の図示

補足:全時刻の95%ベイズ信用区間と中央値

# パラメータを指定
state_name = 'delta'

# 2.5・50・97.5パーセンタイルの値を
# 入れるためのデータフレームを作成
result_llt = pd.DataFrame(np.zeros([len(sales_df_3), 3]))

# 2.5・50・97.5パーセンタイルの値をデータフレームに格納
for i in range(len(sales_df_3)):
    result_llt.iloc[i, :] = np.percentile(mcmc_sample_llt[state_name][:, i],
                                          q=[2.5, 50, 97.5])
# 列名の変更
result_llt.columns = ["lwr", "fit", "upr"]

# 時間情報の追加
result_llt['time'] = sales_df_3['date']

# 実測値を追加
result_llt['obs'] = sales_df_3['sales']

図示

# 描画領域の作成
plt.figure(figsize=(15, 5))
# 50パーセンタイルの値をプロット
plt.plot(result_llt['time'], 
         result_llt['fit'], 
         color='black')
# 2.5-97.5%区間の間を埋める
plt.fill_between(x=result_llt['time'],
                 y1=result_llt['upr'],
                 y2=result_llt['lwr'],
                 color='gray',
                 alpha=0.5)
# x軸ラベルの設定
plt.xlabel('time')
# y軸ラベルの設定
plt.ylabel('sales')
# タイトルの設定
plt.title('ドリフト成分')
# グラフの描画
plt.show()

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

PyStan — pystan 3.9.1 documentation

コメント