一、概述

线性回归的数据如下所示：

D = { ( x 1 , y 1 ) , ( x 2 , y 2 ) , ⋯   , ( x N , y N ) } x i ∈ R p , y i ∈ R , i = 1 , 2 , ⋯   , N X = ( x 1 , x 1 , ⋯   , x N ) T = ( x 1 T x 2 T ⋮ x N T ) = ( x 11 x 12 ⋯ x 1 p x 21 x 22 ⋯ x 2 p ⋮ ⋮ ⋱ ⋮ x N 1 x N 2 ⋯ x N p ) N × p Y = ( y 1 y 2 ⋮ y N ) N × 1 D=\left \{(x_{1},y_{1}),(x_{2},y_{2}),\cdots ,(x_{N},y_{N})\right \}\\ x_{i}\in \mathbb{R}^{p},y_{i}\in \mathbb{R},i=1,2,\cdots ,N\\ X=(x_{1},x_{1},\cdots ,x_{N})^{T}=\begin{pmatrix} x_{1}^{T}\\ x_{2}^{T}\\ \vdots \\ x_{N}^{T} \end{pmatrix}=\begin{pmatrix} x_{11} & x_{12} & \cdots &x_{1p} \\ x_{21} & x_{22}& \cdots &x_{2p} \\ \vdots & \vdots & \ddots &\vdots \\ x_{N1}& x_{N2} & \cdots & x_{Np} \end{pmatrix}_{N \times p}\\ Y=\begin{pmatrix} y_{1}\\ y_{2}\\ \vdots \\ y_{N} \end{pmatrix}_{N \times 1} D={(x1​,y1​),(x2​,y2​),⋯,(xN​,yN​)}xi​∈Rp,yi​∈R,i=1,2,⋯,NX=(x1​,x1​,⋯,xN​)T=⎝⎜⎜⎜⎛​x1T​x2T​⋮xNT​​⎠⎟⎟⎟⎞​=⎝⎜⎜⎜⎛​x11​x21​⋮xN1​​x12​x22​⋮xN2​​⋯⋯⋱⋯​x1p​x2p​⋮xNp​​⎠⎟⎟⎟⎞​N×p​Y=⎝⎜⎜⎜⎛​y1​y2​⋮yN​​⎠⎟⎟⎟⎞​N×1​

对于线性回归问题，假设：

f ( x ) = w T x y = f ( x ) + ε ε ∼ N ( 0 , σ 2 ) f(x)=w^{T}x\\ y=f(x)+\varepsilon \\ \varepsilon \sim N(0,\sigma ^{2}) f(x)=wTxy=f(x)+εε∼N(0,σ2)

其中 w w w是参数， ε \varepsilon ε是噪声。以二维的数据为例，其图像如下：

从几何角度出发求解线性回归问题，可以使用最小二乘估计（Least Square Estimate，LSE），而从概率角度来看也就是噪声为高斯分布的极大似然估计（Maximum Likelihood Estimate，MLE），公式为：

w M L E = a r g m a x w    P ( D a t a ∣ w ) w_{MLE}=\underset{w}{argmax}\; P(Data|w) wMLE​=wargmax​P(Data∣w)

为了防止过拟合的问题，采用了正则化的最小二乘估计的方法，包括Lasso正则化和Ridge正则化两种方式，这种方法从概率角度来看也就是最大后验估计（Maximum a Posteriori Estimation，MAP），其噪声也是高斯分布，这里参数 w w w的后验满足：

P ( w ∣ D a t a ) ∝ P ( D a t a ∣ w ) ⋅ P ( w ) P(w|Data)\propto P(Data|w)\cdot P(w) P(w∣Data)∝P(Data∣w)⋅P(w)

因此最大后验估计的公式也就是：

