贝叶斯线性回归（Bayesian Linear Regression）

标签（空格分隔）： 监督学习

---

@ author : duanxxnj@163.com
@ time : 2015-06-19

原文地址

本文的研究顺序是：

关于参数估计

在很多的机器学习或数据挖掘的问题中，我们所面对的只有数据，但数据中潜在的概率密度函数是不知道的，其概率密度分布需要我们从数据中估计出来。想要确定数据对应的概率密度分布，就需要确定两个东西：概率密度函数的形式 和 概率密度函数的参数。

有时可能知道的是概率密度函数的形式(高斯、瑞利等等)，但是不知道具体的参数，例如均值或者方差；还有的时候可能不知道概率密度的类型，但是知道一些估计的参数，比如均值和方差。

关于上面提到的需要确定的两个东西：概率密度函数的形式和参数，至少在机器学习的教课书上，我所看到的情况都是：给了一堆数据，然后假设其概率密度函数的形式为 高斯分布 ，或者是混合高斯分布，那么，剩下的事情就是对高斯分布的参数，$\mu$ 和 $\sigma^2$ 进行估计。所以，参数估计，便成了极其最重要的问题。

其实，常用的参数估计方法有：极大似然估计、最大后验估计、贝叶斯估计、最大熵估计、混合模型估计。他们之间是有递进关系的，想要理解后一个参数估计方法，最好对前一个参数估计有足够的理解。

要想清晰的说明贝叶斯线性回归，或者叫做贝叶斯参数估计，就必须对极大似然估计、最大后验估计做详细的说明，他们之间是有递进的关系的。

极大似然估计

在之前《多项式回归》的文章中，用最后一小节是线性回归的概率解释，其中就说明了以平方误差维损失函数的最小二乘法和极大似然估计的等价性，在这个基础上，本文更为详细的讨论极大似然估计。

这里先以一个分类问题来说明一般参数估计面对的数据形式。考虑一个$M$类的问题，特征向量服从$p(x|\omega_i),i=1,2...,M$ 分布。这是现实情况中最常见的一种数据存在形式，数据集合$X$是由$M$个类别的数据子集$X_m,m=1,2...,M$ 组成的，第$m$类别的数据子集$X_m$ 对应的概率密度函数是$p(x|\omega_m)$。

前面已经介绍过了，想要确定数据的概率分布，需要知道概率密度函数的 形式 和 参数，这里首先做一个基本假设：概率分布的形式已知，比如假设每个类别的数据都满足高斯分布，那么，似然函数就可以以参数 $\textbf{ $\theta_i$ }$的形式表示，如果是高斯分布，则参数为$\mu_i$和$\sigma_i^2$，即$\theta_i=(\mu_i,\sigma_i^2)$。

为了强调概率分布$p(x|\omega_i)$和 $\theta_i$ 有关，将对应的概率密度函数记为$p(x|\omega_i;\theta_i)$，这种记法属于频率概率学派的记法。这里的极大似然估计对应于一个类条件概率密度函数。

在概率论中一直有两大学派，分别是频率学派和贝叶斯学派。简单点说，频率学派认为，概率是频率的极限，比如投硬币，当实验次数足够大时，正面朝上的频率可以认为是这枚硬币正面朝上的概率，这个是频率学派。但是，如果要预测一些未发生过的事情，比如，北极的冰山在2050年完全融化的概率，由于这个事情完全没有发生过，所以无法用频率来代替概率表示，只能研究过去几十年，北极冰山融化的速率，并将其作为先验条件，来预测北极的冰山在2050年完全融化的概率，这就是概率的贝叶斯学派。上面的问题，如果用贝叶斯学派的记法的话，是：$p(x|\omega_i,\theta_i)$。这两个学派适用的情况不太一样，但是，在我目前所用到的概率论的知识中，貌似这两个学派并没有什么太大的区别，只是记法略有不同，稍微注意下即可。

从上面的描述中可以知道，利用每一个类$X_i$中已知的特征向量集合，可以估计出其对应的参数$\theta_i$。进一步假设每一类中的数据不影响其他类别数据的参数估计，那么上面的$M$个类别的参数估计就可以用下面这个统一的模型，独立的解决：

