ベイズ多変量線形回帰

統計学において、ベイズ多変量線形回帰は、多変量線形回帰、つまり予測結果が単一のスカラーランダム変数ではなく相関ランダム変数のベクトルである線形回帰に対するベイズ的アプローチです 。このアプローチのより一般的な扱いについては、記事MMSE 推定量で説明されています。

予測される従属変数が単一の実数値スカラーではなく、相関する実数の長さmベクトルである回帰問題を考えてみましょう。標準の回帰設定と同様に、 n 個の観測値があり、各観測値iはk −1 個の説明変数で構成され、長さkのベクトルにグループ化されています(切片係数を可能にするために、値が 1 のダミー変数が追加されています)。これは、各観測値iについて、 m 個の関連する回帰問題のセットとして考えることができます。 ここで、誤差のセットはすべて相関しています。同様に、これは、結果が行ベクトルで、回帰係数ベクトルが次のように隣り合って積み重ねられている 単一の回帰問題として考えることもできます。${\displaystyle \mathbf {x} _{i}}$$${\displaystyle {\begin{aligned}y_{i,1}&=\mathbf {x} _{i}^{\mathsf {T}}{\boldsymbol {\beta }}_{1}+\epsilon _{i,1}\\&\;\;\vdots \\y_{i,m}&=\mathbf {x} _{i}^{\mathsf {T}}{\boldsymbol {\beta }}_{m}+\epsilon _{i,m}\end{aligned}}}$$${\displaystyle \{\epsilon _{i,1},\ldots ,\epsilon _{i,m}\}}$ ${\displaystyle \mathbf {y} _{i}^{\mathsf {T}}}$$${\displaystyle \mathbf {y} _{i}^{\mathsf {T}}=\mathbf {x} _{i}^{\mathsf {T}}\mathbf {B} +{\boldsymbol {\epsilon }}_{i}^{\mathsf {T}}.}$$

係数行列Bは、各回帰問題の 係数ベクトルが水平に積み重ねられた行列です。${\displaystyle k\times m}$${\displaystyle {\boldsymbol {\beta }}_{1},\ldots ,{\boldsymbol {\beta }}_{m}}$$${\displaystyle \mathbf {B} ={\begin{bmatrix}{\begin{pmatrix}\\{\boldsymbol {\beta }}_{1}\\\\\end{pmatrix}}\cdots {\begin{pmatrix}\\{\boldsymbol {\beta }}_{m}\\\\\end{pmatrix}}\end{bmatrix}}={\begin{bmatrix}{\begin{pmatrix}\beta _{1,1}\\\vdots \\\beta _{k,1}\end{pmatrix}}\cdots {\begin{pmatrix}\beta _{1,m}\\\vdots \\\beta _{k,m}\end{pmatrix}}\end{bmatrix}}.}$$

各観測値iのノイズ ベクトルは正規分布に従うため、特定の観測値の結果は相関します。 ${\displaystyle {\boldsymbol {\epsilon }}_{i}}$$${\displaystyle {\boldsymbol {\epsilon }}_{i}\sim N(0,{\boldsymbol {\Sigma }}_{\epsilon }).}$$

回帰問題全体を行列形式で記述すると、次のようになります 。 ここで、YとEは行列です。設計行列X は、標準的な線形回帰の設定と同様に、観測値が垂直に積み重ねられた行列です。 $${\displaystyle \mathbf {Y} =\mathbf {X} \mathbf {B} +\mathbf {E} ,}$$${\displaystyle n\times m}$ ${\displaystyle n\times k}$$${\displaystyle \mathbf {X} ={\begin{bmatrix}\mathbf {x} _{1}^{\mathsf {T}}\\\mathbf {x} _{2}^{\mathsf {T}}\\\vdots \\\mathbf {x} _{n}^{\mathsf {T}}\end{bmatrix}}={\begin{bmatrix}x_{1,1}&\cdots &x_{1,k}\\x_{2,1}&\cdots &x_{2,k}\\\vdots &\ddots &\vdots \\x_{n,1}&\cdots &x_{n,k}\end{bmatrix}}.}$$