w M A P = a r g m a x w    P ( w ∣ D a t a ) = a r g m a x w    P ( D a t a ∣ w ) ⋅ P ( w ) w_{MAP}=\underset{w}{argmax}\; P(w|Data)=\underset{w}{argmax}\; P(Data|w)\cdot {\color{Red}{P(w)}} wMAP​=wargmax​P(w∣Data)=wargmax​P(Data∣w)⋅P(w)

这里的 P ( w ) P(w) P(w)也就是参数 w w w的先验，如果 P ( w ) P(w) P(w)服从高斯分布，则对应Lasso正则化，如果 P ( w ) P(w) P(w)服从拉普拉斯分布，则对应Ridge正则化。

对于上述方法的详细解释，可以参考链接：线性回归|机器学习推导系列（三） 。

上述方法都是点估计方法，而在贝叶斯线性回归中我们希望利用贝叶斯学派的方法来求解出后验概率分布 P ( w ∣ D a t a ) P(w|Data) P(w∣Data)。

二、推断与预测

1. 推断

贝叶斯线性回归的推断（Inference）任务也就是求解参数 w w w的后验分布 P ( w ∣ D a t a ) P(w|Data) P(w∣Data)。首先，进行以下变换：

P ( w ∣ D a t a ) = P ( w ∣ X , Y ) = P ( w , Y ∣ X ) P ( Y ∣ X ) = P ( Y ∣ w , X ) ⋅ P ( w ∣ X ) ∫ P ( Y ∣ w , X ) ⋅ P ( w ∣ X ) d w P(w|Data)=P(w|X,Y)\\ =\frac{P(w,Y|X)}{P(Y|X)}\\ =\frac{{\color{Red}{P(Y|w,X)}}\cdot {\color{Blue}{P(w|X)}}}{\int P(Y|w,X)\cdot P(w|X)\mathrm{d}w} P(w∣Data)=P(w∣X,Y)=P(Y∣X)P(w,Y∣X)​=∫P(Y∣w,X)⋅P(w∣X)dwP(Y∣w,X)⋅P(w∣X)​

对于上面得到的式子，分子中 P ( Y ∣ w , X ) P(Y|w,X) P(Y∣w,X)相当于似然， P ( w ∣ X ) P(w|X) P(w∣X)也就是先验 P ( w ) P(w) P(w)（因为 w w w和 X X X没有关系），另外分母是一个常数，因此我们可以得到：

P ( w ∣ D a t a ) ∝ P ( Y ∣ w , X ) ⋅ P ( w ) P(w|Data)\propto P(Y|w,X)\cdot P(w) P(w∣Data)∝P(Y∣w,X)⋅P(w)

对于似然 P ( Y ∣ w , X ) P(Y|w,X) P(Y∣w,X)，由于之前定义的 y y y和 x x x之间的关系是一个线性高斯模型，因此可以得到：

P ( Y ∣ w , X ) = ∏ i = 1 N P ( y i ∣ w , x i ) = ∏ i = 1 N N ( y i ∣ w T x i , σ 2 ) P(Y|w,X)=\prod_{i=1}^{N}P(y_{i}|w,x_{i})=\prod_{i=1}^{N}N(y_{i}|w^{T}x_{i},\sigma ^{2}) P(Y∣w,X)=i=1∏N​P(yi​∣w,xi​)=i=1∏N​N(yi​∣wTxi​,σ2)

对于先验 P ( w ) P(w) P(w)，是人为指定的，只要合理就可以任意给定这个分布，通常我们选择高斯分布，注意这里是一个多维的高斯分布，我们选择的高斯分布均值为 0 0 0，方差矩阵用 Σ p \Sigma _{p} Σp​来表示：

P ( w ) = N ( 0 , Σ p ) P(w)=N(0,\Sigma _{p}) P(w)=N(0,Σp​)

因此也就得到：

