「ベイズ統計モデリングによるデータ分析入門」をPythonとStanで写経 ～ Vol.10 デザイン行列を用いた単回帰モデルの推定
書籍の著者 馬場真哉 先生
この記事は、書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」第３部第４章「デザイン行列を用いた一般化線形モデルの推定」の Python 写経活動記録です。
書籍の第３部は「一般化線形モデル」のベイズ統計モデリングです。
今回も引き続き「単回帰モデル」を堪能します。
書籍によると「formula構文」と「デザイン行列による表現」に慣れる章です。
では書籍を開いてベイズ統計モデリングの旅に出かけましょう🚀
はじめに
このブログシリーズは書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」（講談社、「テキスト」と呼びます）の Python 写経です。
テキストの紹介と引用表記はリンク先の記事に掲載しています。
準備
■ 記事の範囲
この記事はテキスト第３部第４章の以下の節を取り扱います。
4.2 分析の準備
4.4 fromula構文を用いたデザイン行列の作成
4.5 デザイン行列を使うためのStanファイルの修正
4.6 MCMCの実行
ベイズ統計モデリングの理論面が気になった場合には、ぜひテキストの第１部、第３部第１章をご覧いただき、本記事との繋がりをご確認下さいませ。
■ コード記述法
Jupyter Notebook 形式でコードを記述します。
■ 利用データ
テキスト・サポートサイトのデータファイルを引用しています。
Jupyter Notebook ファイルと同一フォルダ内の「data」フォルダにデータファイルを格納しています。
■ ライブラリのインポート
この記事で用いるライブラリをインポートします。
# インポート
# 数値計算
import numpy as np
import pandas as pd
# ベイズ統計モデリング
from cmdstanpy import CmdStanModel # stan
import arviz as az # 分析・可視化
# デザイン行列
from patsy import dmatrices
# 統計処理
import scipy.stats as stats # 正規分布乱数の生成
# ユーティリティ
import os
# 可視化
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_theme() # ggplot風のスタイル
plt.rcParams['font.family'] = 'Meiryo'第４章 デザイン行列を用いた一般化線形モデルの推定
データの読み込み
テキスト 4.2 節に相当します。
テキストの仮想のビールの売り上げデータを引用いたします。
前々回記事から利用しているデータです。
csv ファイルを pandas のデータフレーム形式で変数 file_beer_sales_2 に読み込みます。
# p.180 分析対象データの読み込み
# ファイルの読み込み
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 は気温（おそらく摂氏℃）です。
気温とビール売り上げの関係を単回帰モデルで分析します。
🔵🔵🔵
デザイン行列と係数ベクトル
前回記事までの単回帰モデルでは、線形予測子 $${\mu}$$ の切片と傾きに別々の係数を当てて表現していました。
$$
\mu_i = \underbrace{\beta_0}_{切片} + \underbrace{\beta_1}_{傾き} \underbrace{x_i}_{説明変数}
$$
上式は識別子 $${i}$$ 単位であり、標本サイズ $${n}$$ 個の数式が並びます。
これをベクトル・行列でひとまとめにして表現してみます。
$$
\underbrace{\begin{bmatrix} \mu_1 \\ \mu_2 \\ \vdots \\ \mu_n \end{bmatrix}}_{\bm \mu}
= \begin{bmatrix} \beta_0 + \beta_1 x_1 \\ \beta_0 + \beta_1 x_2 \\ \vdots \\ \beta_0 + \beta_1 x_n\end{bmatrix}
= \underbrace{\begin{bmatrix} 1 & x_1 \\ 1 & x_2 \\ \vdots \\ 1 & x_n \end{bmatrix}}_{\bm X}
\underbrace{\begin{bmatrix} \beta_0 \\ \beta_1 \end{bmatrix} }_{\bm \beta}
= \bm X \bm \beta = \bm \mu \\
$$
$${\bm X}$$ は説明変数の行列であり、書籍は「デザイン行列」と呼んでいます。
１列目の $${1}$$ のみの列ベクトルは定数項です。
切片 $${\beta_0}$$ と掛け算する対象です。２列目は 説明変数ベクトル $${\bm x}$$ です。
単回帰モデルでは説明変数が１個なので１列＝ベクトルですが、重回帰モデルなどの説明変数が複数存在するモデルでは説明変数が複数列＝行列になって $${\bm X}$$ に合流します。
$${\bm \beta}$$ は係数ベクトルです。
単回帰モデルでは切片 $${\beta_0}$$ と傾き $${\beta_1}$$ の２要素で構成されますが、重回帰モデルなどの説明変数が複数存在するモデルでは説明変数の数 $${k}$$ 個の係数要素が合流します。
🔵🔵🔵
デザイン行列とformula構文
テキスト 4.4 節に相当します。
テキストにならい、formula 構文を活用してデザイン行列を作成してみます。
formula 構文のオリジナルは R 言語（だと思うの）です。
Python には R 言語ライクに formula 構文を使ってデザイン行列を作成するライブラリがいくつかあります。
ここでは patsy ライブラリを利用します。
pandas データフレームの列名（＝変数名）を用いて、目的変数と説明変数等の関係を formula 構文 で表現し、デザイン行列等を作成できます。
patsy 公式サイトの Quick Startページのリンクはこちら
今回の例題データは 目的変数：sales、説明変数：temperature であり、formula 構文で単回帰モデルを記述すると：
$$
\mathtt{sales} \sim \mathtt{temperature}
$$
です。
チルダ $${\sim}$$ の左辺に目的変数、右辺に説明変数を書きます。
上式の書き方の場合、定数項を含む formula と同じ意味合いになります。
$$
\underbrace{\mathtt{sales}}_{目的変数} \sim \underbrace{\mathtt{1}}_{定数項} + \underbrace{\mathtt{temperature}}_{説明変数}
$$
実はこの formula 構文は、前々回・前回記事の単回帰分析の際にこっそり使っています。
statsmodels は裏側で patsy ライブラリを用いて formula 構文 ⇒ デザイン行列作成の機能を実現しています。
formula を定義して、デザイン行列を作成しましょう。
patsy の dmatrices 関数を利用します。
目的変数とデザイン行列を生成して、先頭５行を表示します。
# p.182 formulaの作成
# formula構文の設定
formula_lm = 'sales ~ temperature'
# 目的変数Y, デザイン行列（説明変数）Xの作成
Y, X = dmatrices(formula_lm, file_beer_sales_2, return_type='dataframe')
# デザイン行列の先頭5行の表示
X.head()【実行結果】
こちらはデザイン行列 $${\bm X}$$ を体現したデータです。
定数項の変数名が「Intercept」（切片）になっています。
# 目的変数の先頭5行の表示
Y.head()【実行結果】
こちらは目的変数のデータです。
【コード補足】
こちらは formula の定義です。変数 formula_lm に設定しています。
# formula構文の設定
formula_lm = 'sales ~ temperature'こちらは dmatrices 関数に formula を与えて、目的変数 Y_dm、デザイン行列 X_dm を作成しています。
３つの引数は、①formula、②データフレーム名、③戻り値をデータフレームで返す設定、です。
# 目的変数Y, デザイン行列（説明変数）Xの作成
Y_dm, X_dm = dmatrices(formula_lm, file_beer_sales_2, return_type='dataframe')🔵🔵🔵
ベイズモデリング by Stan
テキスト 4.5 節に相当します。
① モデルの概要
デザイン行列と係数ベクトルを用いて、次のモデルを実装します。
$$
\begin{align*}
Y &\sim \text{Normal}\ (\bm \mu,\ \sigma^2) \\
\bm \mu &= \bm X \bm \beta \\
\end{align*}
$$
目的変数 $${Y}$$（sales）は正規分布に従うと仮定しています。
正規分布の平均パラメータ $${\bm \mu}$$ は、デザイン行列 $${\bm X}$$ と係数ベクトル $${\bm \beta}$$ の積です。
係数ベクトル $${\bm \beta}$$、標準偏差パラメータ $${\sigma}$$ には事前分布を明示的に設定しません。
② Stan のモデル設定
Stan ファイル（Stan コード）を作成します。
テキストの Stan ファイル名 および Stan コードを引用いたします。
📑ファイル名：3-4-1-lm-design-matrix.stan
data {
int N; // 標本サイズ
int K; // デザイン行列の列数（説明変数の数＋１）
vector[N] Y; // 目的変数
matrix[N, K] X; // 説明変数：デザイン行列
}
parameters {
vector[K] b; // 切片を含む係数ベクトル
real<lower=0> sigma; // 標準偏差
}
model {
vector[N] mu = X * b; // デザイン行列と係数ベクトルの積
Y ~ normal(mu, sigma);
}【実行結果】なし
【コードの補足】
線形予測子 $${\mathtt{mu}}$$ は デザイン行列と係数ベクトルの積です。
vector[N] mu = X * b; // デザイン行列と係数ベクトルの積【Stan コードでデザイン行列を用いるメリット】
テキスト p.184 では Stan の外から与える「デザイン行列を変えるだけで、さまざまなモデルに対応できるので、Stan ファイルを書き換える必要がなくなります」と、デザイン行列を用いるメリットを説明しています。
③ データセットの作成
標本サイズ・説明変数の数を算出して、Stan に渡すデータセットを辞書にまとめます。
# p.183 データセットの準備
# サンプルサイズ、デザイン行列の列数（説明変数の数＋１）
N, K = X.shape
# 辞書にまとめる ※Y:pd.Series, X:pd.DataFrame
data_dict_design = dict(N=N, K=K, Y=Y['sales'], X=X)【実行結果】なし
④ モデルのコンパイル
モデルのコンパイルを実行します。
%%time
# モデルのコンパイル
# stanプログラムファイルのパス指定
stan_file = '3-4-1-lm-design-matrix.stan'
current_dir = os.path.abspath(os.getcwd())
stan_path = os.path.join(current_dir, 'stan', stan_file)
# モデルオブジェクトの作成（exeの作成）
model = CmdStanModel(stan_file=stan_path) # stanのpathを設定【実行結果】（右側のファイルパスは記載省略）
MCMC の準備が整いました！
🔵🔵🔵
MCMC の実行
テキスト 4.6 節に相当します。
MCMCを実行しましょう。
%%time
# p.183 MCMCの実行
fit = model.sample(
data=data_dict_design, # 対象データ
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.184 要約統計量の表示
az.summary(idata, hdi_prob=0.95)【実行結果】
$${\widehat{R}}$$（r_hat）、有効サンプル数（ess_bulk、ess_tail）に問題はなさそうです。
係数ベクトル b には係数のインデックス $${0,\ 1}$$ が付記されていています。
b[0] は切片、b[1] は傾きです。
前回記事の単回帰モデルの要約表と比べてみます。
（参考：前回記事の単回帰モデルの要約表）
事後平均（mean）、標準誤差（sd）、95% HDI（hdi_2.5%、hdi_97.5%）の値は、若干異なる結果になっています。
結果が変わる原因を Gemini に訊いてみました。
長文ですがぜひお読み下さい！
Stanにおける単回帰モデルの定義方法（個別定義 vs デザイン行列によるベクトル定義）で要約統計量に差が出る理由は、モデルの数学的意味は同じでも、計算アルゴリズムの効率と数値処理のプロセスが異なるためです。
2026年現在のStanの仕様に基づき、主な要因を整理します。
ベクトル化による勾配計算の効率化
StanはHMC（ハミルトニアンモンテカルロ）およびその拡張であるNUTSアルゴリズムを使用しており、対数尤度の「勾配（微分値）」を計算してパラメータ空間を探索します。ベクトル定義：$${\mathtt{y} \sim \mathtt{normal(X * beta, sigma)}}$$ のようにデザイン行列を用いると、Stan内部でベクトル化された演算が行われます。これにより、個別に計算するよりも計算誤差が蓄積しにくく、より滑らかな勾配が得られます。
個別定義：$${\mathtt{y} \sim \mathtt{normal(b0 + b1 * x, sigma)}}$$ では、スカラ演算が繰り返されます。これが計算速度だけでなく、サンプラーが提案する「次の地点」の採択率や探索の軌道に微小な影響を与えます。
有効サンプルサイズ（ESS）と収束の質
統計量の差の多くは、サンプリングの効率（混合の良さ）に起因します。回帰係数 $${b_{0}}$$ と $${b_{1}}$$ は通常、強い相関（共分散）を持ちます。ベクトルとして一括定義されたパラメータは、Stanのサンプラーがその相関構造を「一塊の幾何学的構造」としてより効率的に把握しやすくなる場合があります。
効率が向上すると有効サンプルサイズ（ESS）が大きくなり、結果として平均値の推定精度（標準誤差、MCSE）が向上します。ESSが異なれば、有限のサンプルから得られる要約統計量には必ず微小な差異が生じます。
初期値と乱数の影響
MCMCは乱数に基づくシミュレーションです。定義方法を変えると、内部的なパラメータの並び順やメモリ配置が変わり、同じ乱数シードを使用してもサンプラーが辿る経路が完全に一致しなくなります。
特に収束が不十分な場合や、事後分布の形状が複雑な場合、この経路の差が平均や標準偏差の数値的な「揺らぎ」として現れます。
数値的安定性
デザイン行列を用いた行列演算は、IntelのAVXなどのベクトル命令セットを活用した高度な最適化が行われます。スカラ演算を繰り返す形式に比べ、浮動小数点演算の丸め誤差が抑えられる傾向にあり、これが対数尤度の計算値に極めて僅かな差を生み、サンプリング結果に反映されることがあります。
まとめ
十分な反復回数（Iter）を確保し、収束診断指標（$${\widehat R}$$）が1.01未満であれば、これらの差は統計的な誤差範囲内であり、モデルの本質的な結論に影響を与えるものではありません。より安定した結果を得るには、計算効率に優れるデザイン行列を用いたベクトル化表現が推奨されます。
差が生じるのは仕方がないことのようです。
大きな差が無いことを確認したので、次に進みましょう🚀
トレースプロットを描画します。
# トレースプロットの描画
az.plot_trace(idata, compact=False, backend_kwargs={'tight_layout': True});【実行結果】
右側のチャートがゲジゲジしています。
以上のチェックに基づいて、収束していると考えましょう。
テキスト 第４章 はここまでとなります。
🔵🔵🔵
回帰直線の可視化
実は次の第５章をスキップします…
Python に brms（裏で Stan が動く）に相当するライブラリが見つからないこと、および、分析データはすでに第２・３・４章で取り組み済みであることより、泣く泣くのスキップです。
しかし、第５章最後のチャートを未だ描いていないので、せっかくですから、こちらで描いてしまいましょう。
回帰直線＋95%信用区間、回帰直線＋95%予測区間の２つのチャート描画に取り組みます。
テキスト p.199 ～ p.200 に相当します。
描画に利用するデータは MCMC サンプルを PyMC の外で加工して作成します。
① 回帰直線＋95%信用区間
単回帰モデルの予測平均（＝線形予測子）$${\mu}$$ をMCMCサンプルから算出して、回帰直線と95%信用区間を描画します。
# p.199 図3.5.3 回帰直線95%信用区間付き
## データの作成
# MCMCサンプルからβ0, β1, σを取り出し
beta0s, beta1s, sigmas = fit.draws_pd().iloc[:, -3:].values.T
# x軸の値の設定
x_val = np.linspace(X.iloc[:, 1].min(), X.iloc[:, 1].max(), 1001)
# μの算出
mu_preds = (beta0s + np.outer(x_val, beta1s)).T
# μの95%信用区間の算出
mu_95ci = np.quantile(mu_preds, q=[0.025, 0.975], axis=0).T
# 描画領域の設定
plt.figure(figsize=(8, 4))
# 観測値の散布図の描画
plt.scatter(X['temperature'], Y, color='navy', s=10, label='観測値')
# サンプルの平均値の描画
plt.plot(x_val, mu_preds.mean(axis=0), label='回帰直線')
# サンプルの95%信用区間の塗りつぶし
plt.fill_between(x_val, mu_95ci[:, 0], mu_95ci[:, 1], alpha=0.2,
label='95%信用区間')
# 修飾
plt.xlabel('温度 [℃]', fontsize=12)
plt.ylabel('売上', fontsize=12)
plt.xticks(range(10, 31, 2))
plt.legend();【実行結果】
予測平均 $${\mu}$$ の 95%信用区間は、予測値の平均（期待値）が 95% の確率で収まる区間です。
かなり狭い範囲なんだと感じます。
【コードの補足】
◆ コード１：パラメータの事後分布からのサンプルを取得
# MCMCサンプルからβ0, β1, σを取り出し
beta0s, beta1s, sigmas = fit.draws_pd().iloc[:, -3:].values.Tこちらのコードは Stan の MCMC サンプルを格納した変数 fit に対して draws_pd メソッドを適用して、pandas データフレーム形式の MCMC サンプルを取得します。
このデータフレームから $${\beta_0\, \beta_1,\ \sigma}$$ の MCMC サンプルを取り出します。
◆ コード２：予測平均 $${\mu}$$（＝線形予測子）の算出
# x軸の値の設定
x_val = np.linspace(X.iloc[:, 1].min(), X.iloc[:, 1].max(), 1001)
# μの算出
mu_preds = (beta0s + np.outer(x_val, beta1s)).Tこちらのコードは $${\mu = \beta_0 + \beta_1 x}$$ を計算しています。
◆ コード３：$${\mu}$$ の 95% 信用区間の算出
# μの95%信用区間の算出
mu_95ci = np.quantile(mu_preds, q=[0.025, 0.975], axis=0).Tコード２で計算した $${\mu}$$（shape=(4000, 1001)）のサンプルから 2.5% 点と 95.5% 点を算出します。
🔵
② 回帰直線＋95%予測区間
$${\mu}$$ の（MCMC的）サンプルに観測値の「ばらつき」（正規分布乱数）を加味して「予測分布からのサンプルに相当するデータ」を作成し、回帰直線と95%予測区間を描画します。
# p.200 図3.5.4 回帰直線95%予測区間付き
## MCMCサンプルデータからyの予測値を算出
# yの予測値の算出 ※scipy.statsで正規分布乱数を付加
y_preds = (mu_preds.T + stats.norm.rvs(loc=0, scale=sigmas, random_state=123)).T
# yの予測値の95%予測区間の算出
y_95ci = np.quantile(y_preds, q=[0.025, 0.975], axis=0).T
# 描画領域の設定
plt.figure(figsize=(8, 4))
# 観測値の散布図の描画
plt.scatter(X['temperature'], Y, color='navy', s=10, label='観測値')
# yの予測値の平均値の描画
plt.plot(x_val, y_preds.mean(axis=0), label='回帰直線')
# yの95%予測区間の塗りつぶし
plt.fill_between(x_val, y_95ci[:, 0], y_95ci[:, 1], color='lightpink', alpha=0.3,
label='95%予測区間')
# 修飾
plt.xlabel('温度 [℃]', fontsize=12)
plt.ylabel('売上', fontsize=12)
plt.xticks(range(10, 31, 2))
plt.legend();【実行結果】
ばらつきが加わるので 95% 予測区間の幅はかなり広くなります。
前回記事の事後予測 95% HDI とよく似ています（当たり前）。
今回記事の方が縁が滑らかな理由は、単に x 軸の値の粒度が細かいからです。
（参考：前回記事の事後予測 95% HDI）
【コードの補足】
◆ コード：y の予測値（観測値ベース）の算出
# yの予測値の算出 ※scipy.statsで正規分布乱数を付加
y_preds = (mu_preds.T + stats.norm.rvs(loc=0, scale=sigmas, random_state=123)).T$${\mathtt{scipy.stats}}$$ の $${\mathtt{norm.rvs()}}$$ で正規分布乱数（平均：０、標準偏差：$${\sigma}$$ のMCMCサンプル）を予測平均 $${\mu}$$ に付加して、観測値ベース（ばらつき込み）の予測値を算出します。
以上で デザイン行列を用いた単回帰モデルの学びを終わりにします。
面白かったですね。
今回の記事は以上です。
楽しかったですね！
シリーズの記事
次の記事
前の記事
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の教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
いいなと思ったら応援しよう！
応援ありがとうございます。これからもがんばって記事を作成します！