古典的な頻度主義の線形最小二乗法は、ムーア・ペンローズ擬似逆行列を使用して 回帰係数の行列を単純に推定することです。 ${\displaystyle {\hat {\mathbf {B} }}}$ $${\displaystyle {\hat {\mathbf {B} }}=(\mathbf {X} ^{\mathsf {T}}\mathbf {X} )^{-1}\mathbf {X} ^{\mathsf {T}}\mathbf {Y} .}$$

ベイズ解を得るには、条件付き尤度を指定してから適切な共役事前分布を見つける必要があります。線形ベイズ回帰の単変量の場合と同様に、自然な条件付き共役事前分布（スケール依存）を指定できることがわかります。

条件付き尤度を次のように書き表すと^[1]、 誤差を と で書き表す と、$${\displaystyle \rho (\mathbf {E} |{\boldsymbol {\Sigma }}_{\epsilon })\propto |{\boldsymbol {\Sigma }}_{\epsilon }|^{-n/2}\exp \left(-{\tfrac {1}{2}}\operatorname {tr} \left(\mathbf {E} ^{\mathsf {T}}\mathbf {E} {\boldsymbol {\Sigma }}_{\epsilon }^{-1}\right)\right),}$$${\displaystyle \mathbf {E} }$${\displaystyle \mathbf {Y} ,\mathbf {X} ,}$${\displaystyle \mathbf {B} }$$${\displaystyle \rho (\mathbf {Y} |\mathbf {X} ,\mathbf {B} ,{\boldsymbol {\Sigma }}_{\epsilon })\propto |{\boldsymbol {\Sigma }}_{\epsilon }|^{-n/2}\exp(-{\tfrac {1}{2}}\operatorname {tr} ((\mathbf {Y} -\mathbf {X} \mathbf {B} )^{\mathsf {T}}(\mathbf {Y} -\mathbf {X} \mathbf {B} ){\boldsymbol {\Sigma }}_{\epsilon }^{-1})),}$$

私たちは自然な共役事前分布、 つまり尤度と同じ関数形式の結合密度を求めます。尤度は の 2 次関数なので、尤度を書き直して（古典的なサンプル推定値からの偏差） で正規分布になるようにします。 ${\displaystyle \rho (\mathbf {B} ,\Sigma _{\epsilon })}$${\displaystyle \mathbf {B} }$${\displaystyle (\mathbf {B} -{\hat {\mathbf {B} }})}$

ベイズ線形回帰と同じ手法を使用して、指数項を行列形式の二乗和手法で分解します。ただし、ここでは行列微分積分（クロネッカー積とベクトル化変換）も使用する必要があります。

まず、二乗和を適用して尤度の新しい表現を取得します。 $${\displaystyle \rho (\mathbf {Y} |\mathbf {X} ,\mathbf {B} ,{\boldsymbol {\Sigma }}_{\epsilon })\propto |{\boldsymbol {\Sigma }}_{\epsilon }|^{-(n-k)/2}\exp(-\operatorname {tr} ({\tfrac {1}{2}}\mathbf {S} ^{\mathsf {T}}\mathbf {S} {\boldsymbol {\Sigma }}_{\epsilon }^{-1}))|{\boldsymbol {\Sigma }}_{\epsilon }|^{-k/2}\exp(-{\tfrac {1}{2}}\operatorname {tr} ((\mathbf {B} -{\hat {\mathbf {B} }})^{\mathsf {T}}\mathbf {X} ^{\mathsf {T}}\mathbf {X} (\mathbf {B} -{\hat {\mathbf {B} }}){\boldsymbol {\Sigma }}_{\epsilon }^{-1})),}$$ $${\displaystyle \mathbf {S} =\mathbf {Y} -\mathbf {X} {\hat {\mathbf {B} }}}$$

事前分布の条件付き形式を開発します。 ここで、は逆ウィッシュアート分布で 、 は行列 内の何らかの形式の正規分布です。これは、尤度を行列の関数からベクトルの関数に変換するベクトル化変換を使用して実現されます。 $${\displaystyle \rho (\mathbf {B} ,{\boldsymbol {\Sigma }}_{\epsilon })=\rho ({\boldsymbol {\Sigma }}_{\epsilon })\rho (\mathbf {B} |{\boldsymbol {\Sigma }}_{\epsilon }),}$$${\displaystyle \rho ({\boldsymbol {\Sigma }}_{\epsilon })}$${\displaystyle \rho (\mathbf {B} |{\boldsymbol {\Sigma }}_{\epsilon })}$${\displaystyle \mathbf {B} }$${\displaystyle \mathbf {B} ,{\hat {\mathbf {B} }}}$${\displaystyle {\boldsymbol {\beta }}=\operatorname {vec} (\mathbf {B} ),{\hat {\boldsymbol {\beta }}}=\operatorname {vec} ({\hat {\mathbf {B} }})}$