P ( w ∣ D a t a ) ∝ ∏ i = 1 N N ( y i ∣ w T x i , σ 2 ) ⋅ N ( 0 , Σ p ) P(w|Data)\propto \prod_{i=1}^{N}N(y_{i}|w^{T}x_{i},\sigma ^{2})\cdot N(0,\Sigma _{p}) P(w∣Data)∝i=1∏N​N(yi​∣wTxi​,σ2)⋅N(0,Σp​)

由于似然 P ( Y ∣ w , X ) P(Y|w,X) P(Y∣w,X)和先验 P ( w ) P(w) P(w)都是高斯分布，利用共轭的知识（高斯分布是自共轭的），我们也可以得到 P ( w ∣ D a t a ) P(w|Data) P(w∣Data)也是一个高斯分布，推断任务也就是求解这个高斯分布的均值和方差，我们用 μ w \mu _{w} μw​和 Σ w \Sigma _{w} Σw​来表示这个高斯分布的均值和方差：

P ( w ∣ D a t a ) ∼ N ( μ w , Σ w ) P(w|Data)\sim N(\mu _{w},\Sigma _{w}) P(w∣Data)∼N(μw​,Σw​)

将似然 P ( Y ∣ w , X ) P(Y|w,X) P(Y∣w,X)做进一步的变换，我们可以得到这个概率的均值和方差：

P ( Y ∣ w , X ) = ∏ i = 1 N 1 ( 2 π ) 1 / 2 σ e x p { − 1 2 σ 2 ( y i − w T x i ) 2 } = 1 ( 2 π ) N / 2 σ N e x p { − 1 2 σ 2 ∑ i = 1 N ( y i − w T x i ) 2 } P(Y|w,X)=\prod_{i=1}^{N}\frac{1}{(2\pi )^{1/2}\sigma }exp\left \{-\frac{1}{2\sigma ^{2}}(y_{i}-w^{T}x_{i})^{2}\right \}\\ =\frac{1}{(2\pi )^{N/2}\sigma ^{N}}exp\left \{-\frac{1}{2\sigma ^{2}}{\color{Red}{\sum _{i=1}^{N}(y_{i}-w^{T}x_{i})^{2}}}\right \} P(Y∣w,X)=i=1∏N​(2π)1/2σ1​exp{−2σ21​(yi​−wTxi​)2}=(2π)N/2σN1​exp{−2σ21​i=1∑N​(yi​−wTxi​)2}

上式中红色的部分可以写成矩阵相乘的形式：

∑ i = 1 N ( y i − w T x i ) 2 = ( y 1 − w T x 1 y 2 − w T x 2 ⋯ y N − w T x N ) ⏟ Y T − w T X T ( y 1 − w T x 1 y 2 − w T x 2 ⋮ y N − w T x N ) = ( Y T − w T X T ) ( Y − X w ) = ( Y − X w ) T ( Y − X w ) \sum _{i=1}^{N}(y_{i}-w^{T}x_{i})^{2}=\underset{Y^{T}-w^{T}X^{T}}{\underbrace{\begin{pmatrix} y_{1}-w^{T}x_{1} & y_{2}-w^{T}x_{2} & \cdots & y_{N}-w^{T}x_{N} \end{pmatrix}}}\begin{pmatrix} y_{1}-w^{T}x_{1}\\ y_{2}-w^{T}x_{2}\\ \vdots \\ y_{N}-w^{T}x_{N} \end{pmatrix}\\ =(Y^{T}-w^{T}X^{T})(Y-Xw)\\ =(Y-Xw)^{T}(Y-Xw) i=1∑N​(yi​−wTxi​)2=YT−wTXT (y1​−wTx1​​y2​−wTx2​​⋯​yN​−wTxN​​)​​⎝⎜⎜⎜⎛​y1​−wTx1​y2​−wTx2​⋮yN​−wTxN​​⎠⎟⎟⎟⎞​=(YT−wTXT)(Y−Xw)=(Y−Xw)T(Y−Xw)

因此也就有：

