ベイズ線形回帰

テンプレート:ベイズ統計学 ベイズ線形回帰（ベイズせんけいかいき、テンプレート:Lang-en-short）は、線形回帰モデルの回帰係数や誤差分散を確率変数として扱い、事前分布と観測データに基づいてそれらの事後分布を求めるベイズ推定の手法である。回帰係数を一つの値として求めるだけでなく、その不確実性を事後分布として表し、将来の応答変数については事後予測分布を求めることができる。

正規誤差を仮定した線形回帰では、正規分布などの共役事前分布を用いると事後分布を解析的に求められる。より複雑な事前分布、階層構造、外れ値に頑健な誤差分布などを用いる場合には、マルコフ連鎖モンテカルロ法や変分ベイズ法などの数値計算が利用される。

モデル

観測数を $n$、回帰係数の数を $p$ とする。応答変数の観測値ベクトルを $𝒚\in {ℝ}^n$、説明変数からなる設計行列を $X\in {ℝ}^{n\times p}$、回帰係数ベクトルを $𝜷\in {ℝ}^p$ とすると、通常の正規線形回帰モデルは

$𝒚=X𝜷+𝜺,𝜺\sim N(0,{\sigma }^2{I}_n)$

と表される。したがって尤度は

$𝒚\mid 𝜷,{\sigma }^2,X\sim N(X𝜷,{\sigma }^2{I}_n)$

である。ここで ${\sigma }^2$ は誤差分散、$I_n$ は $n$ 次の単位行列である。切片を含める場合は、設計行列の一列をすべて1とする。通常は説明変数 $X$ を所与として、$𝜷$ と ${\sigma }^2$ の条件付き事後分布を求める。

事前分布と事後分布

誤差分散が既知の場合

誤差分散 ${\sigma }^2$ が既知で、回帰係数に正規事前分布

$𝜷\sim N({𝒎}_0,{\mathrm{\Sigma }}_0)$

を置くと、事後分布も多変量正規分布となる。事後共分散行列 ${\mathrm{\Sigma }}_n$ と事後平均 ${𝒎}_n$ は

${\mathrm{\Sigma }}_n={({\mathrm{\Sigma }}_0^{-1}+\frac{1}{{\sigma }^2}{X}^{𝖳}X)}^{-1},$

${𝒎}_n={\mathrm{\Sigma }}_n({\mathrm{\Sigma }}_0^{-1}{𝒎}_0+\frac{1}{{\sigma }^2}{X}^{𝖳}𝒚)$

であり、

$𝜷\mid 𝒚,X,{\sigma }^2\sim N({𝒎}_n,{\mathrm{\Sigma }}_n)$

となる。事後精度 ${\mathrm{\Sigma }}_n^{-1}$ は、事前精度 ${\mathrm{\Sigma }}_0^{-1}$ とデータが与える精度 $X^{𝖳}X/{\sigma }^2$ の和になっている。

誤差分散が未知の場合

誤差分散も未知の場合、代表的な共役事前分布は

$𝜷\mid {\sigma }^2\sim N({𝒎}_0,{\sigma }^2{V}_0),{\sigma }^2\sim \mathrm{IG}(a_0,b_0)$

である。ここでは逆ガンマ分布を

$p({\sigma }^2)\propto ({\sigma }^2{)}^{-(a_0+1)}\exp (-\frac{b_0}{{\sigma }^2})$

と母数化する。この同時分布は正規逆ガンマ分布と呼ばれる。精度 $\tau =1/{\sigma }^2$ を用いれば、同じ共役構造を正規ガンマ分布で表せる。

事後ハイパーパラメータを

$V_n=(V_0^{-1}+X^{𝖳}X{)}^{-1},$

${𝒎}_n=V_n(V_0^{-1}{𝒎}_0+X^{𝖳}𝒚),$

$a_n=a_0+\frac{n}{2},$

$b_n=b_0+\frac{1}{2}({𝒚}^{𝖳}𝒚+{𝒎}_0^{𝖳}{V}_0^{-1}{𝒎}_0-{𝒎}_n^{𝖳}{V}_n^{-1}{𝒎}_n)$

とすると、

$𝜷\mid {\sigma }^2,𝒚,X\sim N({𝒎}_n,{\sigma }^2{V}_n),{\sigma }^2\mid 𝒚,X\sim \mathrm{IG}(a_n,b_n)$

となる。${\sigma }^2$ を積分消去した回帰係数の周辺事後分布は、自由度 $2a_n$ のテンプレート:仮リンクである。

広く用いられる無情報事前分布の一つは

$p(𝜷,{\sigma }^2)\propto \frac{1}{{\sigma }^2}$

である。この事前分布は非正規化事前分布であるため、事後分布が正規化可能かを確認する必要がある。設計行列の階数が $p$ であり、$n>p$ ならば事後分布は正規化可能となる。

事後予測分布

新しい説明変数ベクトルを ${𝒙}_{\ast }$、対応する未知の応答変数を $y_{\ast }$ とする。誤差分散が既知の場合、回帰係数を事後分布について積分すると

$y_{\ast }\mid 𝒚,X,{𝒙}_{\ast },{\sigma }^2\sim N({𝒙}_{\ast }^{𝖳}{𝒎}_n,{\sigma }^2+{𝒙}_{\ast }^{𝖳}{\mathrm{\Sigma }}_n{𝒙}_{\ast })$