设$x_1,x_2,...,x_N$ 是从概率密度函数$p(x;\theta)$中随机抽取的样本，那么就可以得到联合概率密度函数 $p(X;\theta)$， 其中$X=\{x_1,x_2,...,x_N\}$ 是样本集合。假设不同的样本之间具有统计独立性，那么：

$$p(X;\theta)\equiv p(x_1,x_2,...,x_N;\theta)=\prod_{k=1}^Np(x_k;\theta)$$

注意：这里的$p(x_k;\theta)$ 本来的写法是 $p(x|\omega_i;\theta_i)$ , 是一个类条件概率密度函数，只是因为这里是一个统一的模型，所以可以将$w_i$ 省略。

需要重申一下，想要得到上面这个公式，是做了几个基本的假设的，第一：假设$M$个类别的数据子集的概率密度函数形式一样，只是参数的取值不同；第二：假设类别$i$中的数据和类别$j$中的数据是相互独立抽样的，即类别$j$的参数仅仅根据类别$j$的数据就可以估计出来，类别$i$的数据并不能影响类别$j$的参数估计，反之亦然；第三：每个类别内的样本之间具有统计独立性，即每个类别内的样本之间是独立同分布 ($iid$) 的。

此时，就可以使用最大似然估计(Maximum Likelihood,ML)来估计参数$\theta$了：

$$\hat{\theta}_{ML}=arg\max_{\theta}\prod_{k=1}^Np(x_k;\theta)$$

为了得到最大值，$\hat{\theta}_{ML}$必须满足的必要条件是,似然函数对$\theta$的梯度必须为0，即：

$$\frac{\partial \prod_{k=1}^Np(x_k;\theta)}{\partial\theta}=0$$

一般我们取其对数形式：

$$L(\theta)\equiv ln\prod_{k=1}^Np(x_k;\theta)$$

$$\frac{\partial L(\theta)}{\partial \theta}=\sum_{k=1}^N \frac{\partial ln p(x_k;\theta)}{\partial \theta}=\sum_{k=1}^N\frac{1}{p(x_k;\theta)} \frac{\partial p(x_k;\theta)}{\partial \theta}=0$$

需要注意：极大似然估计对应于似然函数的峰值

极大似然估计有两个非常重要的性质：渐进无偏 和 渐进一致性，有了这两个性质，使得极大似然估计的成为了非常简单而且实用的参数估计方法。这里假设$\theta_0$是密度函数$p(x;\theta)$中未知参数的准确值。

渐进无偏

极大似然估计是渐进无偏的，即：

$$\lim_{N \to \infty}E[\hat{\theta}_{ML}]=\theta_0$$

也就是说，这里认为估计值 $\hat{\theta}_{ML}$ 本身是一个随机变量（因为不同的样本集合$X$会得到不同的 $\hat{\theta}_{ML}$），那么其均值就是未知参数的真实值，这就是渐进无偏。

渐进一致

极大似然估计是渐进一致的，即：

$$\lim_{N \to \infty}prob\{ \lVert \hat{\theta}_{ML}- \theta_0 \rVert \leqslant \epsilon\} = 1$$

这个公式还可以表示为：

$$\lim_{N \to \infty} E \lVert \hat{\theta}_{ML}- \theta_0 \rVert^2 = 0$$

对于一个估计器而言，一致性是非常重要的，因为存在满足无偏性，但是不满足一致性的情况，比如，$\hat{\theta}_{ML}$在 $\theta_0$周围震荡。如果不满足一致性，那么就会出现很大的方差。

注意：以上两个性质，都是在渐进的前提下（$N \to \infty$）才能讨论的，即只有$N$足够大时，上面两个性质才能成立

最大后验估计

在最大似然估计（MAP）中，将$\theta$看做是未知的参数，说的通俗一点，最大似然估计是$\theta$的函数，其求解过程就是找到使得最大似然函数最大的那个参数$\theta$。

从最大后验估计开始，将参数$\theta$ 看成一个随机变量，并在已知样本集$\{x_1,x_2,...,x_N\}$ 的条件下，估计参数$\theta$。

