「ベイズ統計モデリングによるデータ分析入門」をPythonとStanで写経 ～ Vol.9 単回帰モデルを用いた予測

書籍の著者　馬場真哉 先生

この記事は、書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」第３部第３章「モデルを用いた予測」の Python 写経活動記録です。

書籍の第３部は「一般化線形モデル」のベイズ統計モデリングです。
「単回帰モデル」から始めます。
前回記事の単回帰モデルの構築から続いて、今回記事は単回帰モデルによる予測を行います。

では書籍を開いてベイズ統計モデリングの旅に出かけましょう🚀

はじめに

このブログシリーズは書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」（講談社、「テキスト」と呼びます）の Python 写経です。

テキストの紹介と引用表記はリンク先の記事に掲載しています。

準備

■ 記事の範囲
この記事はテキスト第３部第３章の以下の節を取り扱います。

ベイズ統計モデリングの理論面が気になった場合には、ぜひテキストの第１部、第３部第１章をご覧いただき、本記事との繋がりをご確認下さいませ。

■ コード記述法
Jupyter Notebook 形式でコードを記述します。

■ 利用データ
テキスト・サポートサイトのデータファイルを引用しています。
Jupyter Notebook ファイルと同一フォルダ内の「data」フォルダにデータファイルを格納しています。

■ ライブラリのインポート
この記事で用いるライブラリをインポートします。

# インポート

# 数値計算
import numpy as np
import pandas as pd

# ベイズ統計モデリング
from cmdstanpy import CmdStanModel     # stan
import arviz as az                     # 分析・可視化

# 統計モデリング
import statsmodels.formula.api as smf

# ユーティリティ
import os

# 可視化
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_theme()                        # ggplot風のスタイル
plt.rcParams['font.family'] = 'Meiryo'

第３章 単回帰モデルを用いた予測

データの読み込みと外観の確認

テキスト 3.2 節に相当します。
テキストの仮想のビールの売り上げデータを引用いたします。
前回記事と同じデータです。
csv ファイルを pandas のデータフレーム形式で変数 file_beer_sales_2 に読み込みます。

# p.173 分析対象データの読み込み

# ファイルの読み込み
file_beer_sales_2 = pd.read_csv('./data/3-2-1-beer-sales-2.csv')

# 結果の表示
print('file_beer_sales_2.shape: ', file_beer_sales_2.shape)
file_beer_sales_2.head(3)

【実行結果】
標本サイズ 100、変数 sales は 売り上げ（単位：万円）、temperature は気温（おそらく摂氏℃）です。
気温とビール売り上げの関係を単回帰モデルで分析します。

データの要約統計量を確認します。

# データの要約統計量
file_beer_sales_2.describe().T.round(2)

【実行結果】

売り上げは平均 70、最小値 28、最大値 125 です。範囲が広い感じ。
気温は平均 20、最小値 10、最大値 30 です。
最近の真夏の気温と比べると最大値は小さい感じ。

データを可視化しましょう。
テキスト 図 3.2.1 に相当する散布図です（前回記事の散布図と同じです）。
seaborn ライブラリを利用します。

# p.168 図3.2.1 ビールの売上と気温の散布図

# 散布図の描画
plt.figure(figsize=(8, 4))
sns.scatterplot(data=file_beer_sales_2, x='temperature', y='sales')
plt.title('ビールの売上と気温の関係', loc='left')
plt.xticks(range(10, 31, 2));

【実行結果】
気温が高くなるにつれて売り上げが大きくなる傾向が見られます。

🔵🔵🔵

単回帰分析

まず単回帰分析による予測を確認しておき、あとでベイズ統計モデルの結果と比べてみましょう。
Python の統計ライブラリ statsmodels を利用して、次の単回帰モデルを実装します！

$${\varepsilon}$$ は誤差です。

