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

ネイピア DS

2026年4月3日 11:33

書籍の著者　馬場真哉先生

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

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

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

はじめに

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

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

準備

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

3.2 分析の準備
3.4 予測のためのデータの整理
3.5 予測のためのStanファイルの修正
3.6 MCMCの実行
3.7 予測分布の可視化

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

■ コード記述法
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'

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

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

テキスト 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 を利用して、次の単回帰モデルを実装します！

$$
sales = 切片 + 傾き \times temperature + \varepsilon
$$

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

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

$$
\mathtt{sales} \sim \mathtt{temperature}
$$

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

# 単回帰分析 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_vals の作成 ※ shape=(100,)
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 点に対する予測値たちです。

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

$$
\begin{array}{ll}
列名 & 内容 \\
\hline
\text{mean} & 予測値平均（点推定） \\
\text{mean\_se} & 予測値平均の標準誤差 \\
\text{mean\_ci\_lower} & 予測値平均の95\%信頼区間の下限値 \\
\text{mean\_ci\_upper} & 予測値平均の95\%信頼区間の上限値 \\
\text{obs\_ci\_lower} & 個別観測値の95\%予測区間（※）の下限値 \\
\text{obs\_ci\_upper} & 個別観測値の95\%予測区間（※）の上限値 \\
\end{array}
$$

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

データの散布図と予測値を重ね描きしましょう。
予測値平均 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 PyMC

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

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

$$
\begin{align*}
sales &\sim \text{Normal}\ (\mu,\ \sigma^2) \\
\mu &= intercept + \beta \cdot temperature \\
intercept &\sim \text{Normal}\ (0,\ (1e5)^2) \\
\beta &\sim \text{Normal}\ (0,\ (1e5)^2) \\
\sigma &\sim \text{HalfNormal}\ ((1e5)^2) \\
\end{align*}
$$

目的変数 $${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
    mu = pm.Deterministic('mu', intercept + beta * temperature, dims='id')

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

【実行結果】なし

【コードの補足】
変数 mu について前回記事から変更し、決定論的変数にしています。
したがって、mu の MCMCサンプルは生成されます。

【テキストとの相違】
Stan は generated quantities ブロックに予測のためのコードが必要です。
しかし、PyMC はモデルの外側で事後予測サンプルを生成できます。
したがって、上のPyMC モデルには予測値を持たせないようにしています。

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

# モデルの表示
model

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

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

【実行結果】
説明変数（Data）の temperature、尤度 obs、目的変数（Data）の Y は、軸 id を設定した 100 個のデータを持っています。
パラメータ intercept、beta は線形予測子の形式で mu のパラメータになっています。
mu と sigma は尤度 obs のパラメータになっています。

🚀🚀🚀

MCMC の実行

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

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

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

収束の確認をします。
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.4 節、3.7 節に相当します。

① 予測用の説明変数データの作成
予測に用いる説明変数 temperature の値を設定します。

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

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

② 事後予測サンプルの生成
今回記事のクライマックスです（当社比）。
PyMC の sample_posterior_predictive 関数で観測値 obs と mu の事後予測サンプルを生成します。

# 事後予測サンプリング

# 事後予測サンプルの取得
with model:
    # temperatureにx_valsを設定
    pm.set_data({'temperature': temperature_pred}, coords={'id': temperature_pred})
    # Yの要素数を揃えるためダミーを設定
    pm.set_data({'Y': temperature_pred})
    # 事後予測サンプリング
    ypred = pm.sample_posterior_predictive(
        idata, var_names=['obs', 'mu'], random_seed=0)

【実行結果】

【コードの補足】set_dataを用いた事後予測のコード

set_data
モデル内の "Data" 変数の値を変更する関数であり、事後予測サンプリングなどで活用します。
本ケースではまず説明変数 temperature の値を変えます。
同時に、coords の id の値も変えます。要素数が変わるからです。
そして半ば仕方なく、同じ id 軸をもたせた目的変数 Y にも値設定します。意味のないダミー値で埋めます。
sample_posterior_predictive
モデル内の変数 obs と mu を対象にして、事後予測サンプリングを実行します。
生成したサンプルは ypred に格納します。

このコードの選択にあたっては、以下の PyMC コミュニティフォーラムの議論を参照しました。
ありがとうございます！

事後予測サンプル ypred の中身を確認しましょう。

# 事後予測サンプルを確認
ypred

【実行結果】
グループ posterior_predictive（事後予測）に事後予測サンプルが格納されています。
形状は Dimensions に記載の (chain=4, draw=1000, id=21) の３次元です。
気温１度あたり 4000 サンプルあります。

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

# p.178 図3.3.1 95%予測区間 太い線は50%区間 ※インデックスは気温
az.plot_forest(ypred.posterior_predictive, var_names=['obs'], combined=True,
               hdi_prob=0.95);

【実行結果】
[ ] の中に気温の値が表示されます！ヽ(=´▽`=)ﾉﾜｰｲ
coords：id 軸の値に気温をセットしたからです。

④ 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(ypred.posterior_predictive, var_names=['mu', 'obs'],
                    coords={'id': [11]}, combined=True, hdi_prob=0.95)
ax[0].set_title('95%HDI：11度のケース');

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

【コードの補足】
引数「$${\mathtt{coords=\{'id': [11]\}}}$$」によって、id 軸の値が 11（つまり 11 度）のデータを抽出しています。

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

# p.179 図3.3.3 予測分布の図示
az.plot_forest(ypred.posterior_predictive.obs.sel(id=[11, 30]),
               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 を利用します。

# 事後予測の中央値と50%・95% HDI の描画

# 描画領域の設定
plt.figure(figsize=(8, 4))
# 事後予測サンプルの中央値の描画
plt.plot(temperature_pred,
         ypred.posterior_predictive.obs.median(dim=('chain', 'draw')),
         color='tab:red')
# 事後予測サンプルの95% HDIを描画
az.plot_hdi(temperature_pred, ypred.posterior_predictive.obs, hdi_prob=0.95,
            fill_kwargs={'color': 'lightpink', 'alpha': 0.3})
# 事後予測サンプルの50% HDIを描画
az.plot_hdi(temperature_pred, ypred.posterior_predictive.obs, 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 の「幅」が異なっています。

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

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

シリーズの記事

次の記事

前の記事

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

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

いいなと思ったら応援しよう！

応援ありがとうございます。これからもがんばって記事を作成します！