ベイズ線形回帰 (ベイズせんけいかいき、英 : Bayesian linear regression )は、線形回帰 モデルの回帰係数や誤差分散を確率変数として扱い、事前分布と観測データに基づいてそれらの事後分布を求めるベイズ推定 の手法である[ 1] [ 2] 。回帰係数を一つの値として求めるだけでなく、その不確実性を事後分布として表し、将来の応答変数については事後予測分布 を求めることができる。
正規誤差を仮定した線形回帰では、正規分布などの共役事前分布 を用いると事後分布を解析的に求められる。より複雑な事前分布、階層構造、外れ値に頑健な誤差分布などを用いる場合には、マルコフ連鎖モンテカルロ法 や変分ベイズ法 などの数値計算が利用される[ 3] [ 2] 。
モデル
観測数を
n
{\displaystyle n}
、回帰係数の数を
p
{\displaystyle p}
とする。応答変数の観測値ベクトルを
y
∈
R
n
{\displaystyle {\boldsymbol {y}}\in \mathbb {R} ^{n}}
、説明変数からなる設計行列を
X
∈
R
n
×
p
{\displaystyle X\in \mathbb {R} ^{n\times p}}
、回帰係数ベクトルを
β
∈
R
p
{\displaystyle {\boldsymbol {\beta }}\in \mathbb {R} ^{p}}
とすると、通常の正規線形回帰モデルは
y
=
X
β
+
ε
,
ε
∼
N
(
0
,
σ
2
I
n
)
{\displaystyle {\boldsymbol {y}}=X{\boldsymbol {\beta }}+{\boldsymbol {\varepsilon }},\qquad {\boldsymbol {\varepsilon }}\sim N({\boldsymbol {0}},\sigma ^{2}I_{n})}
と表される。したがって尤度は
y
∣
β
,
σ
2
,
X
∼
N
(
X
β
,
σ
2
I
n
)
{\displaystyle {\boldsymbol {y}}\mid {\boldsymbol {\beta }},\sigma ^{2},X\sim N(X{\boldsymbol {\beta }},\sigma ^{2}I_{n})}
である。ここで
σ
2
{\displaystyle \sigma ^{2}}
は誤差分散、
I
n
{\displaystyle I_{n}}
は
n
{\displaystyle n}
次の単位行列である。切片を含める場合は、設計行列の一列をすべて1とする。通常は説明変数
X
{\displaystyle X}
を所与として、
β
{\displaystyle {\boldsymbol {\beta }}}
と
σ
2
{\displaystyle \sigma ^{2}}
の条件付き事後分布を求める[ 2] [ 4] 。
事前分布と事後分布
誤差分散が既知の場合
誤差分散
σ
2
{\displaystyle \sigma ^{2}}
が既知で、回帰係数に正規事前分布
β
∼
N
(
m
0
,
Σ
0
)
{\displaystyle {\boldsymbol {\beta }}\sim N({\boldsymbol {m}}_{0},\Sigma _{0})}
を置くと、事後分布も多変量正規分布となる。事後共分散行列
Σ
n
{\displaystyle \Sigma _{n}}
と事後平均
m
n
{\displaystyle {\boldsymbol {m}}_{n}}
は
Σ
n
=
(
Σ
0
−
1
+
1
σ
2
X
T
X
)
−
1
,
{\displaystyle \Sigma _{n}=\left(\Sigma _{0}^{-1}+{\frac {1}{\sigma ^{2}}}X^{\mathsf {T}}X\right)^{-1},}
m
n
=
Σ
n
(
Σ
0
−
1
m
0
+
1
σ
2
X
T
y
)
{\displaystyle {\boldsymbol {m}}_{n}=\Sigma _{n}\left(\Sigma _{0}^{-1}{\boldsymbol {m}}_{0}+{\frac {1}{\sigma ^{2}}}X^{\mathsf {T}}{\boldsymbol {y}}\right)}
であり、
β
∣
y
,
X
,
σ
2
∼
N
(
m
n
,
Σ
n
)
{\displaystyle {\boldsymbol {\beta }}\mid {\boldsymbol {y}},X,\sigma ^{2}\sim N({\boldsymbol {m}}_{n},\Sigma _{n})}
となる[ 3] [ 4] 。事後精度
Σ
n
−
1
{\displaystyle \Sigma _{n}^{-1}}
は、事前精度
Σ
0
−
1
{\displaystyle \Sigma _{0}^{-1}}
とデータが与える精度
X
T
X
/
σ
2
{\displaystyle X^{\mathsf {T}}X/\sigma ^{2}}
の和になっている。
誤差分散が未知の場合
誤差分散も未知の場合、代表的な共役事前分布は
β
∣
σ
2
∼
N
(
m
0
,
σ
2
V
0
)
,
σ
2
∼
IG
(
a
0
,
b
0
)
{\displaystyle {\boldsymbol {\beta }}\mid \sigma ^{2}\sim N({\boldsymbol {m}}_{0},\sigma ^{2}V_{0}),\qquad \sigma ^{2}\sim \operatorname {IG} (a_{0},b_{0})}
である。ここでは逆ガンマ分布を
p
(
σ
2
)
∝
(
σ
2
)
−
(
a
0
+
1
)
exp
(
−
b
0
σ
2
)
{\displaystyle p(\sigma ^{2})\propto (\sigma ^{2})^{-(a_{0}+1)}\exp \left(-{\frac {b_{0}}{\sigma ^{2}}}\right)}
と母数化する。この同時分布は正規逆ガンマ分布 と呼ばれる。精度
τ
=
1
/
σ
2
{\displaystyle \tau =1/\sigma ^{2}}
を用いれば、同じ共役構造を正規ガンマ分布 で表せる[ 3] 。
事後ハイパーパラメータを
V
n
=
(
V
0
−
1
+
X
T
X
)
−
1
,
{\displaystyle V_{n}=(V_{0}^{-1}+X^{\mathsf {T}}X)^{-1},}
m
n
=
V
n
(
V
0
−
1
m
0
+
X
T
y
)
,
{\displaystyle {\boldsymbol {m}}_{n}=V_{n}(V_{0}^{-1}{\boldsymbol {m}}_{0}+X^{\mathsf {T}}{\boldsymbol {y}}),}
a
n
=
a
0
+
n
2
,
{\displaystyle a_{n}=a_{0}+{\frac {n}{2}},}
b
n
=
b
0
+
1
2
(
y
T
y
+
m
0
T
V
0
−
1
m
0
−
m
n
T
V
n
−
1
m
n
)
{\displaystyle b_{n}=b_{0}+{\frac {1}{2}}\left({\boldsymbol {y}}^{\mathsf {T}}{\boldsymbol {y}}+{\boldsymbol {m}}_{0}^{\mathsf {T}}V_{0}^{-1}{\boldsymbol {m}}_{0}-{\boldsymbol {m}}_{n}^{\mathsf {T}}V_{n}^{-1}{\boldsymbol {m}}_{n}\right)}
とすると、
β
∣
σ
2
,
y
,
X
∼
N
(
m
n
,
σ
2
V
n
)
,
σ
2
∣
y
,
X
∼
IG
(
a
n
,
b
n
)
{\displaystyle {\boldsymbol {\beta }}\mid \sigma ^{2},{\boldsymbol {y}},X\sim N({\boldsymbol {m}}_{n},\sigma ^{2}V_{n}),\qquad \sigma ^{2}\mid {\boldsymbol {y}},X\sim \operatorname {IG} (a_{n},b_{n})}
となる。
σ
2
{\displaystyle \sigma ^{2}}
を積分消去した回帰係数の周辺事後分布は、自由度
2
a
n
{\displaystyle 2a_{n}}
の多変量t分布 (英語版 ) である[ 3] 。
広く用いられる無情報事前分布の一つは
p
(
β
,
σ
2
)
∝
1
σ
2
{\displaystyle p({\boldsymbol {\beta }},\sigma ^{2})\propto {\frac {1}{\sigma ^{2}}}}
である。この事前分布は非正規化事前分布であるため、事後分布が正規化可能かを確認する必要がある。設計行列の階数が
p
{\displaystyle p}
であり、
n
>
p
{\displaystyle n>p}
ならば事後分布は正規化可能となる[ 2] 。
事後予測分布
新しい説明変数ベクトルを
x
∗
{\displaystyle {\boldsymbol {x}}_{*}}
、対応する未知の応答変数を
y
∗
{\displaystyle y_{*}}
とする。誤差分散が既知の場合、回帰係数を事後分布について積分すると
y
∗
∣
y
,
X
,
x
∗
,
σ
2
∼
N
(
x
∗
T
m
n
,
σ
2
+
x
∗
T
Σ
n
x
∗
)
{\displaystyle y_{*}\mid {\boldsymbol {y}},X,{\boldsymbol {x}}_{*},\sigma ^{2}\sim N\left({\boldsymbol {x}}_{*}^{\mathsf {T}}{\boldsymbol {m}}_{n},\,\sigma ^{2}+{\boldsymbol {x}}_{*}^{\mathsf {T}}\Sigma _{n}{\boldsymbol {x}}_{*}\right)}
を得る[ 3] 。予測分散の第1項は新しい観測に固有のばらつき、第2項は回帰係数が確定していないことによる不確実性を表す。
誤差分散も未知で正規逆ガンマ事前分布を用いた場合、事後予測分布は自由度
2
a
n
{\displaystyle 2a_{n}}
、位置
x
∗
T
m
n
{\displaystyle {\boldsymbol {x}}_{*}^{\mathsf {T}}{\boldsymbol {m}}_{n}}
、尺度の二乗
b
n
a
n
(
1
+
x
∗
T
V
n
x
∗
)
{\displaystyle {\frac {b_{n}}{a_{n}}}\left(1+{\boldsymbol {x}}_{*}^{\mathsf {T}}V_{n}{\boldsymbol {x}}_{*}\right)}
のt分布 となる[ 3] [ 2] 。このように、将来の観測値に関する不確実性を確率分布として表すことができる。
最小二乗法および正則化との関係
設計行列が列フルランクで、回帰係数に平坦な事前分布を置いた極限では、事後平均は最小二乗法 による推定量
β
^
O
L
S
=
(
X
T
X
)
−
1
X
T
y
{\displaystyle {\hat {\boldsymbol {\beta }}}_{\mathrm {OLS} }=(X^{\mathsf {T}}X)^{-1}X^{\mathsf {T}}{\boldsymbol {y}}}
に一致する。また、標準的な無情報事前分布のもとでは、回帰係数の周辺信用区間が通常の線形回帰における信頼区間と数値的に一致する場合があるが、両者の確率解釈は異なる[ 2] [ 3] 。
誤差分散が既知で、
β
∼
N
(
0
,
τ
2
I
p
)
{\displaystyle {\boldsymbol {\beta }}\sim N({\boldsymbol {0}},\tau ^{2}I_{p})}
とすると、事後平均と最大事後確率推定値は
m
n
=
(
X
T
X
+
λ
I
p
)
−
1
X
T
y
,
λ
=
σ
2
τ
2
{\displaystyle {\boldsymbol {m}}_{n}=\left(X^{\mathsf {T}}X+\lambda I_{p}\right)^{-1}X^{\mathsf {T}}{\boldsymbol {y}},\qquad \lambda ={\frac {\sigma ^{2}}{\tau ^{2}}}}
となり、リッジ回帰 の推定量と一致する[ 3] 。このように正規事前分布は回帰係数を事前平均へ縮小する働きをもつ。事前共分散を
g
σ
2
(
X
T
X
)
−
1
{\displaystyle g\sigma ^{2}(X^{\mathsf {T}}X)^{-1}}
とするg-事前分布 (英語版 ) も、回帰モデルの変数選択やモデル比較で用いられる。
回帰係数や回帰式の集合に交換可能性を仮定し、それらを共通の分布から生じたものとして扱うと、複数の係数や群の間で情報を共有する階層ベイズモデル になる。LindleyとSmithは、この考え方を多段階の線形モデルとして定式化し、推定値が全体平均などへ縮小されることや、リッジ回帰に近い推定式が得られることを示した[ 5] 。
モデル評価と拡張
モデルの適合度は、残差の図示に加え、事後分布から反復データを生成して観測データと比較する事後予測チェックによって調べられる。例えば、残差と当てはめ値の関係、分散の不均一性、外れ値、誤差の相関などを検討する[ 2] 。
説明変数の組合せや異なる回帰モデルを比較する方法には、周辺尤度 、ベイズ因子 、予測性能に基づく基準などがある。各モデルの不確実性を残したまま予測を平均する方法はベイズモデル平均 (英語版 ) と呼ばれる[ 1] [ 3] 。
通常の正規線形回帰は、誤差が独立で同じ分散をもつことを仮定する。この仮定を緩めることで、不均一分散、相関のある誤差、ロバスト回帰 (英語版 ) 、階層線形モデルなどへ拡張できる。また、回帰係数にラプラス分布やスパイク・アンド・スラブ分布などの事前分布を置くことで、縮小推定や変数選択を行うこともできる[ 2] 。
計算とソフトウェア
共役事前分布を用いる基本モデルでは、行列分解と既知の確率分布からの乱数生成だけで事後分布を計算できる。数値的な安定性のため、
(
X
T
X
)
−
1
{\displaystyle (X^{\mathsf {T}}X)^{-1}}
を直接計算する代わりに、QR分解 やコレスキー分解 を用いて連立方程式を解くことが多い[ 2] 。非共役事前分布や階層構造を含むモデルでは、汎用的なベイズ計算ソフトウェアが用いられる。
環境
主な機能
Python ・Sscikit-learn
sklearn.linear_model.BayesianRidge は、回帰係数と誤差の精度に事前分布を置き、周辺尤度の最大化によって正則化のハイパーパラメータを推定するベイズリッジ回帰を実装している[ 6] 。任意の事前分布を指定する汎用ベイズ推論とは目的が異なる。
Python・PyMC (英語版 )
Pythonコードで確率モデルを記述し、マルコフ連鎖モンテカルロ法や変分推論によって線形回帰、階層回帰、ロバスト回帰などを推定できる[ 7] 。
Stan (英語版 )
尤度と事前分布をモデルとして記述し、勾配を用いるマルコフ連鎖モンテカルロ法などで事後分布を計算する。単純な線形回帰から多階層の一般化線形モデルまで記述できる[ 8] 。PythonからはCmdStanPyなどのインターフェースを利用できる[ 9] 。
R
brms と rstanarm はStanを計算基盤として用い、Rの数式形式に近い記法でベイズ回帰モデルを指定する[ 10] [ 11] 。
MATLAB
Econometrics Toolboxの bayeslm は、回帰係数と誤差分散の事前モデルを作成し、事後分布の推定や予測を行う機能を提供する[ 12] 。
歴史
RaiffaとSchlaiferは1961年の『Applied Statistical Decision Theory』で、正規回帰過程に対する自然共役事前分布と事後更新を体系的に扱った[ 4] 。1972年、LindleyとSmithは通常の線形モデルをベイズ的に再構成し、回帰係数に交換可能性と多段階の事前分布を導入した。彼らの枠組みは、複数の回帰式や実験計画における情報共有、縮小推定、階層線形モデルの発展に影響を与えた[ 5] 。
脚注
1 2 ピーター・D・ホフ『標準ベイズ統計学』入江薫・菅澤翔之助・橋本真太郎 訳、朝倉書店、2022年、166–190頁。ISBN 978-4-254-12267-1 。
1 2 3 4 5 6 7 8 9 Gelman, Andrew; Carlin, John B.; Stern, Hal S.; Dunson, David B.; Vehtari, Aki; Rubin, Donald B. (2013) (英語). Bayesian Data Analysis (3rd ed.). CRC Press. pp. 353–379. ISBN 978-1-4398-4095-5
1 2 3 4 5 6 7 8 9 Murphy, Kevin P. (2012) (英語). Machine Learning: A Probabilistic Perspective . MIT Press. pp. 231–239. ISBN 978-0-262-01802-9
1 2 3 Raiffa, Howard; Schlaifer, Robert (1961) (英語). Applied Statistical Decision Theory . Division of Research, Graduate School of Business Administration, Harvard University. pp. 334–348
1 2 Lindley, D. V.; Smith, A. F. M. (1972). “Bayes Estimates for the Linear Model” (英語). Journal of the Royal Statistical Society: Series B (Methodological) 34 (1): 1–18. doi :10.1111/j.2517-6161.1972.tb00885.x .
↑ “BayesianRidge ” (英語). scikit-learn documentation . 2026年7月16日閲覧。
↑ “GLM: Linear regression ” (英語). PyMC documentation . 2026年7月16日閲覧。
↑ “Regression Models ” (英語). Stan User's Guide . 2026年7月16日閲覧。
↑ “cmdstanpy – Python interface to CmdStan ” (英語). CmdStanPy documentation . 2026年7月16日閲覧。
↑ “brms: Bayesian Regression Models using 'Stan' ” (英語). The Comprehensive R Archive Network . 2026年7月16日閲覧。
↑ “Bayesian Applied Regression Modeling via Stan ” (英語). rstanarm . 2026年7月16日閲覧。
↑ “bayeslm - Create Bayesian linear regression model object ” (英語). MathWorks . 2026年7月16日閲覧。