P ( Y ∣ w , X ) = 1 ( 2 π ) N / 2 σ N e x p { − 1 2 σ 2 ∑ i = 1 N ( y i − w T x i ) 2 } = 1 ( 2 π ) N / 2 σ N e x p { − 1 2 σ 2 ( Y − X w ) T ( Y − X w ) } = 1 ( 2 π ) N / 2 σ N e x p { − 1 2 ( Y − X w ) T σ − 2 I ( Y − X w ) } = N ( X w , σ 2 I ) P(Y|w,X)=\frac{1}{(2\pi )^{N/2}\sigma ^{N}}exp\left \{-\frac{1}{2\sigma ^{2}}\sum _{i=1}^{N}(y_{i}-w^{T}x_{i})^{2}\right \}\\ =\frac{1}{(2\pi )^{N/2}\sigma ^{N}}exp\left \{-\frac{1}{2\sigma ^{2}}(Y-Xw)^{T}(Y-Xw)\right \}\\ =\frac{1}{(2\pi )^{N/2}\sigma ^{N}}exp\left \{-\frac{1}{2}(Y-Xw)^{T}\sigma ^{-2}I(Y-Xw)\right \}\\ =N(Xw,\sigma ^{2}I) P(Y∣w,X)=(2π)N/2σN1​exp{−2σ21​i=1∑N​(yi​−wTxi​)2}=(2π)N/2σN1​exp{−2σ21​(Y−Xw)T(Y−Xw)}=(2π)N/2σN1​exp{−21​(Y−Xw)Tσ−2I(Y−Xw)}=N(Xw,σ2I)

至此似然 P ( Y ∣ w , X ) P(Y|w,X) P(Y∣w,X)就表示成了标准的高斯分布的形式，代入后验中我们可以得到：

P ( w ∣ D a t a ) ∝ N ( X w , σ 2 I ) ⋅ N ( 0 , Σ p ) ∝ e x p { − 1 2 ( Y − X w ) T σ − 2 I ( Y − X w ) } ⋅ e x p { − 1 2 w T Σ p − 1 w } = e x p { − 1 2 σ 2 ( Y − X w ) T ( Y − X w ) − 1 2 w T Σ p − 1 w } = e x p { − 1 2 σ 2 ( Y T Y − 2 Y T X w + w T X T X w ) − 1 2 w T Σ p − 1 w } P(w|Data)\propto N(Xw,\sigma ^{2}I)\cdot N(0,\Sigma _{p})\\ \propto exp\left \{-\frac{1}{2}(Y-Xw)^{T}\sigma ^{-2}I(Y-Xw)\right \}\cdot exp\left \{-\frac{1}{2}w^{T}\Sigma _{p}^{-1}w\right \}\\ =exp\left \{-\frac{1}{2\sigma ^{2}}(Y-Xw)^{T}(Y-Xw)-\frac{1}{2}w^{T}\Sigma _{p}^{-1}w\right \}\\ =exp\left \{-\frac{1}{2\sigma ^{2}}(Y^{T}Y-2Y^{T}Xw+w^{T}X^{T}Xw)-\frac{1}{2}w^{T}\Sigma _{p}^{-1}w\right \} P(w∣Data)∝N(Xw,σ2I)⋅N(0,Σp​)∝exp{−21​(Y−Xw)Tσ−2I(Y−Xw)}⋅exp{−21​wTΣp−1​w}=exp{−2σ21​(Y−Xw)T(Y−Xw)−21​wTΣp−1​w}=exp{−2σ21​(YTY−2YTXw+wTXTXw)−21​wTΣp−1​w}

对于上面的式子我们可以写出它的指数的关于 w w w的二次项和一次项：

