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

2016年06月21日 09:50:40 Duanxx 阅读数 54254更多

分类专栏：监督学习

本文链接：https://blog.csdn.net/daunxx/article/details/51725086

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

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

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

原文地址

贝叶斯线性回归Bayesian Linear Regression

本文的研究顺序是：

极大似然估计最大后验估计贝叶斯估计贝叶斯线性回归

关于参数估计

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

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

关于上面提到的需要确定的两个东西：概率密度函数的形式和参数，至少在机器学习的教课书上，我所看到的情况都是：给了一堆数据，然后假设其概率密度函数的形式为高斯分布，或者是混合高斯分布，那么，剩下的事情就是对高斯分布的参数， $μ$

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

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

极大似然估计

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

这里先以一个分类问题来说明一般参数估计面对的数据形式。考虑一个 $M$

前面已经介绍过了，想要确定数据的概率分布，需要知道概率密度函数的形式和参数，这里首先做一个基本假设：概率分布的形式已知，比如假设每个类别的数据都满足高斯分布，那么，似然函数就可以以参数 $θ_{i}$

为了强调概率分布 $p (x | ω_{i})$

在概率论中一直有两大学派，分别是频率学派和贝叶斯学派。简单点说，频率学派认为，概率是频率的极限，比如投硬币，当实验次数足够大时，正面朝上的频率可以认为是这枚硬币正面朝上的概率，这个是频率学派。但是，如果要预测一些未发生过的事情，比如，北极的冰山在2050年完全融化的概率，由于这个事情完全没有发生过，所以无法用频率来代替概率表示，只能研究过去几十年，北极冰山融化的速率，并将其作为先验条件，来预测北极的冰山在2050年完全融化的概率，这就是概率的贝叶斯学派。上面的问题，如果用贝叶斯学派的记法的话，是： $p (x | ω_{i}, θ_{i})$

从上面的描述中可以知道，利用每一个类 $X_{i}$

设 $x_{1}, x_{2}, . . ., x_{N}$

p (X; θ) \equiv p (x 1, x 2, . . ., x N; θ) = \prod k = 1 N p (x k; θ)

注意：这里的 $p (x_{k}; θ)$

需要重申一下，想要得到上面这个公式，是做了几个基本的假设的，第一：假设 $M$

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

θ^M L = a r g max θ \prod k = 1 N p (x k; θ)

为了得到最大值， ${\hat{θ}}_{M L}$

\partial \prod N k = 1 p ( x k ; θ ) \partial θ = 0

一般我们取其对数形式：

L (θ) \equiv l n \prod k = 1 N p (x k; θ)

\partial L ( θ ) \partial θ = \sum k = 1 N \partial l n p ( x k ; θ ) \partial θ = \sum k = 1 N

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

极大似然估计有两个非常重要的性质：渐进无偏和渐进一致性，有了这两个性质，使得极大似然估计的成为了非常简单而且实用的参数估计方法。这里假设 $θ_{0}$

渐进无偏

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

lim N \to \infty E [θ^M L] = θ 0

也就是说，这里认为估计值 ${\hat{θ}}_{M L}$

渐进一致

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

limN→∞prob{∥θ^ML−θ0∥⩽ϵ}=1limN→∞prob{‖θ^ML−θ0‖⩽ϵ}=1

这个公式还可以表示为：

lim N \to \infty E ∥ θ^M L - θ 0 ∥ 2 = 0

对于一个估计器而言，一致性是非常重要的，因为存在满足无偏性，但是不满足一致性的情况，比如， ${\hat{θ}}_{M L}$

注意：以上两个性质，都是在渐进的前提下（ $N \to \infty$

最大后验估计

在最大似然估计（MAP）中，将 $θ$

从最大后验估计开始，将参数 $θ$

这里一定要注意，在最大似然估计中，参数 $θ$

根据贝叶斯理论：

p (θ | X) = p ( θ ) p ( X | θ ) p ( X )

这就是参数 $θ$

