一文讲懂普通线性回归与贝叶斯线性回归的区别

将均方误差应用到线性回归问题中

首先我们先回顾频率派的解决方法，即常规解法——利用最小二乘法解决问题。损失函数采用均方误差函数。
关于模型容量问题：
在《深度学习》中专门指出，模型的容量也是一个需要用户指定的超参数，即拟合的多项式次数。此处为了简化操作，暂定为一次多项式。
解决代码如下：

import numpy as np
x = np.random.randn(10) + 5
# 初始参数
k = 2.
b = 6.
y = x * k + b
# 以均方误差形式解决线性回归
k_hat = 1
b_hat = 0
alpha = 0.01
for i in range(10000):
    # 采用梯度下降法进行迭代，更新参数k、偏置项b
    y_ = k_hat * x + b_hat
    # 均方误差函数
    loss = 1 / 2 * (y_ - y) * (y_ - y)
    grand = (y_ - y) * x
    if i % 1000 == 0:
        print("loss:" + str(np.mean(loss)))
    k_hat -= alpha * np.mean(grand)
    b_hat -= alpha * np.mean(y_ - y)
print('result:')
print(k_hat, b_hat)

输出结果如下：
loss:58.678592587395315
loss:0.29681809805254167
loss:0.11683356934778104
loss:0.04598804121481371
loss:0.018101817367916613
loss:0.007125239157084194
loss:0.0028046373473878386
loss:0.0011039616322971802
loss:0.00043454148776818417
loss:0.0001710442637384755
result:
2.0105150388524025 5.9472884150823475

可以看出模型最终拟合程度较好，损失不断下降，最终的参数与初试参数差别不大。

将最大似然估计应用到回归问题上

最小二乘法是将线性回归作为学习从输入x映射到输出$\hat{y}$的算法。现在，我们希望模型能够得到P(y$\mid$x)，而不是得到一个单独的预测$\hat{y}$。

我们定义P(y$\mid$x,w)=N(y;$\hat{y}$(x;w),${\sigma }^2$)，函数$\hat{y}$(x;k)预测高斯的均值

• 首先我们指定一个方差${\sigma }^2$
• 其次我们需要列出似然函数
• 最后求导，做最大似然估计

ps：在和用最大似然估计做解决分类问题时，有一些差别。分类问题中，我们不需要指定固定的方差。例如在逻辑回归问题中，我们仅需要通过非线性映射将概率整合进分类任务中即可。

似然函数定义如下：
$$L(\theta )=log(\prod _{i=1}^{n}P(y\mid x;\theta ))=\sum _{i=1}^{n}logP(y\mid x;\theta )=-nlog\sigma -\frac{1}{n}log(2\pi )-\sum _{i=1}^n\frac{\mid \mid {\hat{y}}^{(i)}-y^{(i)}\mid {\mid }^2}{2{\sigma }^2}$$
最大化该似然函数，等价于最小化均方误差函数，该方法其实将得到和最小二乘法类似的解法。

贝叶斯线性回归

在贝叶斯学派的观点中，我们将参数w视为一个随机变量。我们要先确定好w的先验分布。对比于最大似然估计法将P(y$\mid$x)视为正态分布求解问题，在贝叶斯线性回归下，我们将w的先验分布视为正态分布，从而进一步确定其后验分布。
贝叶斯的基本模型如下：$$y\underset{服从}{\to }N(w^{T}X,{\sigma }^{2}I)$$其中I是单位矩阵。

为什么y服从的正态分布均值为$\hat{y}$?

• 因为均值一旦为$\hat{y}$，那么最终的表达式将会和均值误差成正比，符合我们前面在极大似然估计处的推理。

为什么将$\theta$的先验分布指定为正态分布?

• 因为为了令$\theta$具有一定的广泛性，所以指定一个高熵的分布。

定义P(w)=N(w;$\mu$,${\sigma }^2$)，这样我们就确定了w的先验分布，虽然该分布可能会出错（在样本较少的情况下），但是在理论上如果样本点无限多那么我们可以得出无限收敛于正确参数的估计模型，这也是贝叶斯的基本思想。