を得る。予測分散の第1項は新しい観測に固有のばらつき、第2項は回帰係数が確定していないことによる不確実性を表す。

誤差分散も未知で正規逆ガンマ事前分布を用いた場合、事後予測分布は自由度 $2a_n$、位置 ${𝒙}_{\ast }^{𝖳}{𝒎}_n$、尺度の二乗

$\frac{b_n}{a_n}(1+{𝒙}_{\ast }^{𝖳}{V}_n{𝒙}_{\ast })$

のt分布となる。このように、将来の観測値に関する不確実性を確率分布として表すことができる。

最小二乗法および正則化との関係

設計行列が列フルランクで、回帰係数に平坦な事前分布を置いた極限では、事後平均は最小二乗法による推定量

${\widehat{𝜷}}_{\mathrm{O}\mathrm{L}\mathrm{S}}=(X^{𝖳}X{)}^{-1}{X}^{𝖳}𝒚$

に一致する。また、標準的な無情報事前分布のもとでは、回帰係数の周辺信用区間が通常の線形回帰における信頼区間と数値的に一致する場合があるが、両者の確率解釈は異なる。

誤差分散が既知で、

$𝜷\sim N(0,{\tau }^2{I}_p)$

とすると、事後平均と最大事後確率推定値は

${𝒎}_n={(X^{𝖳}X+\lambda I_p)}^{-1}{X}^{𝖳}𝒚,\lambda =\frac{{\sigma }^2}{{\tau }^2}$

となり、リッジ回帰の推定量と一致する。このように正規事前分布は回帰係数を事前平均へ縮小する働きをもつ。事前共分散を $g{\sigma }^2(X^{𝖳}X{)}^{-1}$ とするテンプレート:仮リンクも、回帰モデルの変数選択やモデル比較で用いられる。

回帰係数や回帰式の集合に交換可能性を仮定し、それらを共通の分布から生じたものとして扱うと、複数の係数や群の間で情報を共有する階層ベイズモデルになる。LindleyとSmithは、この考え方を多段階の線形モデルとして定式化し、推定値が全体平均などへ縮小されることや、リッジ回帰に近い推定式が得られることを示した。

モデル評価と拡張

モデルの適合度は、残差の図示に加え、事後分布から反復データを生成して観測データと比較する事後予測チェックによって調べられる。例えば、残差と当てはめ値の関係、分散の不均一性、外れ値、誤差の相関などを検討する。

説明変数の組合せや異なる回帰モデルを比較する方法には、周辺尤度、ベイズ因子、予測性能に基づく基準などがある。各モデルの不確実性を残したまま予測を平均する方法はテンプレート:仮リンクと呼ばれる。

通常の正規線形回帰は、誤差が独立で同じ分散をもつことを仮定する。この仮定を緩めることで、不均一分散、相関のある誤差、テンプレート:仮リンク、階層線形モデルなどへ拡張できる。また、回帰係数にラプラス分布やスパイク・アンド・スラブ分布などの事前分布を置くことで、縮小推定や変数選択を行うこともできる。

計算とソフトウェア

共役事前分布を用いる基本モデルでは、行列分解と既知の確率分布からの乱数生成だけで事後分布を計算できる。数値的な安定性のため、$(X^{𝖳}X{)}^{-1}$ を直接計算する代わりに、QR分解やコレスキー分解を用いて連立方程式を解くことが多い。非共役事前分布や階層構造を含むモデルでは、汎用的なベイズ計算ソフトウェアが用いられる。

環境	主な機能
Python・scikit-learn	sklearn.linear_model.BayesianRidge は、回帰係数と誤差の精度に事前分布を置き、周辺尤度の最大化によって正則化のハイパーパラメータを推定するベイズリッジ回帰を実装している。任意の事前分布を指定する汎用ベイズ推論とは目的が異なる。
Python・テンプレート:仮リンク	Pythonコードで確率モデルを記述し、マルコフ連鎖モンテカルロ法や変分推論によって線形回帰、階層回帰、ロバスト回帰などを推定できる。
テンプレート:仮リンク	尤度と事前分布をモデルとして記述し、勾配を用いるマルコフ連鎖モンテカルロ法などで事後分布を計算する。単純な線形回帰から多階層の一般化線形モデルまで記述できる。PythonからはCmdStanPyなどのインターフェースを利用できる。
R	brms と rstanarm はStanを計算基盤として用い、Rの数式形式に近い記法でベイズ回帰モデルを指定する。
MATLAB	Econometrics Toolboxの bayeslm は、回帰係数と誤差分散の事前モデルを作成し、事後分布の推定や予測を行う機能を提供する。

歴史

RaiffaとSchlaiferは1961年の『Applied Statistical Decision Theory』で、正規回帰過程に対する自然共役事前分布と事後更新を体系的に扱った。1972年、LindleyとSmithは通常の線形モデルをベイズ的に再構成し、回帰係数に交換可能性と多段階の事前分布を導入した。彼らの枠組みは、複数の回帰式や実験計画における情報共有、縮小推定、階層線形モデルの発展に影響を与えた。

関連項目

• 線形回帰
• ベイズ推定
• 共役事前分布
• 事後予測分布
• リッジ回帰
• 階層ベイズモデル
• ベイズ因子
• ガウス過程