書く $${\displaystyle \operatorname {tr} ((\mathbf {B} -{\hat {\mathbf {B} }})^{\mathsf {T}}\mathbf {X} ^{\mathsf {T}}\mathbf {X} (\mathbf {B} -{\hat {\mathbf {B} }}){\boldsymbol {\Sigma }}_{\epsilon }^{-1})=\operatorname {vec} (\mathbf {B} -{\hat {\mathbf {B} }})^{\mathsf {T}}\operatorname {vec} (\mathbf {X} ^{\mathsf {T}}\mathbf {X} (\mathbf {B} -{\hat {\mathbf {B} }}){\boldsymbol {\Sigma }}_{\epsilon }^{-1})}$$

ここで、は 行列AとBのクロネッカー積を表します。クロネッカー積は、行列を行列で乗算して2 つの行列の要素の積のあらゆる組み合わせからなる行列を 生成する外積の一般化です。$${\displaystyle \operatorname {vec} (\mathbf {X} ^{\mathsf {T}}\mathbf {X} (\mathbf {B} -{\hat {\mathbf {B} }}){\boldsymbol {\Sigma }}_{\epsilon }^{-1})=({\boldsymbol {\Sigma }}_{\epsilon }^{-1}\otimes \mathbf {X} ^{\mathsf {T}}\mathbf {X} )\operatorname {vec} (\mathbf {B} -{\hat {\mathbf {B} }}),}$$${\displaystyle \mathbf {A} \otimes \mathbf {B} }$${\displaystyle m\times n}$${\displaystyle p\times q}$${\displaystyle mp\times nq}$

すると 、 では通常の確率になります。 $${\displaystyle {\begin{aligned}&\operatorname {vec} (\mathbf {B} -{\hat {\mathbf {B} }})^{\mathsf {T}}({\boldsymbol {\Sigma }}_{\epsilon }^{-1}\otimes \mathbf {X} ^{\mathsf {T}}\mathbf {X} )\operatorname {vec} (\mathbf {B} -{\hat {\mathbf {B} }})\\&=({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})^{\mathsf {T}}({\boldsymbol {\Sigma }}_{\epsilon }^{-1}\otimes \mathbf {X} ^{\mathsf {T}}\mathbf {X} )({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})\end{aligned}}}$$${\displaystyle ({\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}})}$

尤度をより扱いやすい形式で使用することで、自然な（条件付き）共役事前分布を見つけることができます。

ベクトル化された変数を用いた自然共役事前分布は、次の形式となる。^[1] ここで 、 ${\displaystyle {\boldsymbol {\beta }}}$ $${\displaystyle \rho ({\boldsymbol {\beta }},{\boldsymbol {\Sigma }}_{\epsilon })=\rho ({\boldsymbol {\Sigma }}_{\epsilon })\rho ({\boldsymbol {\beta }}|{\boldsymbol {\Sigma }}_{\epsilon }),}$$$${\displaystyle \rho ({\boldsymbol {\Sigma }}_{\epsilon })\sim {\mathcal {W}}^{-1}(\mathbf {V} _{0},{\boldsymbol {\nu }}_{0})}$$$${\displaystyle \rho ({\boldsymbol {\beta }}|{\boldsymbol {\Sigma }}_{\epsilon })\sim N({\boldsymbol {\beta }}_{0},{\boldsymbol {\Sigma }}_{\epsilon }\otimes {\boldsymbol {\Lambda }}_{0}^{-1}).}$$

