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

書籍の著者　馬場真哉 先生

この記事は、書籍「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'

第２章 単回帰モデル

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

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

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

# ファイルの読み込み
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));

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

相関係数を確認しましょう。

# 相関係数
file_beer_sales_2.corr().round(2)

【実行結果】
相関係数は 0.66。
気温と売り上げとの間には正の相関があります。
単回帰モデルの回帰直線は右上がりの直線になりそうです。

おまけの可視化を１つ。
散布図の外側に１変数のヒストグラムを同時描画します。
seaborn の jointplot を利用します。

# ヒストグラムと散布図の同時描画

# 散布図の描画
sns.jointplot(data=file_beer_sales_2, x='temperature', y='sales', ratio=2)
plt.xticks(range(10, 31, 2));

【実行結果】

上部に配置した気温のヒストグラムを見ると、一様に分布している感じです。
右部に配置した売り上げのヒストグラムを見ると、平均値あたりを頂点とする山型 $${\approx}$$ 正規分布のように見えます。
目的変数 sales は正規分布に従うと仮定…、ただの独り言です…

🔵🔵🔵

単回帰分析

まず単回帰分析を確認しておき、あとでベイズ統計モデルの結果と比べてみましょう。
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()

【実行結果】
こちらは（お馴染みの？）単回帰分析のサマリーです。
決定係数（R-squared）は 0.434 であり、当てはまりがいい、とは言い切れない感じです。
F 値（F-statistic）は 75.15 と大きな値であり、モデルは有効な感じはします。

係数の推定値にフォーカスしてみます。

# 要約表から係数の推定値の部分を取り出し
res_sm.summary().tables[1]

【実行結果】

Intercept：切片の推定値（coef）は 21.17、temperature：傾きの推定値は 2.46 です。
傾きに注目すると、p 値（P>|t|）がとても小さく、また、95% 信頼区間（[0.025  0.975]）は０を含んでいないので、傾きは０ではないと考えられ、気温は売り上げへの影響力がある、と解釈できます。

傾きの推定値を解釈しましょう。
「温度が１度高くなると売り上げは 2.46 万円増える」

誤差 $${\varepsilon}$$ の標準偏差を確認しましょう。
単回帰分析の結果 res_sm の scale 属性で誤差分散を取り出して正の平方根をとります。

# 誤差の標準偏差
print('誤差の標準偏差:', np.sqrt(res_sm.scale))

【実行結果】
誤差の標準偏差は 16.85 です。

最後にデータの散布図と回帰直線を重ね描きしましょう。
単回帰分析の結果 res_sm に対して predict メソッドを適用して予測値を得ます。

# 散布図と回帰直線の描画

## 予測
# 予測に用いるx:気温データの作成
x_vals = np.linspace(*np.sort(file_beer_sales_2.temperature)[[0, -1]], 100)
# 回帰分析の結果を用いて予測を実行
y_preds = res_sm.predict(dict(temperature=x_vals))

## 描画
# 描画領域の設定
plt.figure(figsize=(8, 4))
# 回帰直線の描画
plt.plot(x_vals, y_preds, color='tab:red')
# 観測値の散布図の描画
sns.scatterplot(data=file_beer_sales_2, x='temperature', y='sales')
# 修飾
a, b = res_sm.params
plt.title(f'ビールの売上と気温の関係：回帰直線 y = {a:.3f} + {b:.3f}x', loc='left')
plt.xticks(range(10, 31, 2));

【実行結果】
回帰直線は右上がりの形状をしています。
グラフタイトルに回帰式を付記しています。

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

🔵🔵🔵

ベイズモデリング by Stan

テキスト 2.4 節、2.5 節に相当します。

① モデルの概要
次のモデルを実装します。

目的変数 $${sales}$$ は正規分布に従うと仮定しています（どこかで見たフレーズが…）。
パラメータ $${Intercept,\ \beta,\ \sigma}$$ には事前分布を明示的に設定しません。

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

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

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

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

【実行結果】なし