この単回帰モデルを statsmodels の最小二乗法 olsにあたえる formula 構文に変換します。
$${\sim}$$ を挟んで、左辺は目的変数、右辺は説明変数です。

では単回帰分析を実行します！

# 単回帰分析 by statsmodels
formula = 'sales ~ temperature'
res_sm = smf.ols(formula=formula, data=file_beer_sales_2).fit()
res_sm.summary()

【実行結果】
こちらは（お馴染みの？）単回帰分析のサマリーです。

予測をしましょう。
単回帰分析の結果 res_sm に対して get_prediction メソッドを適用して、いろんな予測値を算出します！

# 回帰直線と95%予測区間の描画

## 予測
# 予測に用いるx:気温データの作成
x_min, x_max = np.sort(file_beer_sales_2.temperature)[[0, -1]]
x_vals = np.linspace(x_min, x_max, 100)

# 気温データを辞書化
x_dict = dict(temperature=x_vals)

# 回帰分析の結果を用いて予測を実行
preds = res_sm.get_prediction(exog=x_dict).summary_frame()
preds

【実行結果】
x_vals の 100 点に対する予測値たちです。

予測値の内容を表にまとめました。

（※）データ点のばらつきを含む予測値に基づく

データの散布図と予測値を重ね描きしましょう。
予測値平均 mean と 95% 予測区間を描画します。

# 描画

# 描画領域の設定
plt.figure(figsize=(8, 4))
# 回帰直線の描画
plt.plot(x_vals, preds['mean'], color='tab:red')
# 95% 予測区間の塗りつぶし描画
plt.fill_between(x_vals, preds.obs_ci_lower, preds.obs_ci_upper,
                 color='lightpink', alpha=0.4)
# 観測値の散布図の描画
sns.scatterplot(data=file_beer_sales_2, x='temperature', y='sales')
# 修飾
plt.title('ビールの売上と気温の関係：回帰直線・95%予測区間', loc='left')
plt.xticks(range(10, 31, 2));

【実行結果】
赤実線が予測値平均を結んだ回帰直線です。
薄赤色の塗りつぶし区間が 95% 予測区間です。
各観測値のデータ点は 95% 予測区間に含まれていますね！
95% 予測区間の幅は 60 くらいあり、かなり広い感じがいたします。

ではベイズ流の単回帰モデルへ進みます。

🔵🔵🔵

ベイズモデリング by Stan

テキスト 3.4 節、3.5 節に相当します。

① モデルの概要
次のモデルを model ブロックに実装します。

目的変数 $${sales}$$ は正規分布に従うと仮定しています。
パラメータ $${Intercept,\ \beta,\ \sigma}$$ には事前分布を明示的に設定しません。

続いて、予測分布からの乱数生成に必要な変数を定義して、generated quantities ブロックに実装します。

$${mu\_pred_i}$$ は予測平均値、$${sales\_pred_i}$$ は予測分布からのサンプルです。
$${temperature\_pred_i }$$ は予測に用いる気温データです。

② Stan のモデル設定
Stan ファイル（Stan コード）を作成します。
テキストの Stan ファイル名 および Stan コードを引用いたします。
model ブロックには「ベクトル化されたコード」を用いています。
📑ファイル名：3-3-1-simple-lm-pred.stan

data {
    int N;                  // 標本サイズ
    vector[N] sales;        // 売上データ
    vector[N] temperature;  // 気温データ

    int N_pred;                       // 予測対象データの大きさ
    vector[N_pred] temperature_pred;  // 予測対象となる気温
}

parameters {
    real Intercept;         // 切片
    real beta;              // 係数
    real<lower=0> sigma;    // 標準偏差
}

model {
    // 平均 Intercept + beta*temperature
    // 標準偏差 sigma の正規分布に従ってデータが得られたと仮定
    sales ~ normal(Intercept + beta * temperature, sigma);
}