贝叶斯公式如下：
$$P(w\mid X,y)=\frac{P(y\mid X,w)P(w\mid X)}{P(y\mid x)}=\frac{P(y\mid X,w)P(w)}{P(y)}$$
在该式子中，分母并不是我们需要关注的点，因为对于样本来说分母的值相同，且在后期计算梯度时，将会按常数处理，故我们仅仅研究分子。所以可以得出来参数w的后验分布与$P(y\mid X,w)P(w)$成正比。

接下来我们通过最大后验估计（MAP），即最大化$P(w\mid X,y)$从而得到w的最优解，直观的理解便是令我们后验概率最大化。求解过程类似于最大似然估计依旧需要取对数，构建似然函数。
$$\begin{gathered} argma{x}_{w}log(P(y\mid X,w)P(w))=argma{x}_{w}logP(y\mid X,w)+logP(w) \\ logP(y\mid X,w)+logP(w)=\log \left[\frac{1}{\sqrt{2\pi }{\sigma }_n}\exp \left(-\frac{{\left|y-X^{\top }w\right|}^2}{2{\sigma }_n^2}\right)\right]+\log \left[\frac{1}{\sqrt{2\pi }{\sigma }_w}\exp \left(-\frac{w^{\top }w}{2{\sigma }_w^2}\right)\right] \\ =-\frac{1}{2{\sigma }_n^2}{\left(y-X^{\top }w\right)}^{\top }\left(y-X^{\top }w\right)-\frac{1}{2{\sigma }_w^2}{w}^{\top }w \\ =-\frac{1}{2{\sigma }_n^2}{\Vert y-{\mathbf{X}}^{\top }w\Vert }^2-\frac{1}{2{\sigma }_w^2}\Vert w{\Vert }^2 \end{gathered}$$
很有趣，我们最终得到的式子竟然是均方误差与正则惩罚项的和（差）。这也回到了一开始的最小二乘法。
MAP为我们提供了对正则惩罚项的直观理解！所以我们可以得出结论，利用贝叶斯方法解决线性回归问题，可以有效的减少模型的过拟合程度。
由于上面的例子过于简单，维度为1。所以为了更好地举例，我们创建具有两个特征值（维度为2）的样本数据，在计算过程中将牵涉到矩阵运算。
解决代码如下：

# X代表个有10个样本值，每个样本有两个特征值
# X随机生成
X = np.random.random([10, 2]) * 10+ 2
W = np.array([
    [1],
    [2]
])
B = 1
Y = np.dot(X, W) + B
# 初始化
k_hat = np.array([
    [0.6],
    [0.5]
])
b_hat = 0.5
# 以贝叶斯解决线性回归
# 多维情况下
alpha = 0.001
lamt = 0.01 # 代表正则化惩罚程度
# 迭代一千次
for i in range(1000):
    y_ = np.dot(X, k_hat) + b_hat
    # 带正则化惩罚项的损失函数
    loss = 1 / 2 * np.dot((y_ - Y).T,  (y_ - Y)) + lamt * np.dot(k_hat.T, k_hat)
    grand = np.dot(X.T,  (y_ - Y)) + lamt * k_hat
    if i % 100 == 0:
        print("loss:" + str(np.mean(loss)))
    k_hat -= alpha * (grand)
    b_hat -= alpha * (y_ - Y)
print('result:')
print(k_hat)

输出结果如下：

loss:0.052098026507421846
loss:0.0520601388908503
loss:0.052030256081297147
loss:0.0520057579403434
loss:0.05198566785068535
loss:0.05196918640742312
loss:0.05195565917870616
loss:0.051944550414666016
loss:0.051935421524299995
loss:0.05192791345559716
result:
[[1.05673593]
[2.01810083]]

可以看出最终拟合的权重参数k_hat与初始的W基本相同，且损失值不断降低。