pystanでモデルを用いた予測

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

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

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

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

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

DataAnalOji.hatena.sample/python_samples/stan/3-3-モデルを用いた予測.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.

分析の準備
予測のためのデータの整理
1. 気温を11度から30度まで変化させて、その時の売り上げを予測する
2. dictにまとめる
MCMCの実行
1. MCMCの実行
2. 結果の表示
予測分布の図示

分析の準備

パッケージの読み込み

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'

分析対象のデータ

file_beer_sales_2 = pd.read_csv('3-2-1-beer-sales-2.csv')

サンプルサイズ

sample_size = len(file_beer_sales_2)

予測のためのデータの整理

気温を11度から30度まで変化させて、その時の売り上げを予測する

temperature_pred = np.linspace(start=11, stop=30, num=20)
print(temperature_pred)


[11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28.
 29. 30.]

dictにまとめる

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

data_list_pred = dict(N=sample_size,
                      sales=file_beer_sales_2['sales'],
                      temperature=file_beer_sales_2['temperature'],
                      N_pred=len(temperature_pred),
                      temperature_pred=temperature_pred)

MCMCの実行

mcmc_result_pred = pystan.stan(file='3-3-1-simple-lm-pred.stan',
                               data=data_list_pred,
                               seed=1,
                               n_jobs=1)

結果の表示

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

Inference for Stan model: anon_model_b9ba45baf41873b2dc42b9a64b7e03ff.
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
Intercept       20.85    0.16   5.94    8.8   20.9  32.12   1431    1.0
beta             2.48  7.8e-3   0.29   1.92   2.47   3.07   1395    1.0
sigma           17.08    0.03   1.22  14.94  17.01  19.66   1995    1.0
mu_pred[1]      48.09    0.07   3.03  42.09  48.15  53.77   1668    1.0
省略・・・
mu_pred[20]     95.16    0.08   3.46  88.39  95.16 101.94   1756    1.0
sales_pred[1]   48.69    0.27   17.2  14.83  48.56  83.47   3979    1.0
省略・・・
sales_pred[20]  95.19    0.29  17.46  60.65  95.29  129.4   3602    1.0
lp__           -330.1    0.03   1.22 -333.4 -329.8 -328.7   1356    1.0

Samples were drawn using NUTS at Tue Jul 14 12:22:40 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).

予測分布の図示

MCMCサンプルの抽出

ただし，bayseplotの代わりにarvizを利用するので抽出したMCMCサンプルは利用しない。

mcmc_sample_pred = mcmc_result_pred.extract()

気温が11度～30度まで1度ずつ変えたの時の予測売り上げの95%予測区間の図示

arviz.plot_forest(mcmc_result_pred,
                  var_names=["sales_pred"],
                  hdi_prob=0.95,
                  combined=True,
                  colors='b');

95%区間の比較

arviz.plot_forest(mcmc_result_pred,
                  var_names=['mu_pred', 'sales_pred'],
                  coords={'mu_pred_dim_0': 0,
                          'sales_pred_dim_0': 0},
                  hdi_prob=0.95,
                  combined=True,
                  colors='b');

気温が11度と30度の時の、売り上げの予測分布

arviz.plot_forest(mcmc_result_pred,
                  var_names=['sales_pred'],
                  coords={'sales_pred_dim_0': [0, 19]},
                  kind='ridgeplot',
                  hdi_prob=0.95,
                  combined=True,
                  ridgeplot_alpha=0.8,
                  ridgeplot_quantiles=[.25, .5, .75],
                  colors='blue');