③ データの整理
Stan に渡すデータセットを辞書にまとめます。

# p.170 データセットの準備

# 辞書にまとめる
data_dict = dict(
    N=len(file_beer_sales_2),
    sales=file_beer_sales_2['sales'],
    temperature=file_beer_sales_2['temperature'],
)

【実行結果】なし

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

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

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

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

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

MCMC の準備が整いました！

🔵🔵🔵

MCMC の実行

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

%%time
# p.170 乱数の生成
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.170 結果はこちら（要約統計量の表示）
az.summary(idata, hdi_prob=0.95)

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

【考察】
パラメータの事後分布の平均（mean）を単回帰分析の推定値と比べると、両者は近い値になっています。

（参考：単回帰分析の係数の推定値、誤差の標準偏差の推定値）

単回帰分析と同様に係数の事後平均を解釈すると…
「気温が１度高くなると売り上げは 2.46 万円増える」

🔵🔵🔵

事後分布の可視化

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

① トレースプロットの描画
トレースプロットを描画します。
テキスト 図 3.2.2 に相当します。

# p.173 図3.2.2 事後分布とトレースプロットの図示
az.plot_trace(idata, compact=False,
              backend_kwargs={'tight_layout': True});

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

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

② 事後分布プロットの活用例
パラメータの事後分布を少々掘り下げてみましょう。
「単回帰分析」では難しいけども「ベイズ統計による単回帰モデル」では容易に解釈できることを試してみます。

事後分布プロットを活用して、パラメータの事後分布に関する基準値の下の確率、上の確率を可視化しましょう。
arviz の plot_posterior を利用します。

引数 ref_val で調べたい「基準値」を設定します。
今回は intercept：11、beta：1.5、sigma：15 を基準値にします。

# 事後分布プロット
az.plot_posterior(
    idata, hdi_prob=0.95, textsize=20, round_to=3,
    ref_val={'Intercept': [{'ref_val': 11}],
             'beta': [{'ref_val': 1.5}],
             'sigma': [{'ref_val': 15}]}
);

【実行結果】
オレンジの垂直線が「基準値」、オレンジ文字が下の確率・上の確率です。

【見方】
ひとまず一覧表形式でまとめます。

例えば傾き beta が 1.5 を下回る確率は $${0\%}$$ です。
このベイズ統計モデルが適切ならば、傾きが 1.5 を下回ることは無い ⇒ 1.5  を下回るリスクは０と見込まれます。

ベイズ統計では「パラメータは確率変数」と扱えるので、MCMC サンプルに基づく事後分布を用いて、パラメータの確率を見積もることができるのです。
一方で通常の（最小二乗法による）単回帰分析の場合、パラメータは「定数」扱いなので、ベイズ統計モデルのような「確率」を語ることが困難なのです。

③ 回帰直線の可視化
最後にデータの散布図と回帰直線を重ね描きします。
事後平均を用いて回帰直線の切片と傾きとします。

# 散布図と回帰直線の描画

## 予測
# MCMC サンプルの平均値を算出して切片a、傾きbを取得
a = idata.posterior.Intercept.mean().data
b = idata.posterior.beta.mean().data
# yの予測値を算出
y_hats = a + b * x_vals

## 描画
# 描画領域の設定
plt.figure(figsize=(8, 4))
# 回帰直線の描画
plt.plot(x_vals, y_hats, color='tab:red')
# 観測値の散布図の描画
sns.scatterplot(data=file_beer_sales_2, x='temperature', y='sales')
# 修飾
plt.title(f'ビールの売上と気温の関係：回帰直線 y = {a:.3f} + {b:.3f}x', loc='left')
plt.xticks(range(10, 31, 2));

【実行結果】

このコードでは Stan の範囲外で目的変数の予測値を計算しました。
次回記事では Stan を用いて、つまりベイズ統計モデリングの範疇で、目的変数の予測を行います。

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

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

シリーズの記事

次の記事

前の記事

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の教科書です。
よかったらぜひ、お試しくださいませ。

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