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

書籍の著者　馬場真哉 先生

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

書籍の第３部は「一般化線形モデル」のベイズ統計モデリングです。
「単回帰モデル」から始めます。
今回・次回記事で単回帰モデルの構築と予測に取り組みます。

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

はじめに

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

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

準備

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

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

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

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

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

# インポート

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

# ベイズ統計モデリング
import pymc as pm                        # pymc
import arviz as az                       # 分析・可視化

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

# 可視化
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 PyMC

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

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

目的変数 $${sales}$$ は正規分布に従うと仮定しています（どこかで見たフレーズが…）。
正規分布の平均パラメータ $${\mu}$$ は、切片 $${intercept}$$ と傾き $${\beta}$$ の単回帰式を示す変数です。
標準偏差パラメータ $${\sigma}$$ は標準偏差 100000 の半正規分布に従うとします。
単回帰式の切片・傾きパラメータ $${intercept,\ \beta}$$ は平均０、標準偏差 100000 の正規分布に従うとします。

② PyMC のモデル設定
上式のモデルを PyMC で記述します。

# モデリング

# coordsの設定
coords = {'id': file_beer_sales_2.index.values}

# モデルの定義
with pm.Model(coords=coords) as model:
    
    ## dataの設定
    # 目的変数: 売上データ
    Y = pm.Data('Y', value=file_beer_sales_2['sales'].values, dims='id')
    # 説明変数: 気温データ
    temperature = pm.Data(
        'temperature',value=file_beer_sales_2['temperature'].values, dims='id')

    ## 事前分布: 無情報事前分布的な分布
    # 切片
    intercept = pm.Normal('intercept', mu=0, sigma=1e5)
    # 係数
    beta = pm.Normal('beta', mu=0, sigma=1e5)
    # 標準偏差
    sigma = pm.HalfNormal('sigma', sigma=1e5)

    ## 線形予測子（式を見やすくする目的で別出し、MCMCサンプルには残らない）
    mu = intercept + beta * temperature

    ## 尤度関数: 正規分布を仮定
    obs = pm.Normal('obs', mu=mu, sigma=sigma, observed=Y, dims='id')

【実行結果】なし

【コードの補足】
このモデルの変数 mu の書き方をすると、PyMC内部変数が定義されないため、確率変数でも決定論的変数でもない扱いになります。
したがって、mu の MCMCサンプルは生成されません。
（次回記事では別の扱いをしますので、お楽しみに🍀）

モデルの数式ライク表示とグラフィカルモデル描画をします。

# モデルの表示
model

【実行結果】
最下行の obs が尤度、上の３行が事前分布です。

# モデルの可視化
pm.model_to_graphviz(model)

【実行結果】
説明変数（Data）の temperature、尤度 obs、目的変数（Data）の Y は、軸 id を設定した 100 個のデータを持っています。
事後分布を推定するパラメータ intercept、beta、sigma は（mu を経由して）尤度 obs のパラメータになっています。

🚀🚀🚀

MCMC の実行

テキスト 2.6 節に相当します。
MCMCを実行しましょう。
NUTS サンプラーに nutpie を利用します。

%%time
# p.170 MCMCの実行
with model:
    idata = pm.sample(draws=1000, tune=1000, chains=4, random_seed=1,
                      nuts_sampler='nutpie')

【実行結果】
Divergences（ダイバージェンス）は０件です。

収束の確認をします。
MCMC サンプルの要約表を表示します。

# p.170 結果はこちら（要約統計量の表示）
var_names = ['intercept', 'beta', 'sigma']
az.summary(idata, var_names=var_names, 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, var_names=var_names, compact=False,
              backend_kwargs={'tight_layout': True});

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

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

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

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

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

# 事後分布プロット
az.plot_posterior(
    idata, var_names=var_names, 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));

【実行結果】

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

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

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

シリーズの記事

次の記事

前の記事

Stan版

目次

ブログの紹介

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

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