二 次 项 ： − 1 2 σ 2 w T X T X w − 1 2 w T Σ p − 1 w = − 1 2 w T ( σ − 2 X T X + Σ p − 1 ) w 一 次 项 ： − 1 2 σ 2 ( − 2 Y T X w ) = σ − 2 Y T X w 二次项：-\frac{1}{2\sigma ^{2}}w^{T}X^{T}Xw-\frac{1}{2}w^{T}\Sigma _{p}^{-1}w=-\frac{1}{2}w^{T}(\sigma ^{-2}X^{T}X+\Sigma _{p}^{-1})w\\ 一次项：-\frac{1}{2\sigma ^{2}}(-2Y^{T}Xw)=\sigma ^{-2}Y^{T}Xw 二次项：−2σ21​wTXTXw−21​wTΣp−1​w=−21​wT(σ−2XTX+Σp−1​)w一次项：−2σ21​(−2YTXw)=σ−2YTXw

对于一个多维的高斯分布的指数部分，将其展开可以得到：

e x p { − 1 2 ( X − μ ) T Σ − 1 ( X − μ ) } = e x p { − 1 2 ( X T Σ − 1 − μ T Σ − 1 ) ( X − μ ) } = e x p { − 1 2 ( X T Σ − 1 X ⏟ 二 次 项 − 2 μ T Σ − 1 X ⏟ 一 次 项 + μ T Σ − 1 μ ) } exp\left \{-\frac{1}{2}(X-\mu)^{T}\Sigma ^{-1}(X-\mu)\right \}\\ =exp\left \{-\frac{1}{2}(X^{T}\Sigma ^{-1}-\mu^{T}\Sigma ^{-1})(X-\mu)\right \}\\ =exp\left \{-\frac{1}{2}(\underset{二次项}{\underbrace{X^{T}\Sigma ^{-1}X}}\underset{一次项}{\underbrace{-2\mu^{T}\Sigma ^{-1}X}}+\mu^{T}\Sigma ^{-1}\mu)\right \} exp{−21​(X−μ)TΣ−1(X−μ)}=exp{−21​(XTΣ−1−μTΣ−1)(X−μ)}=exp⎩⎨⎧​−21​(二次项 XTΣ−1X​​一次项 −2μTΣ−1X​​+μTΣ−1μ)⎭⎬⎫​

而这里的有关于 X X X的二次项和一次项为：

二 次 项 ： − 1 2 X T Σ − 1 X 一 次 项 ： μ T Σ − 1 X 二次项：-\frac{1}{2}X^{T}\Sigma ^{-1}X\\ 一次项：\mu^{T}\Sigma ^{-1}X 二次项：−21​XTΣ−1X一次项：μTΣ−1X

通过二次项也就可以得到后验 P ( w ∣ D a t a ) P(w|Data) P(w∣Data)的精度矩阵 Λ w \Lambda _{w} Λw​（这是一个对称矩阵）：

Λ w = Σ w − 1 = σ − 2 X T X + Σ p − 1 \Lambda _{w}=\Sigma _{w}^{-1}=\sigma ^{-2}X^{T}X+\Sigma _{p}^{-1} Λw​=Σw−1​=σ−2XTX+Σp−1​

由一次项可以得出：

μ w T Σ w − 1 = μ w T Λ w = σ − 2 Y T X ⇒ Λ w μ w = σ − 2 X T Y ⇒ μ w = σ − 2 Λ w − 1 X T Y \mu _{w}^{T}\Sigma _{w}^{-1}=\mu _{w}^{T}\Lambda _{w}=\sigma ^{-2}Y^{T}X\\ \Rightarrow \Lambda _{w}\mu _{w}=\sigma ^{-2}X^{T}Y\\ \Rightarrow \mu _{w}=\sigma ^{-2}\Lambda _{w}^{-1}X^{T}Y μwT​Σw−1​=μwT​Λw​=σ−2YTX⇒Λw​μw​=σ−2XTY⇒μw​=σ−2Λw−1​XTY

如此我们最终就求得了后验 P ( w ∣ D a t a ) P(w|Data) P(w∣Data)的均值和方差：