θ^M A P = a r g max θ p (θ | X) = a r g max θ p (θ) p (X |

为了得到参数 $θ$

p ( θ | X ) \partial θ = p ( θ ) p ( X | θ ） \partial θ = 0

注意：这里 $p (X | θ)$

这里需要说明，虽然从公式上来看 $M A P = M L * p (θ)$

在MAP中， $p (θ)$

贝叶斯估计

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

贝叶斯估计核心问题

为了防止标号混淆，这里定义已有的样本集合为 $D$

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

虽然 $p (x)$

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

由于参数 $θ$

现在，将贝叶斯估计核心问题 $p (x | D)$

p (x | D) = \int p (x, θ | D) d θ = \int p (x | θ, D) p (θ | D) d θ

上面式子中， $x$

p (x | D) = \int p (x | θ) p (θ | D) d θ

下面这句话一定要理解：这里 $p (x | θ)$

p (θ | D) = p ( D | θ ) p ( θ ) p ( D ) = p ( D | θ ) p ( θ ) \int

p (D | θ) = \prod k = 1 N p (x k | θ)

上面这个式子就是贝叶斯估计最核心的公式，它把类条件概率密度 $p (x | D)$

贝叶斯估计的增量学习

为了明确的表示样本集合 $D$

p (D n | θ) = p (x n | θ) p (D n - 1 | θ)

可以很容易得到：

p (θ | D n) = p ( x n | θ ) p ( D n - 1 | θ ) p ( θ ) \int p (

当没有观测样本时，定义 $p (θ | D^{0}) = p (θ)$

贝叶斯线性回归

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

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

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

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

这里面对的模型是线性回归模型，其详细的介绍可以参见前面的文章《线性回归》,线性回归模型是一组输入变量 $x$

y (x, w) = w 0 + \sum j = 1 M ω j ϕ j (x)

这里 $ϕ_{j} (x)$

y (x, w) = \sum j = 0 M ω j ϕ j (x) = w T ϕ (x)

w = (w 0, w 1, w 2, . . ., w M)

ϕ = (ϕ 0, ϕ 1, ϕ 2, . . ., ϕ M)

则线性模型的概率表示如下：

p (t | x, w, β) = N (t | y (x, w), β - 1 I)

假设参数 $w$

p (w) = N (w | 0, α - 1 I)

一般来说，我们称 $p (w)$

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

l n p (θ | D) = l n p (w | T) = - β 2 \sum n = 1 N { y ( x n , w ) -

这里 $T$

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

根据前面关于贝叶斯估计的增量学习可以很容易得到下面这个式子，这个就是贝叶斯学习过程：在前一个训练集合 $D^{n - 1}$

p (θ | D n) \propto p (x n | θ) p (θ | D n - 1)

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

image_1almvfrg71b51132bg166ndsbs9.png-442.7kB

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

y (x, w) = w 0 + w 1 x

第一行：
第一行是初始状态，此时只有关于 $w$

第二行：
此时有了第一个样本点 $x_{1}$

w 1 = 1 9 - 10 9 w 0

近似为左边那幅图的画法。由于第二行的先验分布是第一行的后验分布，也就是第一行的中间那幅图。则，第二行的后验分布的求法就是：将第二行的第左边那幅图和第一行的中间那幅图相乘，就可以得到第二行中间那幅图。第二行最右边那幅图就是从第二行中间那幅图中随机抽取一些点 $(w_{0}, w_{1})$

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

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

贝叶斯回归的优缺点

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

贝叶斯脊回归（Bayesian Ridge Regression）

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

p (w) = N (w | 0, α - 1 I)

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

l n p (θ | D) = l n p (w | T) = - β 2 \sum n = 1 N { y ( x n , w ) -

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

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

#!/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()

figure_1-1.png-67.4kB

其运行结果为：

贝叶斯脊回归参数：
[ -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

posted on 2019-09-07 21:06 曹明阅读(4103) 评论(0) 编辑收藏举报