generated quantities {
    vector[N_pred] mu_pred;     // ビールの売上の期待値
    vector[N_pred] sales_pred;  // ビールの売上の予測値

    for (i in 1:N_pred) {
        mu_pred[i] = Intercept + beta * temperature_pred[i];
        sales_pred[i] = normal_rng(mu_pred[i], sigma);
    }
}

③ データの整理
予測に用いる説明変数 temperature の値を設定します。

# p.174 気温を11度から30度まで変化させて、その時の売上を予測する
temperature_pred = np.arange(10, 31)  # 10度からスタート
temperature_pred

【実行結果】
テキスト 3.4 節とちょっぴり変えて、10度から30度まで１度刻みで用意します。
予測分布の可視化のときに左端のもやもやを回避したいからです。

標本サイズを算出して、Stan に渡すデータセットを辞書にまとめます。

# p.173 標本サイズ
sample_size = len(file_beer_sales_2)

# p.175 辞書にまとめる
data_dict = 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,
)

【実行結果】なし

④ モデルのコンパイル
モデルのコンパイルを実行します。

%%time
# モデルのコンパイル

# stanプログラムファイルのパス指定
current_dir = os.path.abspath(os.getcwd())
stan_file = os.path.join(current_dir, 'stan', '3-3-1-simple-lm-pred.stan')

# モデルオブジェクトの作成（exeの作成）
model = CmdStanModel(stan_file=stan_file)  # stanファイルを指定

【実行結果】（右側のファイルパスは記載省略）

MCMC の準備が整いました！

🔵🔵🔵

MCMC の実行

テキスト 3.6 節に相当します。
MCMCを実行しましょう。

%%time
# p.176 乱数の生成
fit = model.sample(
    data=data_dict,    # 対象データ
    seed=1,            # 乱数の種
    sig_figs=18,       # 出力CSV等に適用する数値精度
)

【実行結果】

収束の確認をします。
診断メソッド diagnose を利用します。

# 事後分布の診断
print(fit.diagnose())

【実行結果】（１行目のファイルパスは記載省略）
問題は検出されませんでした（no problems detected.）。

Stan の MCMC サンプルを arviz の idata 形式に変換して、具体的な値を見ていきましょう。

# arvizのidataに変換
idata = az.from_cmdstanpy(posterior=fit, log_likelihood='lp__')

【実行結果】なし

MCMC サンプルの要約表を表示します。

# p.176 要約統計量の表示
var_names = ['Intercept', 'beta', 'sigma']
az.summary(idata, var_names=var_names, hdi_prob=0.95)

【実行結果】
$${\widehat{R}}$$（r_hat）、有効サンプル数（ess_bulk、ess_tail）に問題はなさそうです。

トレースプロットを描画します。

# トレースプロットの図示
az.plot_trace(idata, compact=False, var_names=var_names,
              backend_kwargs={'tight_layout': True});

【実行結果】
右側のチャートがゲジゲジしています。

以上のチェックに基づいて、収束していると考えましょう。

🔵🔵🔵

予測分布の可視化

テキスト 3.7 節に相当します。

① 事後予測の 95% HDI の可視化
テキスト 図 3.3.1 に類似するチャートです。
arviz の plot_forest を利用して、気温 10 度から 30 度までの事後予測の 95% HDI の区間を描画します。
（分位点 quantile を用いるテキストの 95% 予測区間とは異なります）

# p.178 図3.3.1 95%予測区間 太い線は50%区間
az.plot_forest(idata.posterior.sales_pred, combined=True, hdi_prob=0.95);

【実行結果】
[0]：10度 ～ [20]：30 度までの 95% HDI が並んでいます。

② mu の 95% HDI と obs の 95% HDI の対比
テキスト 図 3.3.2 に類似するチャートです。
予測平均値 mu の 95% HDI よりも、観測値の誤差を含む予測値 obs の 95% HDI の方が幅が広くなる（誤差を含むので）ことを可視化で確かめます。
テキストに合わせて 11 度のケースを比べます。
引き続き plot_forest を利用します。

