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.

分析の準備

パッケージの読み込み

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の実行

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');

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

PyStan — pystan 3.9.1 documentation
Redirecting to new ArviZ documentation host: ReadTheDocs

コメント