上記の事前分布と尤度を用いると、事後分布は次のように表される: ^[1] ここで。 を含む項は次のようにグループ化できる ( ) : および $${\displaystyle {\begin{aligned}\rho ({\boldsymbol {\beta }},{\boldsymbol {\Sigma }}_{\epsilon }|\mathbf {Y} ,\mathbf {X} )\propto {}&|{\boldsymbol {\Sigma }}_{\epsilon }|^{-({\boldsymbol {\nu }}_{0}+m+1)/2}\exp {(-{\tfrac {1}{2}}\operatorname {tr} (\mathbf {V} _{0}{\boldsymbol {\Sigma }}_{\epsilon }^{-1}))}\\&\times |{\boldsymbol {\Sigma }}_{\epsilon }|^{-k/2}\exp {(-{\tfrac {1}{2}}\operatorname {tr} ((\mathbf {B} -\mathbf {B} _{0})^{\mathsf {T}}{\boldsymbol {\Lambda }}_{0}(\mathbf {B} -\mathbf {B} _{0}){\boldsymbol {\Sigma }}_{\epsilon }^{-1}))}\\&\times |{\boldsymbol {\Sigma }}_{\epsilon }|^{-n/2}\exp {(-{\tfrac {1}{2}}\operatorname {tr} ((\mathbf {Y} -\mathbf {XB} )^{\mathsf {T}}(\mathbf {Y} -\mathbf {XB} ){\boldsymbol {\Sigma }}_{\epsilon }^{-1}))},\end{aligned}}}$$${\displaystyle \operatorname {vec} (\mathbf {B} _{0})={\boldsymbol {\beta }}_{0}}$${\displaystyle \mathbf {B} }$${\displaystyle {\boldsymbol {\Lambda }}_{0}=\mathbf {U} ^{\mathsf {T}}\mathbf {U} }$$${\displaystyle {\begin{aligned}&\left(\mathbf {B} -\mathbf {B} _{0}\right)^{\mathsf {T}}{\boldsymbol {\Lambda }}_{0}\left(\mathbf {B} -\mathbf {B} _{0}\right)+\left(\mathbf {Y} -\mathbf {XB} \right)^{\mathsf {T}}\left(\mathbf {Y} -\mathbf {XB} \right)\\={}&\left({\begin{bmatrix}\mathbf {Y} \\\mathbf {U} \mathbf {B} _{0}\end{bmatrix}}-{\begin{bmatrix}\mathbf {X} \\\mathbf {U} \end{bmatrix}}\mathbf {B} \right)^{\mathsf {T}}\left({\begin{bmatrix}\mathbf {Y} \\\mathbf {U} \mathbf {B} _{0}\end{bmatrix}}-{\begin{bmatrix}\mathbf {X} \\\mathbf {U} \end{bmatrix}}\mathbf {B} \right)\\={}&\left({\begin{bmatrix}\mathbf {Y} \\\mathbf {U} \mathbf {B} _{0}\end{bmatrix}}-{\begin{bmatrix}\mathbf {X} \\\mathbf {U} \end{bmatrix}}\mathbf {B} _{n}\right)^{\mathsf {T}}\left({\begin{bmatrix}\mathbf {Y} \\\mathbf {U} \mathbf {B} _{0}\end{bmatrix}}-{\begin{bmatrix}\mathbf {X} \\\mathbf {U} \end{bmatrix}}\mathbf {B} _{n}\right)+\left(\mathbf {B} -\mathbf {B} _{n}\right)^{\mathsf {T}}\left(\mathbf {X} ^{\mathsf {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}\right)\left(\mathbf {B} -\mathbf {B} _{n}\right)\\={}&\left(\mathbf {Y} -\mathbf {X} \mathbf {B} _{n}\right)^{\mathsf {T}}\left(\mathbf {Y} -\mathbf {X} \mathbf {B} _{n}\right)+\left(\mathbf {B} _{0}-\mathbf {B} _{n}\right)^{\mathsf {T}}{\boldsymbol {\Lambda }}_{0}\left(\mathbf {B} _{0}-\mathbf {B} _{n}\right)+\left(\mathbf {B} -\mathbf {B} _{n}\right)^{\mathsf {T}}\left(\mathbf {X} ^{\mathsf {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}\right)\left(\mathbf {B} -\mathbf {B} _{n}\right),\end{aligned}}}$$$${\displaystyle \mathbf {B} _{n}=\left(\mathbf {X} ^{\mathsf {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}\right)^{-1}\left(\mathbf {X} ^{\mathsf {T}}\mathbf {X} {\hat {\mathbf {B} }}+{\boldsymbol {\Lambda }}_{0}\mathbf {B} _{0}\right)=\left(\mathbf {X} ^{\mathsf {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}\right)^{-1}\left(\mathbf {X} ^{\mathsf {T}}\mathbf {Y} +{\boldsymbol {\Lambda }}_{0}\mathbf {B} _{0}\right).}$$

これにより、事後分布をより便利な形式で記述できるようになりました。 $${\displaystyle {\begin{aligned}\rho ({\boldsymbol {\beta }},{\boldsymbol {\Sigma }}_{\epsilon }|\mathbf {Y} ,\mathbf {X} )\propto {}&|{\boldsymbol {\Sigma }}_{\epsilon }|^{-({\boldsymbol {\nu }}_{0}+m+n+1)/2}\exp {(-{\tfrac {1}{2}}\operatorname {tr} ((\mathbf {V} _{0}+(\mathbf {Y} -\mathbf {XB_{n}} )^{\mathsf {T}}(\mathbf {Y} -\mathbf {XB_{n}} )+(\mathbf {B} _{n}-\mathbf {B} _{0})^{\mathsf {T}}{\boldsymbol {\Lambda }}_{0}(\mathbf {B} _{n}-\mathbf {B} _{0})){\boldsymbol {\Sigma }}_{\epsilon }^{-1}))}\\&\times |{\boldsymbol {\Sigma }}_{\epsilon }|^{-k/2}\exp {(-{\tfrac {1}{2}}\operatorname {tr} ((\mathbf {B} -\mathbf {B} _{n})^{\mathsf {T}}(\mathbf {X} ^{T}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})(\mathbf {B} -\mathbf {B} _{n}){\boldsymbol {\Sigma }}_{\epsilon }^{-1}))}.\end{aligned}}}$$

これは逆ウィシャート分布と行列正規分布を掛け合わせた 形 をとります。 $${\displaystyle \rho ({\boldsymbol {\Sigma }}_{\epsilon }|\mathbf {Y} ,\mathbf {X} )\sim {\mathcal {W}}^{-1}(\mathbf {V} _{n},{\boldsymbol {\nu }}_{n})}$$$${\displaystyle \rho (\mathbf {B} |\mathbf {Y} ,\mathbf {X} ,{\boldsymbol {\Sigma }}_{\epsilon })\sim {\mathcal {MN}}_{k,m}(\mathbf {B} _{n},{\boldsymbol {\Lambda }}_{n}^{-1},{\boldsymbol {\Sigma }}_{\epsilon }).}$$

この事後分布のパラメータは次のように与えられます。 $${\displaystyle \mathbf {V} _{n}=\mathbf {V} _{0}+(\mathbf {Y} -\mathbf {XB_{n}} )^{\mathsf {T}}(\mathbf {Y} -\mathbf {XB_{n}} )+(\mathbf {B} _{n}-\mathbf {B} _{0})^{\mathsf {T}}{\boldsymbol {\Lambda }}_{0}(\mathbf {B} _{n}-\mathbf {B} _{0})}$$ $${\displaystyle {\boldsymbol {\nu }}_{n}={\boldsymbol {\nu }}_{0}+n}$$ $${\displaystyle \mathbf {B} _{n}=(\mathbf {X} ^{\mathsf {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0})^{-1}(\mathbf {X} ^{\mathsf {T}}\mathbf {Y} +{\boldsymbol {\Lambda }}_{0}\mathbf {B} _{0})}$$ $${\displaystyle {\boldsymbol {\Lambda }}_{n}=\mathbf {X} ^{\mathsf {T}}\mathbf {X} +{\boldsymbol {\Lambda }}_{0}}$$

• ベイズ線形回帰
• 行列正規分布

• Box, GEP ; Tiao, GC (1973). "8".統計分析におけるベイズ推論。Wiley。ISBN 0-471-57428-7。
• ガイサー、S. (1965) 「多変量解析におけるベイズ推定」数理統計年報36 ( 1): 150–159. JSTOR  2238083.
• Tiao, GC; Zellner, A. ( 1964). 「多変量回帰のベイズ推定について」。英国王立統計学会誌。シリーズ B (方法論)。26 (2): 277–285。JSTOR  2984424。