# p.179 図3.3.2 mu_predとsales_predの比較
ax = az.plot_forest(
    [idata.posterior.sales_pred.sel(sales_pred_dim_0=1),
     idata.posterior.mu_pred.sel(mu_pred_dim_0=1)],
    model_names=['sales', 'mu'], combined=True, hdi_prob=0.95)
ax[0].set_title('95%HDI：11度のケース');

【実行結果】
誤差を含む obs が、なんて幅広なのでしょう。

③ 異なる気温の予測分布を比較
テキスト 図 3.3.3 に相当します。
11 度と 30 度の予測分布を比べます。
引き続き plot_forest を利用します。

# p.179 図3.3.3 予測分布の図示
az.plot_forest(idata.posterior.sales_pred[:, :, [0, 19]],
               kind='ridgeplot',
               combined=True,
               hdi_prob=0.6,
               ridgeplot_truncate=False,
               ridgeplot_quantiles=[.005, .5, .995],
               ridgeplot_overlap=0.5,
               colors='lightblue',
               figsize=(6, 4));

【実行結果】
青い領域が 60 % HDI、ひし形が 0.5%点、50%点、99.5%点（両端で 99% HDI を構成）です。
11 度の上側 50% を 30 度の下側 50% が被っている感じがします。

④ 散布図と事後予測の重ね描き
単回帰分析のときに描いた散布図＋95%予測区間に類似するチャートを描きます。
事後予測の HDI の塗りつぶしには arviz の plot_hdi を利用します。

# 予測分布の区間表示

# 描画領域の設定
plt.figure(figsize=(8, 4))
# 観測値の散布図の描画
# 事後予測サンプルの中央値の描画
plt.plot(temperature_pred,
         idata.posterior.sales_pred.stack(sample=('chain', 'draw')).median(axis=1),
         color='tab:red')
# 事後予測サンプルの95% HDI の塗りつぶし
az.plot_hdi(temperature_pred, idata.posterior.sales_pred, hdi_prob=0.95,
            fill_kwargs={'color': 'lightpink', 'alpha': 0.3})
# 事後予測サンプルの50% HDI の塗りつぶし
az.plot_hdi(temperature_pred, idata.posterior.sales_pred, hdi_prob=0.50,
            fill_kwargs={'color': 'lightpink', 'alpha': 0.6})
# 観測値の散布図の描画
sns.scatterplot(data=file_beer_sales_2, x='temperature', y='sales', legend=True)
# 修飾
plt.title('ビールの売上と気温の関係：事後予測の中央値と50% および 95% HDI', loc='left')
plt.xticks(range(10, 31, 2));

【実行結果】
赤実線が事後予測の中央値、濃赤色の塗りつぶしが事後予測の 50% HDI、薄赤色の塗りつぶしが事後予測の 95% HDI です。

（参考：単回帰分析の予測区間）

単回帰分析の予測区間は直線的ですが、ベイズ単回帰モデルの方は観測値のデータ点の「幅」（縦方向の散らばり度合い）によって、HDI の「幅」が異なっています。

以上で ベイズ流の単回帰モデル（後編）を終わりにします。
面白かったですね。

今回の記事は以上です。
楽しかったですね！

シリーズの記事

次の記事

前の記事

PyMC版

目次

ブログの紹介

note で８つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．統計・データ分析とつながる

シリーズ「統計・データ分析とつながる」は、統計・データ分析との「つながり」を発掘して、コラム風に仕立てたブログシリーズです。
生成 AI の力を借りながら、統計・データ分析の入り口をイメージして、自由気ままに書きました。
たとえば…
・日常生活と統計のつながり
・統計検定２級からその先へのつながり
気楽にお読みいただけたら嬉しいです🍀

３．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

４．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

５．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

６．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

７．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

８．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。