μ w = σ − 2 Λ w − 1 X T Y Σ w = Λ w − 1 其 中 Λ w = σ − 2 X T X + Σ p − 1 \mu _{w}=\sigma ^{-2}\Lambda _{w}^{-1}X^{T}Y\\ \Sigma _{w}=\Lambda _{w}^{-1}\\ 其中\Lambda _{w}=\sigma ^{-2}X^{T}X+\Sigma _{p}^{-1} μw​=σ−2Λw−1​XTYΣw​=Λw−1​其中Λw​=σ−2XTX+Σp−1​

2. 预测

预测（Prediction）问题是利用已经求解出来的贝叶斯线性回归模型，在给定 x ∗ x^{*} x∗的前提下，求解 y ∗ y^{*} y∗。要注意，由于这里的 w w w服从后验分布 P ( w ∣ D a t a ) P(w|Data) P(w∣Data)，因此需要求解的是分布 P ( y ∗ ∣ D a t a , x ∗ ) P(y^{*}|Data,x^{*}) P(y∗∣Data,x∗)。

首先对于 f ( x ∗ ) f(x^{*}) f(x∗)，有 f ( x ∗ ) = ( x ∗ ) T w f(x^{*})=(x^{*})^{T}w f(x∗)=(x∗)Tw，而 w w w服从后验分布 N ( μ w , Σ w ) N(\mu _{w},\Sigma _{w}) N(μw​,Σw​)， x ∗ x^{*} x∗可以看做一个常数，因此就有：

f ( x ∗ ) = ( x ∗ ) T w ∼ N ( ( x ∗ ) T μ w , ( x ∗ ) T Σ w x ∗ ) f(x^{*})=(x^{*})^{T}w\sim N((x^{*})^{T}\mu _{w},(x^{*})^{T}\Sigma _{w}x^{*}) f(x∗)=(x∗)Tw∼N((x∗)Tμw​,(x∗)TΣw​x∗)

另外有 y ∗ = f ( x ∗ ) + ε y^{*}=f(x^{*})+\varepsilon y∗=f(x∗)+ε， ε ∼ N ( 0 , σ 2 ) \varepsilon \sim N(0,\sigma ^{2}) ε∼N(0,σ2)，因此就有（两个高斯分布随机变量相加后的概率分布相当于将两个高斯分布的均值和方差均相加得到的高斯分布）：

P ( y ∗ ∣ D a t a , x ∗ ) = N ( ( x ∗ ) T μ w , ( x ∗ ) T Σ w x ∗ + σ 2 ) P(y^{*}|Data,x^{*})=N((x^{*})^{T}\mu _{w},(x^{*})^{T}\Sigma _{w}x^{*}+\sigma ^{2}) P(y∗∣Data,x∗)=N((x∗)Tμw​,(x∗)TΣw​x∗+σ2)

综合一下，具体的过程为：

P ( y ∗ ∣ D a t a , x ∗ ) = ∫ w P ( y ∗ ∣ w , X , Y , x ∗ ) P ( w ∣ X , Y , x ∗ ) d w = ∫ w P ( y ∗ ∣ w , x ∗ ) P ( w ∣ X , Y ) d w = N ( ( x ∗ ) T μ w , ( x ∗ ) T Σ w x ∗ + σ 2 ) P(y^{*}|Data,x^{*})=\int _{w}P(y^{*}|w,X,Y,x^{*})P(w|X,Y,x^{*})\mathrm{d}w\\ =\int _{w}P(y^{*}|w,x^{*})P(w|X,Y)\mathrm{d}w\\ =N((x^{*})^{T}\mu _{w},(x^{*})^{T}\Sigma _{w}x^{*}+\sigma ^{2}) P(y∗∣Data,x∗)=∫w​P(y∗∣w,X,Y,x∗)P(w∣X,Y,x∗)dw=∫w​P(y∗∣w,x∗)P(w∣X,Y)dw=N((x∗)Tμw​,(x∗)TΣw​x∗+σ2)