这里一定要注意，在最大似然估计中，参数$\theta$是一个定值，只是这个值未知，最大似然函数是$\theta$的函数，这里$\theta$是没有概率意义的。但是，在最大后验估计中，$\theta$是有概率意义的，$\theta$有自己的分布，而这个分布函数，需要通过已有的样本集合$X$得到，即最大后验估计需要计算的是 $p(\theta | X)$

根据贝叶斯理论：

$$p(\theta|X)=\frac{p(\theta)p(X|\theta)}{p(X)}$$

这就是参数$\theta$关于已有数据集合$X$的后验概率，要使得这个后验概率最大，和极大似然估计一样，这里需要对后验概率函数求导。由于分子中的$p(X)$相对于$\theta$是独立的，随意可以直接忽略掉$p(X)$。

$$\hat{\theta}_{MAP}=arg\max_{\theta}p(\theta|X)=arg\max_{\theta}p(\theta)p(X|\theta)$$

为了得到参数$\theta$，和ML一样，需要对$p(\theta|X)$求梯度，并使其等于0：

$$\frac{p(\theta|X)}{\partial\theta}=\frac{p(\theta)p(X|\theta）}{\partial\theta}=0$$

注意：这里$p(X|\theta)$和极大似然估计中的似然函数$p(X;\theta)$是一样的，只是记法不一样。MAP和ML的区别是：MAP是在ML的基础上加上了$p(\theta)$

这里需要说明，虽然从公式上来看 $MAP=ML*p(\theta)$，但是这两种算法有本质的区别，ML将$\theta$视为一个确定未知的值，而MAP则将$\theta$视为一个随机变量。

在MAP中，$p(\theta)$称为$\theta$的先验，假设其服从均匀分布，即对于所有$\theta$取值，$p(\theta)$都是同一个常量，则MAP和ML会得到相同的结果。当然了，如果$p(\theta)$的方差非常的小，也就是说，$p(\theta)$是近似均匀分布的话，MAP和ML的结果自然也会非常的相似。

贝叶斯估计

注意：以下所有的概率分布表述方式均为贝叶斯学派的表述方式。

贝叶斯估计核心问题

为了防止标号混淆，这里定义已有的样本集合为$D$，而不是之前的$X$。样本集合$D$中的样本都是从一个 固定但是未知 的概率密度函数$p(x)$中独立抽取出来的，要求根据这些样本估计$x$的概率分布，记为$p(x|D)$，并且使得$p(x|D)$尽量的接近$p(x)$，这就是贝叶斯估计的核心问题。

贝叶斯估计第一个重要元素

虽然$p(x)$是未知的，但是前面提到过，一个密度分布的两个要素为：形式和参数，我们可以假设$p(x)$的形式已知，但是参数$\theta$的取值未知。这里就有了贝叶斯估计的第一个重要元素 $p(x|\theta)$，这是一个条件概率密度函数，准确的说，是一个类条件概率密度函数（具体原因参见本文前面关于极大似然估计的说明）。强调一下：$p(x|\theta)$的形式是已知的，只是参数$\theta$的取值未知。由于这里的$x$可以看成一个测试样本，所以这个条件密度函数，从本质上讲，是$\theta$ 在点 $x$ 处的似然估计。

贝叶斯估计第二个重要元素

由于参数$\theta$的取值未知，且，我们将$\theta$看成一个随机变量，那么，在观察到具体的训练样本之前，关于$\theta$的全部知识，可以用一个先验概率密度函数$p(\theta)$来表示，对于训练样本的观察，使得我们能够把这个先验概率密度转化成为后验概率密度函数$p(\theta|D)$，根据后验概率密度相关的论述知道，我们希望$p(\theta|D)$在$\theta$的真实值附近有非常显著的尖峰。这里的这个后验概率密度，就是贝叶斯估计的第二个主要元素。

现在，将贝叶斯估计核心问题$p(x|D)$，和贝叶斯估计的两个重要元素：$p(x|\theta)$、$p(\theta|D)$联系起来：

$$p(x|D)=\int p(x,\theta|D) d\theta=\int p(x|\theta,D)p(\theta|D)d\theta$$

上面式子中，$x$是测试样本，$D$是训练集，$x$和$D$的选取是独立进行的，因此，$p(x|\theta,D)$可以写成$p(x|\theta)$。所以，贝叶斯估计的核心问题就是下面这个公式：

$$p(x|D)=\int p(x|\theta)p(\theta|D)d\theta$$

下面这句话一定要理解：这里$p(x|\theta)$是$\theta$关于测试样本$x$这一个点的似然估计，而$p(\theta|D)$则是$\theta$在已有样本集合上的后验概率。这就是为什么本文一开始并没有直接讲贝叶斯估计，而是先说明极大似然估计和最大后验估计的原因。其中，后验概率$p(\theta|D)$为：

$$p(\theta|D)=\frac{p(D|\theta)p(\theta)}{p(D)}=\frac{p(D|\theta)p(\theta)}{\int p(D|\theta)p(\theta)d\theta}$$

$$p(D|\theta)=\prod_{k=1}^N p(x_k|\theta)$$

上面这个式子就是贝叶斯估计最核心的公式，它把类条件概率密度$p(x|D)$ (这里一定要理解为什么是类条件概率密度，其实这个的准确写法可以是$p(x|D_m)$,或者 $p(x|w_m,D_m)$，具体原因参见本文前面关于极大似然估计的部分) 和未知参数向量$\theta$的后验概率密度$p(\theta|D)$联系在了一起。如果后延概率密度$p(\theta|D)$在某一个值$\hat{\theta}$附近形成显著的尖峰，那么就有$p(x|D)\approx p(x|\hat{\theta})$，就是说，可以用估计值 $\hat{\theta}$近似代替真实值所得的结果。

贝叶斯估计的增量学习

为了明确的表示样本集合$D$中有$n$个样本，这里采用记号：$D^n=\{x_1,x_2,...,x_n\}$。根据前一个公式，在$n > 1$的情况下有：

$$p(D^n|\theta)=p(x_n|\theta)p(D^{n-1}|\theta)$$

可以很容易得到：

$$p(\theta|D^n)=\frac{p(x_n|\theta)p(D^{n-1}|\theta)p(\theta)}{\int p(x_n|\theta)p(D^{n-1}|\theta)p(\theta)d\theta}=\frac{p(x_n|\theta)p(\theta |D^{n-1})}{\int p(x_n|\theta)p(\theta |D^{n-1})d\theta}$$

当没有观测样本时，定义$p(\theta|D^0)=p(\theta)$，为参数$\theta$的初始估计。然后让样本集合依次进入上述公式，就可以得到一系列的概率密度函数：$p(\theta|D^0)$、$p(\theta|D^1)$、$p(\theta|D^2)$、…、$p(\theta|D^n)$，这一过程称为参数估计贝叶斯递归法，也叫贝叶斯估计的增量学习。这是一个在线学习算法，它和随机梯度下降法有很多相似之处。

贝叶斯线性回归

根据之前的文章《线性回归》、《多项式回归》中关于极大似然估计的说明，以及本文前面关于极大似然估计的论述，可以很容易知道，如果要将极大似然估计应用到线性回归模型中，模型的复杂度会被两个因素所控制：基函数的数目和样本的数目。尽管为对数极大似然估计加上一个正则项（或者是参数的先验分布），在一定程度上可以限制模型的复杂度，防止过拟合，但基函数的选择对模型的性能仍然起着决定性的作用。

上面说了那么大一段，就是想说明一个问题：由于极大似然估计总是会使得模型过于的复杂以至于产生过拟合的现象，所以单纯的适用极大似然估计并不是特别的有效。

当然，交叉验证是一种有效的限制模型复杂度，防止过拟合的方法，但是交叉验证需要将数据分为训练集合测试集，对数据样本的浪费也是非常的严重的。

基于上面的讨论，这里就可以引出本文的核心内容：贝叶斯线性回归。贝叶斯线性回归不仅可以解决极大似然估计中存在的过拟合的问题，而且，它对数据样本的利用率是100%，仅仅使用训练样本就可以有效而准确的确定模型的复杂度。

这里面对的模型是线性回归模型，其详细的介绍可以参见前面的文章《线性回归》,线性回归模型是一组输入变量$x$的基函数的线性组合，在数学上其形式如下：

$$y(x,w)=w_{0}+\sum_{j=1}^{M}\omega_{j}\phi_{j}(x)$$

这里$\phi_{j}(x)$就是前面提到的基函数，总共的基函数的数目为$M$个，如果定义$\phi_{0}(x)=1$的话，那个上面的式子就可以简单的表示为：

$$y(x,w)=\sum_{j=0}^{M}\omega_{j}\phi_{j}(x)=w^T\phi(x)$$

$$w=(w_{0},w_{1},w_{2},...,w_{M})$$

$$\phi=(\phi_{0},\phi_{1},\phi_{2},...,\phi_{M})$$
则线性模型的概率表示如下：

$$p(t|x,w,\beta)=N(t|y(x,w),\beta^{-1}I)$$

假设参数$w$满足高斯分布，这是一个先验分布：

$$p(w)=N(w|0,\alpha^{-1}I)$$

一般来说，我们称$p(w)$为共轭先验(conjugate prior)。这里$t$是$x$对应的目标输出，$\beta^{-1}$和$\alpha^{-1}$分别对应于样本集合和$w$的高斯分布的方差，$w$是参数，

那么，线性模型的对数后验概率函数：

$$ln p(\theta|D)=ln p(w|T)=-\frac{\beta}{2}\sum_{n=1}^N\{y(x_{n},w)-t_{n}\}^2+\frac{\alpha}{2}w^Tw+const$$

这里$T$是数据样本的目标值向量，$T=\{t_1,t_2,...,t_n\}$，const是和参数$w$无关的量。关于这个式子的具体推导，可参见前面的文章《多项式回归》。

贝叶斯线性回归的学习过程

根据前面关于贝叶斯估计的增量学习可以很容易得到下面这个式子，这个就是贝叶斯学习过程：在前一个训练集合$D^{n-1}$的后验概率$p(\theta |D^{n-1})$上，乘以新的测试样本点$x_n$的似然估计，得到新的集合$D^n$的后验概率$p(\theta|D^n)$，这样，相当于$p(\theta |D^{n-1})$成为了$p(\theta|D^n)$的先验概率分布：

$$p(\theta|D^n) \propto p(x_n|\theta)p(\theta |D^{n-1})$$

有了上面的基础知识，这里就着重的讲下面这幅图，这个图是从RMPL第155页截取下来的，这幅图清晰的描述了贝叶斯线性回归的学习过程，下面结合这幅图，详细的说明一下贝叶斯学习过程。

首先，说一下这里的模型：

$$y(x,w)=w_0+w_1x$$

第一行：
第一行是初始状态，此时只有关于$w$的先验信息，即：$p(\theta|D^0)=p(\theta)=N(w|0,\alpha^{-1}I)$。先看中间这幅图，这幅图是关于$w$的先验分布，由于我们假设$w$初始为高斯分布$N(w|0,\alpha^{-1}I)$，故其图形是以$(0,0)$为中心的圆组成的。由于此时还没有样本点进入，所以没有关于样本的似然估计，故第一行中左边likelihood没有图。第一行右边data space的那幅图显示的是从第二幅图prior/posterior中随机抽取一些点$(w_0,w_1)$，并以$(w_0,w_1)$为参数，画出来的直线，此时这些直线是随机的。

第二行：
此时有了第一个样本点$x_1$，那么根据$x_1$就可以得到第二行中，关于$x_1$的似然估计，由于$y=w_0+w_1x$，似然估计的结果其实是这个式子的对偶式，即$w_1=\frac{1}{x}y-\frac{1}{x}w_0$。从第二行的最右边data space中的图中可以估计出，第一个样本点的坐标大概为:$(0.9,0.1)$，所以其第一幅图中，似然估计的中心线的方程为:

$$w_1=\frac{1}{9}-\frac{10}{9}w_0$$

近似为左边那幅图的画法。由于第二行的先验分布是第一行的后验分布，也就是第一行的中间那幅图。则，第二行的后验分布的求法就是：将第二行的第左边那幅图和第一行的中间那幅图相乘，就可以得到第二行中间那幅图。第二行最右边那幅图就是从第二行中间那幅图中随机抽取一些点$(w_0,w_1)$，并以$(w_0,w_1)$为参数，画出来的直线。

第三行之后，就可以一次类推了。

上面就是贝叶斯学习过程的完整描述。

贝叶斯回归的优缺点

优点：
1. 贝叶斯回归对数据有自适应能力，可以重复的利用实验数据，并防止过拟合
2. 贝叶斯回归可以在估计过程中引入正则项
缺点：
1. 贝叶斯回归的学习过程开销太大

贝叶斯脊回归（Bayesian Ridge Regression）

前面已经证明过了，如果贝叶斯线性回归的先验分布为

$$p(w)=N(w|0,\alpha^{-1}I)$$

那么，其最终的后验分布公式为：

$$ln p(\theta|D)=ln p(w|T)=-\frac{\beta}{2}\sum_{n=1}^N\{y(x_{n},w)-t_{n}\}^2+\frac{\alpha}{2}w^Tw+const$$

这个相当于脊回归，所以将这种特殊情况称为贝叶斯脊回归，它拥有脊回归的所有特性，具体可以参见前面的文章《脊回归》。

下面这份代码提供了贝叶斯回归的用法，以及其和最小二乘法的比较。

#!/usr/bin/python
# -*- coding: utf-8 -*-

"""
author ： duanxxnj@163.com
time : 2016-06-21-09-21

贝叶斯脊回归
这里在一个自己生成的数据集合上测试贝叶斯脊回归

贝叶斯脊回归和最小二乘法(OLS)得到的线性模型的参数是有一定的差别的
相对于最小二乘法(OLS)二样，贝叶斯脊回归得到的参数比较接近于0

贝叶斯脊回归的先验分布是参数向量的高斯分布

"""
print(__doc__)

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
import time

from sklearn.linear_model import BayesianRidge, LinearRegression

###############################################################################
# 随机函数的种子
np.random.seed(int(time.time()) % 100)
# 样本数目为100，特征数目也是100
n_samples, n_features = 100, 100
# 生成高斯分布
X = np.random.randn(n_samples, n_features)
# 首先使用alpha为4的先验分布.
alpha_ = 4.
w = np.zeros(n_features)
# 随机提取10个特征出来作为样本特征
relevant_features = np.random.randint(0, n_features, 10)
# 基于先验分布，产生特征对应的初始权值
for i in relevant_features:
    w[i] = stats.norm.rvs(loc=0, scale=1. / np.sqrt(alpha_))

# 产生alpha为50的噪声
alpha_ = 50.
noise = stats.norm.rvs(loc=0, scale=1. / np.sqrt(alpha_), size=n_samples)
# 产生目标数据
y = np.dot(X, w) + noise

###############################################################################
# 使用贝叶斯脊回归拟合数据
clf = BayesianRidge(compute_score=True)
clf.fit(X, y)

# 使用最小二乘法拟合数据
ols = LinearRegression()
ols.fit(X, y)

###############################################################################
# 作图比较两个方法的结果
plt.figure(figsize=(6, 5))
plt.title("Weights of the model")
plt.plot(clf.coef_, ‘b-‘, label="Bayesian Ridge estimate")
plt.plot(w, ‘g-‘, label="Ground truth")
plt.plot(ols.coef_, ‘r--‘, label="OLS estimate")
plt.xlabel("Features")
plt.ylabel("Values of the weights")
plt.legend(loc="best", prop=dict(size=12))

plt.show()

其运行结果为：

贝叶斯脊回归参数：
[ -7.16688614e-02   3.73638195e-02  -4.04171217e-02   7.28338457e-03
 ...
   2.60774221e-01   4.26079127e-02   1.01660304e-02   6.79853349e-02]
最小二乘法参数：
[-1.77270023 -0.38832798  0.58907738 -1.61514115  0.58202424  0.09483505
...
 -1.37056305  2.81533169  0.02429617  0.90196961]

可以很容易的看出，相对于最小二乘法(OLS)二样，贝叶斯脊回归得到的参数比较接近于0