线性回归小结（转） - 开发者知识库

[原创] 线性回归小结（转）

秋水 2018/03/02 24 线性回归/ 回归

线性回归小结
　　线性回归可以说是机器学习中最基本的问题类型了，这里就对线性回归的原理和算法做一个小结。

1. 线性回归的模型函数和损失函数

线性回归遇到的问题一般是这样的。我们有m个样本，每个样本对应于n维特征和一个结果输出，如下：

我们的问题是，对于一个新的(), 他所对应的是多少呢？ 如果这个问题里面的y是连续的，则是一个回归问题，否则是一个分类问题。

对于n维特征的样本数据，如果我们决定使用线性回归，那么对应的模型是这样的：

, 其中 (i = 0,1,2... n)为模型参数， (i = 0,1,2... n)为每个样本的n个特征值。这个表示可以简化，我们增加一个特征 ，这样。

进一步用矩阵形式表达更加简洁如下：

其中， 假设函数为的向量,θ为nx1的向量，里面有n个代数法的模型参数。为mxn维的矩阵。代表样本的个数，代表样本的特征数。

得到了模型，我们需要求出需要的损失函数，一般线性回归我们用均方误差作为损失函数。损失函数的代数法表示如下：

进一步用矩阵形式表达损失函数：

由于矩阵法表达比较的简洁，后面我们将统一采用矩阵方式表达模型函数和损失函数。

2. 线性回归的算法

对于线性回归的损失函数，我们常用的有两种方法来求损失函数最小化时候的θ参数：一种是梯度下降法，一种是最小二乘法。由于已经在其它篇中单独介绍了梯度下降法和最小二乘法，可以点链接到对应的文章链接去阅读。

如果采用梯度下降法，则θ的迭代公式是这样的：

通过若干次迭代后，我们可以得到最终的θ的结果

如果采用最小二乘法，则θ的结果公式如下：

当然线性回归，还有其他的常用算法，比如牛顿法和拟牛顿法，这里不详细描述。

3. 线性回归的推广：多项式回归

回到我们开始的线性模型，, 如果这里不仅仅是x的一次方，比如增加二次方，那么模型就变成了多项式回归。这里写一个只有两个特征的p次方多项式回归的模型：

我们令 ,这样我们就得到了下式：

可以发现，我们又重新回到了线性回归，这是一个五元线性回归，可以用线性回归的方法来完成算法。对于每个二元样本特征(),我们得到一个五元样本特征()，通过这个改进的五元样本特征，我们重新把不是线性回归的函数变回线性回归。

4. 线性回归的推广：广义线性回归

在上一节的线性回归的推广中，我们对样本特征端做了推广，这里我们对于特征y做推广。比如我们的输出不满足和X的线性关系，但是 和满足线性关系，模型函数如下：

这样对与每个样本的输入y，我们用 lny去对应， 从而仍然可以用线性回归的算法去处理这个问题。我们把 Iny一般化，假设这个函数是单调可微函数,则一般化的广义线性回归形式是：

或者

这个函数我们通常称为联系函数。

5. 线性回归的正则化

为了防止模型的过拟合，我们在建立线性模型的时候经常需要加入正则化项。一般有L1正则化和L2正则化。

线性回归的L1正则化通常称为Lasso回归，它和一般线性回归的区别是在损失函数上增加了一个L1正则化的项，L1正则化的项有一个常数系数来调节损失函数的均方差项和正则化项的权重，具体Lasso回归的损失函数表达式如下：

其中n为样本个数，为常数系数，需要进行调优。为L1范数。

Lasso回归可以使得一些特征的系数变小，甚至还是一些绝对值较小的系数直接变为0。增强模型的泛化能力。

Lasso回归的求解办法一般有坐标轴下降法（coordinate descent）和最小角回归法（ Least Angle Regression），由于它们比较复杂，在我的这篇文章单独讲述： 线程回归的正则化-Lasso回归小结

线性回归的L2正则化通常称为Ridge回归，它和一般线性回归的区别是在损失函数上增加了一个L2正则化的项，和Lasso回归的区别是Ridge回归的正则化项是L2范数，而Lasso回归的正则化项是L1范数。具体Ridge回归的损失函数表达式如下：

其中为常数系数，需要进行调优。为L2范数。

Ridge回归在不抛弃任何一个特征的情况下，缩小了回归系数，使得模型相对而言比较的稳定，但和Lasso回归比，这会使得模型的特征留的特别多，模型解释性差。

Ridge回归的求解比较简单，一般用最小二乘法。这里给出用最小二乘法的矩阵推导形式，和普通线性回归类似。

令的导数为0，得到下式：

整理即可得到最后的θ的结果：

其中为单位矩阵。

除了上面这两种常见的线性回归正则化，还有一些其他的线性回归正则化算法，区别主要就在于正则化项的不同，和损失函数的优化方式不同，这里就不累述了。

Lasso回归算法： 坐标轴下降法与最小角回归法
　　前面的文章对线性回归做了一个小结，文章在这： 线性回归原理小结。里面对线程回归的正则化也做了一个初步的介绍。提到了线程回归的L2正则化-Ridge回归，以及线程回归的L1正则化-Lasso回归。但是对于Lasso回归的解法没有提及，本文是对该文的补充和扩展。以下都用矩阵法表示，如果对于矩阵分析不熟悉，推荐学习张贤达的《矩阵分析与应用》。

1. 回顾线性回归

首先我们简要回归下线性回归的一般形式：

需要极小化的损失函数是：

如果用梯度下降法求解，则每一轮θ迭代的表达式是：

其中为步长。

如果用最小二乘法，则θ的结果是：

2. 回顾Ridge回归

由于直接套用线性回归可能产生过拟合，我们需要加入正则化项，如果加入的是L2正则化项，就是Ridge回归，有时也翻译为脊回归。它和一般线性回归的区别是在损失函数上增加了一个L2正则化的项，和一个调节线性回归项和正则化项权重的系数。损失函数表达式如下：

其中为常数系数，需要进行调优。为L2范数。

Ridge回归的解法和一般线性回归大同小异。如果采用梯度下降法，则每一轮θ迭代的表达式是：

其中β为步长。

如果用最小二乘法，则θ的结果是：

其中为单位矩阵。

Ridge回归在不抛弃任何一个变量的情况下，缩小了回归系数，使得模型相对而言比较的稳定，但这会使得模型的变量特别多，模型解释性差。有没有折中一点的办法呢？即又可以防止过拟合，同时克服Ridge回归模型变量多的缺点呢？有，这就是下面说的Lasso回归。

3. 初识Lasso回归

Lasso回归有时也叫做线性回归的L1正则化，和Ridge回归的主要区别就是在正则化项，Ridge回归用的是L2正则化，而Lasso回归用的是L1正则化。Lasso回归的损失函数表达式如下：

其中n为样本个数，为常数系数，需要进行调优。为L1范数。

Lasso回归使得一些系数变小，甚至还是一些绝对值较小的系数直接变为0，因此特别适用于参数数目缩减与参数的选择，因而用来估计稀疏参数的线性模型。

但是Lasso回归有一个很大的问题，导致我们需要把它单独拎出来讲，就是它的损失函数不是连续可导的，由于L1范数用的是绝对值之和，导致损失函数有不可导的点。也就是说，我们的最小二乘法，梯度下降法，牛顿法与拟牛顿法对它统统失效了。那我们怎么才能求有这个L1范数的损失函数极小值呢？

OK，本章主角，两种全新的求极值解法坐标轴下降法（coordinate descent）和最小角回归法（ Least Angle Regression， LARS）该隆重出场了。

4. 用坐标轴下降法求解Lasso回归

坐标轴下降法顾名思义，是沿着坐标轴的方向去下降，这和梯度下降不同。梯度下降是沿着梯度的负方向下降。不过梯度下降和坐标轴下降的共性就都是迭代法，通过启发式的方式一步步迭代求解函数的最小值。

坐标轴下降法的数学依据主要是这个结论（此处不做证明）：一个可微的凸函数, 其中θ是nx1的向量，即有n个维度。如果在某一点，使得J(θ)在每一个坐标轴(i = 1,2,...n)上都是最小值，那么就是一个全局的最小值。

于是我们的优化目标就是在θ的n个坐标轴上(或者说向量的方向上)对损失函数做迭代的下降，当所有的坐标轴上的都达到收敛时，我们的损失函数最小，此时的θ即为我们要求的结果。

下面我们看看具体的算法过程：

首先，我们把θ向量随机取一个初值。记为 ，上面的括号里面的数字代表我们迭代的轮数，当前初始轮数为0.

对于第k轮的迭代。我们从开始，到为止，依次求。的表达式如下：

也就是说是使最小化时候的的值。此时只有是变量，其余均为常量，因此最小值容易通过求导求得。

如果上面这个式子不好理解，我们具体一点，在第k轮，θ向量的n个维度的迭代式如下：

...

检查向量和向量在各个维度上的变化情况，如果在所有维度上变化都足够小，那么即为最终结果，否则转入2，继续第k+1轮的迭代。

以上就是坐标轴下降法的求极值过程，可以和梯度下降做一个比较：

a) 坐标轴下降法在每次迭代中在当前点处沿一个坐标方向进行一维搜索 ，固定其他的坐标方向，找到一个函数的局部极小值。而梯度下降总是沿着梯度的负方向求函数的局部最小值。

b) 坐标轴下降优化方法是一种非梯度优化算法。在整个过程中依次循环使用不同的坐标方向进行迭代，一个周期的一维搜索迭代过程相当于一个梯度下降的迭代。

c) 梯度下降是利用目标函数的导数来确定搜索方向的，该梯度方向可能不与任何坐标轴平行。而坐标轴下降法法是利用当前坐标方向进行搜索，不需要求目标函数的导数，只按照某一坐标方向进行搜索最小值。

d) 两者都是迭代方法，且每一轮迭代，都需要O(mn)的计算量(m为样本数，n为系数向量的维度)

5. 用最小角回归法求解Lasso回归

第四节介绍了坐标轴下降法求解Lasso回归的方法，此处再介绍另一种常用方法， 最小角回归法(Least Angle Regression， LARS)。

在介绍最小角回归前，我们先看看两个预备算法，好吧，这个算法真没有那么好讲。

5.1 前向选择（Forward Selection）算法

第一个预备算法是前向选择（Forward Selection）算法。

前向选择算法的原理是是一种典型的贪心算法。要解决的问题是对于:

这样的线性关系，如何求解系数向量θ的问题。其中为 mx1的向量，为mxn的矩阵，θ为nx1的向量。m为样本数量，n为特征维度。

把 矩阵看做n个mx1的向量(i=1,2,...n)，在的变量(i =1,2,...m)中，选择和目标最为接近(余弦距离最大)的一个变量，用来逼近,得到下式：

其中：

即： 是 在 上的投影。那么，可以定义残差(residual):   。由于是投影，所以很容易知道 和是正交的。再以为新的因变量，去掉后，剩下的自变量的集合,i=1,2,3...k−1,k+1,...n}为新的自变量集合，重复刚才投影和残差的操作，直到残差为0，或者所有的自变量都用完了，才停止算法。

img

当只有2维时，例子如上图，和Y最接近的是，首先在上面投影，残差如上图长虚线。此时模拟了，模拟了θ(仅仅模拟了一个维度)。接着发现最接近的是，此时用残差接着在投影，残差如图中短虚线。由于没有其他自变量了，此时模拟了,对应的模拟了两个维度的θ即为最终结果，此处θ计算设计较多矩阵运算，这里不讨论。

此算法对每个变量只需要执行一次操作，效率高，速度快。但也容易看出，当自变量不是正交的时候，由于每次都是在做投影，所有算法只能给出一个局部近似解。因此，这个简单的算法太粗糙，还不能直接用于我们的Lasso回归。

5.2 前向梯度（Forward Stagewise）算法

第二个预备算法是前向梯度（Forward Stagewise）算法。

前向梯度算法和前向选择算法有类似的地方，也是在的变量(i =1,2,...n)中，选择和目标最为接近(余弦距离最大)的一个变量，用来逼近，但是前向梯度算法不是粗暴的用投影，而是每次在最为接近的自变量的方向移动一小步，然后再看残差和哪个(i =1,2,...n)最为接近。此时我们也不会把 去除，因为我们只是前进了一小步，有可能下面最接近的自变量还是。如此进行下去，直到残差减小到足够小，算法停止。

当只有2维时，例子如上图，和最接近的是，首先在上面走一小段距离，此处ε为一个较小的常量，发现此时的残差还是和最接近。那么接着沿走，一直走到发现残差不是和最接近，而是和最接近，此时残差如上图长虚线。接着沿着走一小步，发现残差此时又和最接近，那么开始沿着走，走完一步后发现残差为0，那么算法停止。此时由刚才所有的所有步相加而模拟，对应的算出的系数θ即为最终结果。此处θ计算涉及较多矩阵运算，这里不讨论。

当算法在ε很小的时候，可以很精确的给出最优解，当然，其计算的迭代次数也是大大的增加。和前向选择算法相比，前向梯度算法更加精确，但是更加复杂。

有没有折中的办法可以综合前向梯度算法和前向选择算法的优点，做一个折中呢？有！这就是终于要出场的最小角回归法。

5.3  最小角回归(Least Angle Regression， LARS)算法

好吧，最小角回归(Least Angle Regression， LARS)算法终于出场了。最小角回归法对前向梯度算法和前向选择算法做了折中，保留了前向梯度算法一定程度的精确性，同时简化了前向梯度算法一步步迭代的过程。具体算法是这样的：

首先，还是找到与因变量最接近或者相关度最高的自变量，使用类似于前向梯度算法中的残差计算方法，得到新的目标，此时不用和前向梯度算法一样小步小步的走。而是直接向前走直到出现一个，使得和的相关度和与的相关度是一样的，此时残差就在和的角平分线方向上，此时我们开始沿着这个残差角平分线走，直到出现第三个特征和的相关度足够大的时候，即到当前残差的相关度和，与的一样。将其也加入到Y的逼近特征集合中，并用的逼近特征集合的共同角分线，作为新的逼近方向。以此循环，直到足够的小，或者说所有的变量都已经取完了，算法停止。此时对应的系数即为最终结果。

img

当θ只有2维时，例子如上图，和Y最接近的是，首先在上面走一段距离，一直到残差在和的角平分线上，此时沿着角平分线走，直到残差最够小时停止，此时对应的系数β即为最终结果。此处θ计算设计较多矩阵运算，这里不讨论。

最小角回归法是一个适用于高维数据的回归算法，其主要的优点有：

1）特别适合于特征维度n 远高于样本数m的情况。

2）算法的最坏计算复杂度和最小二乘法类似，但是其计算速度几乎和前向选择算法一样

3）可以产生分段线性结果的完整路径，这在模型的交叉验证中极为有用

主要的缺点是：

由于LARS的迭代方向是根据目标的残差而定，所以该算法对样本的噪声极为敏感。

6. 总结

Lasso回归是在ridge回归的基础上发展起来的，如果模型的特征非常多，需要压缩，那么Lasso回归是很好的选择。一般的情况下，普通的线性回归模型就够了。

另外，本文对最小角回归法怎么求具体的θ参数值没有提及，仅仅涉及了原理，如果对具体的算计推导有兴趣，可以参考Bradley Efron的论文《Least Angle Regression》，网上很容易找到。

MCMC(一)蒙特卡罗方法
　　作为一种随机采样方法，马尔科夫链蒙特卡罗（Markov Chain Monte Carlo，以下简称MCMC）在机器学习,深度学习以及自然语言处理等领域都有广泛的应用，是很多复杂算法求解的基础。比如我们前面讲到的分解机(Factorization Machines)推荐算法，还有前面讲到的受限玻尔兹曼机（RBM）原理总结，都用到了MCMC来做一些复杂运算的近似求解。下面我们就对MCMC的原理做一个总结。

1. MCMC概述

从名字我们可以看出，MCMC由两个MC组成，即蒙特卡罗方法（Monte Carlo Simulation，简称MC）和马尔科夫链（Markov Chain ，也简称MC）。要弄懂MCMC的原理我们首先得搞清楚蒙特卡罗方法和马尔科夫链的原理。我们将用三篇来完整学习MCMC。在本篇，我们关注于蒙特卡罗方法。

2. 蒙特卡罗方法引入

蒙特卡罗原来是一个赌场的名称，用它作为名字大概是因为蒙特卡罗方法是一种随机模拟的方法，这很像赌博场里面的扔骰子的过程。最早的蒙特卡罗方法都是为了求解一些不太好求解的求和或者积分问题。比如积分：

如果我们很难求解出的原函数，那么这个积分比较难求解。当然我们可以通过蒙特卡罗方法来模拟求解近似值。如何模拟呢？假设我们函数图像如下图:

img

则一个简单的近似求解方法是在[a,b]之间随机的采样一个点。比如,然后用代表在[a,b]区间上所有的的值。那么上面的定积分的近似求解为:

当然，用一个值代表[a,b]区间上所有的的值，这个假设太粗糙。那么我们可以采样[a,b]区间的n个值：,用它们的均值来代表[a,b]区间上所有的的值。这样我们上面的定积分的近似求解为:

虽然上面的方法可以一定程度上求解出近似的解，但是它隐含了一个假定，即在[a,b]之间是均匀分布的，而绝大部分情况，在[a,b]之间不是均匀分布的。如果我们用上面的方法，则模拟求出的结果很可能和真实值相差甚远。

怎么解决这个问题呢？ 如果我们可以得到在[a,b]的概率分布函数p(x)，那么我们的定积分求和可以这样进行：

上式最右边的这个形式就是蒙特卡罗方法的一般形式。当然这里是连续函数形式的蒙特卡罗方法，但是在离散时一样成立。

可以看出，最上面我们假设xx在[a,b]之间是均匀分布的时候，p(xi)=1/(b−a)p(xi)=1/(b−a)，带入我们有概率分布的蒙特卡罗积分的上式，可以得到：

也就是说，我们最上面的均匀分布也可以作为一般概率分布函数在均匀分布时候的特例。那么我们现在的问题转到了如何求出的分布的若干和样本上来。

3. 概率分布采样

上一节我们讲到蒙特卡罗方法的关键是得到的概率分布。如果求出了的概率分布，我们可以基于概率分布去采样基于这个概率分布的n个的样本集，代入蒙特卡罗求和的式子即可求解。但是还有一个关键的问题需要解决，即如何基于概率分布去采样基于这个概率分布的n个的样本集。

对于常见的均匀分布uniform(0,1)是非常容易采样样本的，一般通过线性同余发生器可以很方便的生成(0,1)之间的伪随机数样本。而其他常见的概率分布，无论是离散的分布还是连续的分布，它们的样本都可以通过uniform(0,1)的样本转换而得。比如二维正态分布的样本可以通过通过独立采样得到的uniform(0,1)样本对通过如下的式子转换而得：

其他一些常见的连续分布，比如t分布，F分布，Beta分布，Gamma分布等，都可以通过类似的方式从uniform(0,1)得到的采样样本转化得到。在python的numpy，scikit-learn等类库中，都有生成这些常用分布样本的函数可以使用。

不过很多时候，我们的xx的概率分布不是常见的分布，这意味着我们没法方便的得到这些非常见的概率分布的样本集。那这个问题怎么解决呢？

4. 接受-拒绝采样

对于概率分布不是常见的分布，一个可行的办法是采用接受-拒绝采样来得到该分布的样本。既然  太复杂在程序中没法直接采样，那么我设定一个程序可采样的分布  比如高斯分布，然后按照一定的方法拒绝某些样本，以达到接近 分布的目的，其中叫做 proposal distribution。

img

具体采用过程如下，设定一个方便采样的常用概率分布函数 ，以及一个常量 k，使得 总在  的下方。如上图。

首先，采样得到的一个样本，采样方法如第三节。然后，从均匀分布（)中采样得到一个值u。如果u落在了上图中的灰色区域，则拒绝这次抽样，否则接受这个样本。重复以上过程得到n个接受的样本,则最后的蒙特卡罗方法求解结果为：

整个过程中，我们通过一系列的接受拒绝决策来达到用模拟概率分布的目的。

5. 蒙特卡罗方法小结

使用接受-拒绝采样，我们可以解决一些概率分布不是常见的分布的时候，得到其采样集并用蒙特卡罗方法求和的目的。但是接受-拒绝采样也只能部分满足我们的需求，在很多时候我们还是很难得到我们的概率分布的样本集。比如：

1）对于一些二维分布，有时候我们只能得到条件分布和,却很难得到二维分布一般形式，这时我们无法用接受-拒绝采样得到其样本集。

2）对于一些高维的复杂非常见分布，我们要找到一个合适的和非常困难。

从上面可以看出，要想将蒙特卡罗方法作为一个通用的采样模拟求和的方法，必须解决如何方便得到各种复杂概率分布的对应的采样样本集的问题。而我们下一篇要讲到的马尔科夫链就是帮助找到这些复杂概率分布的对应的采样样本集的白衣骑士。下一篇我们来总结马尔科夫链的原理。

MCMC(二)马尔科夫链
　　在MCMC(一)蒙特卡罗方法中，我们讲到了如何用蒙特卡罗方法来随机模拟求解一些复杂的连续积分或者离散求和的方法，但是这个方法需要得到对应的概率分布的样本集，而想得到这样的样本集很困难。因此我们需要本篇讲到的马尔科夫链来帮忙。

1. 马尔科夫链概述

马尔科夫链定义本身比较简单，它假设某一时刻状态转移的概率只依赖于它的前一个状态。举个形象的比喻，假如每天的天气是一个状态的话，那个今天是不是晴天只依赖于昨天的天气，而和前天的天气没有任何关系。当然这么说可能有些武断，但是这样做可以大大简化模型的复杂度，因此马尔科夫链在很多时间序列模型中得到广泛的应用，比如循环神经网络RNN，隐式马尔科夫模型HMM等，当然MCMC也需要它。

如果用精确的数学定义来描述，则假设我们的序列状态是，那么我们的在时刻的状态的条件概率仅仅依赖于时刻，即：

既然某一时刻状态转移的概率只依赖于它的前一个状态，那么我们只要能求出系统中任意两个状态之间的转换概率，这个马尔科夫链的模型就定了。我们来看看下图这个马尔科夫链模型的具体的例子(来源于维基百科)。

img

这个马尔科夫链是表示股市模型的，共有三种状态：牛市（Bull market）, 熊市（Bear market）和横盘（Stagnant market）。每一个状态都以一定的概率转化到下一个状态。比如，牛市以0.025的概率转化到横盘的状态。这个状态概率转化图可以以矩阵的形式表示。如果我们定义矩阵阵某一位置的值为,即从状态转化到状态的概率，并定义牛市为状态0， 熊市为状态1, 横盘为状态2. 这样我们得到了马尔科夫链模型的状态转移矩阵为：

讲了这么多，那么马尔科夫链模型的状态转移矩阵和我们蒙特卡罗方法需要的概率分布样本集有什么关系呢？这需要从马尔科夫链模型的状态转移矩阵的性质讲起。

2. 马尔科夫链模型状态转移矩阵的性质

得到了马尔科夫链模型的状态转移矩阵，我们来看看马尔科夫链模型的状态转移矩阵的性质。

仍然以上面的这个状态转移矩阵为例。假设我们当前股市的概率分布为：[0.3,0.4,0.3],即30%概率的牛市，40%概率的熊盘与30%的横盘。然后这个状态作为序列概率分布的初始状态t0t0，将其带入这个状态转移矩阵计算t1,t2,t3...t1,t2,t3...的状态。代码如下：

import numpy as np
matrix = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]], dtype=float)
vector1 = np.matrix([[0.3,0.4,0.3]], dtype=float)
for i in range(100):
    vector1 = vector1*matrix
    print "Current round:" , i+1
    print vector1
部分输出结果如下：

Current round: 1
[[ 0.405   0.4175  0.1775]]
Current round: 2
[[ 0.4715   0.40875  0.11975]]
Current round: 3
[[ 0.5156  0.3923  0.0921]]
Current round: 4
[[ 0.54591   0.375535  0.078555]]
。。。。。。
Current round: 58
[[ 0.62499999  0.31250001  0.0625    ]]
Current round: 59
[[ 0.62499999  0.3125      0.0625    ]]
Current round: 60
[[ 0.625   0.3125  0.0625]]
。。。。。。
Current round: 99
[[ 0.625   0.3125  0.0625]]
Current round: 100
[[ 0.625   0.3125  0.0625]]
　　可以发现，从第60轮开始，我们的状态概率分布就不变了，一直保持在[0.625   0.3125  0.0625]，即62.5%的牛市，31.25%的熊市与6.25%的横盘。那么这个是巧合吗？

我们现在换一个初始概率分布试一试，现在我们用[0.7,0.1,0.2]作为初始概率分布，然后这个状态作为序列概率分布的初始状态，将其带入这个状态转移矩阵计算t1,t2,t3...的状态。代码如下：

matrix = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]], dtype=float)
vector1 = np.matrix([[0.7,0.1,0.2]], dtype=float)
for i in range(100):
    vector1 = vector1*matrix
    print "Current round:" , i+1
    print vector1
部分输出结果如下：

Current round: 1
[[ 0.695   0.1825  0.1225]]
Current round: 2
[[ 0.6835   0.22875  0.08775]]
Current round: 3
[[ 0.6714  0.2562  0.0724]]
Current round: 4
[[ 0.66079   0.273415  0.065795]]
。。。。。。。
Current round: 55
[[ 0.62500001  0.31249999  0.0625    ]]
Current round: 56
[[ 0.62500001  0.31249999  0.0625    ]]
Current round: 57
[[ 0.625   0.3125  0.0625]]
。。。。。。。
Current round: 99
[[ 0.625   0.3125  0.0625]]
Current round: 100
[[ 0.625   0.3125  0.0625]]
　　可以看出，尽管这次我们采用了不同初始概率分布，最终状态的概率分布趋于同一个稳定的概率分布[0.625   0.3125  0.0625]， 也就是说我们的马尔科夫链模型的状态转移矩阵收敛到的稳定概率分布与我们的初始状态概率分布无关。这是一个非常好的性质，也就是说，如果我们得到了这个稳定概率分布对应的马尔科夫链模型的状态转移矩阵，则我们可以用任意的概率分布样本开始，带入马尔科夫链模型的状态转移矩阵，这样经过一些序列的转换，最终就可以得到符合对应稳定概率分布的样本。

这个性质不光对我们上面的状态转移矩阵有效，对于绝大多数的其他的马尔科夫链模型的状态转移矩阵也有效。同时不光是离散状态，连续状态时也成立。

同时，对于一个确定的状态转移矩阵，它的n次幂在当n大于一定的值的时候也可以发现是确定的，我们还是以上面的例子为例，计算代码如下：

matrix = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]], dtype=float)
for i in range(10):
    matrix = matrix*matrix
    print "Current round:" , i+1
    print matrix
输出结果如下：

Current round: 1
[[ 0.8275   0.13375  0.03875]
 [ 0.2675   0.66375  0.06875]
 [ 0.3875   0.34375  0.26875]]
Current round: 2
[[ 0.73555   0.212775  0.051675]
 [ 0.42555   0.499975  0.074475]
 [ 0.51675   0.372375  0.110875]]
。。。。。。
Current round: 5
[[ 0.62502532  0.31247685  0.06249783]
 [ 0.6249537   0.31254233  0.06250397]
 [ 0.62497828  0.31251986  0.06250186]]
Current round: 6
[[ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]]
Current round: 7
[[ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]]
。。。。。。
Current round: 9
[[ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]]
Current round: 10
[[ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]]　　
　　我们可以发现，在以后，的值稳定不再变化，而且每一行都为[0.625   0.3125  0.0625]，这和我们前面的稳定分布是一致的。这个性质同样不光是离散状态，连续状态时也成立。

好了，现在我们可以用数学语言总结下马尔科夫链的收敛性质了：

如果一个非周期的马尔科夫链有状态转移矩阵, 并且它的任何两个状态是连通的，那么与无关，我们有：

1）

3）

4) 是方程的唯一非负解，其中：

上面的性质中需要解释的有：

1）非周期的马尔科夫链：这个主要是指马尔科夫链的状态转化不是循环的，如果是循环的则永远不会收敛。幸运的是我们遇到的马尔科夫链一般都是非周期性的。用数学方式表述则是：对于任意某一状态i，d为集合{}的最大公约数，如果 d=1 ，则该状态为非周期的

2）任何两个状态是连通的：这个指的是从任意一个状态可以通过有限步到达其他的任意一个状态，不会出现条件概率一直为0导致不可达的情况。

3）马尔科夫链的状态数可以是有限的，也可以是无限的。因此可以用于连续概率分布和离散概率分布。

4）通常称为马尔科夫链的平稳分布。

3. 基于马尔科夫链采样

如果我们得到了某个平稳分布所对应的马尔科夫链状态转移矩阵，我们就很容易采用出这个平稳分布的样本集。

假设我们任意初始的概率分布是, 经过第一轮马尔科夫链状态转移后的概率分布是，。。。第i轮的概率分布是。假设经过n轮后马尔科夫链收敛到我们的平稳分布，即：

对于每个分布，我们有：

现在我们可以开始采样了，首先，基于初始任意简单概率分布比如高斯分布采样得到状态值，基于条件概率分布采样状态值，一直进行下去，当状态转移进行到一定的次数时，比如到次时，我们认为此时的采样集()即是符合我们的平稳分布的对应样本集，可以用来做蒙特卡罗模拟求和了。

总结下基于马尔科夫链的采样过程：

1）输入马尔科夫链状态转移矩阵，设定状态转移次数阈值，需要的样本个数

2）从任意简单概率分布采样得到初始状态值

3）for t=0 to : 从条件概率分布中采样得到样本

样本集(即为我们需要的平稳分布对应的样本集。

4. 马尔科夫链采样小结

如果假定我们可以得到我们需要采样样本的平稳分布所对应的马尔科夫链状态转移矩阵，那么我们就可以用马尔科夫链采样得到我们需要的样本集，进而进行蒙特卡罗模拟。但是一个重要的问题是，随意给定一个平稳分布ππ,如何得到它所对应的马尔科夫链状态转移矩阵PP呢？这是个大问题。我们绕了一圈似乎还是没有解决任意概率分布采样样本集的问题。

幸运的是，MCMC采样通过迂回的方式解决了上面这个大问题，我们在下一篇来讨论MCMC的采样，以及它的使用改进版采样: M-H采样和Gibbs采样。

MCMC(三)MCMC采样和M-H采样
　　在MCMC(二)马尔科夫链中我们讲到给定一个概率平稳分布π, 很难直接找到对应的马尔科夫链状态转移矩阵P。而只要解决这个问题，我们就可以找到一种通用的概率分布采样方法，进而用于蒙特卡罗模拟。本篇我们就讨论解决这个问题的办法：MCMC采样和它的易用版M-H采样。

1. 马尔科夫链的细致平稳条件

在解决从平稳分布π, 找到对应的马尔科夫链状态转移矩阵P之前，我们还需要先看看马尔科夫链的细致平稳条件。定义如下：

如果非周期马尔科夫链的状态转移矩阵和概率分布对于所有的满足：

则称概率分布是状态转移矩阵的平稳分布。

证明很简单,由细致平稳条件有：

将上式用矩阵表示即为：

即满足马尔可夫链的收敛性质。也就是说，只要我们找到了可以使概率分布满足细致平稳分布的矩阵即可。这给了我们寻找从平稳分布, 找到对应的马尔科夫链状态转移矩阵的新思路。

不过不幸的是，仅仅从细致平稳条件还是很难找到合适的矩阵。比如我们的目标平稳分布是,随机找一个马尔科夫链状态转移矩阵,它是很难满足细致平稳条件的，即：

那么如何使这个等式满足呢？下面我们来看MCMC采样如何解决这个问题。

2. MCMC采样

由于一般情况下，目标平稳分布π(x)和某一个马尔科夫链状态转移矩阵不满足细致平稳条件，即

我们可以对上式做一个改造，使细致平稳条件成立。方法是引入一个,使上式可以取等号，即：

问题是什么样的可以使等式成立呢？其实很简单，只要满足下两式即可：

这样，我们就得到了我们的分布对应的马尔科夫链状态转移矩阵，满足：

也就是说，我们的目标矩阵可以通过任意一个马尔科夫链状态转移矩阵乘以得到。我们又一般称之为接受率。取值在[0,1]之间，可以理解为一个概率值。即目标矩阵可以通过任意一个马尔科夫链状态转移矩阵以一定的接受率获得。这个很像我们在[MCMC(一)蒙特卡罗方法]第4节讲到的接受-拒绝采样，那里是以一个常用分布通过一定的接受-拒绝概率得到一个非常见分布，这里是以一个常见的马尔科夫链状态转移矩阵通过一定的接受-拒绝概率得到目标转移矩阵PP,两者的解决问题思路是类似的。

好了，现在我们来总结下MCMC的采样过程。

1）输入我们任意选定的马尔科夫链状态转移矩阵，平稳分布，设定状态转移次数阈值，需要的样本个数

2）从任意简单概率分布采样得到初始状态值

3）for t=0 to :

a) 从条件概率分布中采样得到样本

b) 从均匀分布采样u∼uniform[0,1]

c) 如果, 则接受转移，即

d) 否则不接受转移，

样本集()即为我们需要的平稳分布对应的样本集。

上面这个过程基本上就是MCMC采样的完整采样理论了，但是这个采样算法还是比较难在实际中应用，为什么呢？问题在上面第三步的c步骤，接受率这儿。由于可能非常的小，比如0.1，导致我们大部分的采样值都被拒绝转移，采样效率很低。有可能我们采样了上百万次马尔可夫链还没有收敛，也就是上面这个要非常非常的大，这让人难以接受，怎么办呢？这时就轮到我们的M-H采样出场了。

3. M-H采样

M-H采样是Metropolis-Hastings采样的简称，这个算法首先由Metropolis提出，被Hastings改进，因此被称之为Metropolis-Hastings采样或M-H采样

M-H采样解决了我们上一节MCMC采样接受率过低的问题。

我们回到MCMC采样的细致平稳条件：

我们采样效率低的原因是太小了，比如为0.1，而为0.2。即：

这时我们可以看到，如果两边同时扩大五倍，接受率提高到了0.5，但是细致平稳条件却仍然是满足的，即：

这样我们的接受率可以做如下改进，即：

通过这个微小的改造，我们就得到了可以在实际应用中使用的M-H采样算法过程如下：

1）输入我们任意选定的马尔科夫链状态转移矩阵，平稳分布，设定状态转移次数阈值，需要的样本个数

2）从任意简单概率分布采样得到初始状态值

3）for  to :

a) 从条件概率分布中采样得到样本

b) 从均匀分布采样u∼uniform[0,1]

c) 如果, 则接受转移，即

d) 否则不接受转移，

样本集()即为我们需要的平稳分布对应的样本集。

很多时候，我们选择的马尔科夫链状态转移矩阵如果是对称的，即满足,这时我们的接受率可以进一步简化为：

4. M-H采样实例

为了更容易理解，这里给出一个M-H采样的实例。

例子里，我们的目标平稳分布是一个均值3，标准差2的正态分布，而选择的马尔可夫链状态转移矩阵的条件转移概率是以ii为均值,方差1的正态分布在位置jj的值。这个例子仅仅用来让大家加深对M-H采样过程的理解。毕竟一个普通的一维正态分布用不着去用M-H采样来获得样本。代码如下：

import random
import math
from scipy.stats import norm
import matplotlib.pyplot as plt
%matplotlib inline

def norm_dist_prob(theta):
    y = norm.pdf(theta, loc=3, scale=2)
    return y

T = 5000
pi = [0 for i in range(T)]
sigma = 1
t = 0
while t δ

对应的负梯度误差为：

r(yi,f(xi))={yi−f(xi)δsign(yi−f(xi))|yi−f(xi)|≤δ|yi−f(xi)|>δ

d) 分位数损失。它对应的是分位数回归的损失函数，表达式为

L(y,f(x))=∑y≥f(x)θ|y−f(x)|+∑y<f(x)(1−θ)|y−f(x)|

其中θθ为分位数，需要我们在回归前指定。对应的负梯度误差为：

r(yi,f(xi))={θθ−1yi≥f(xi)yi<f(xi)

对于Huber损失和分位数损失，主要用于健壮回归，也就是减少异常点对损失函数的影响。

6. GBDT的正则化

和Adaboost一样，我们也需要对GBDT进行正则化，防止过拟合。GBDT的正则化主要有三种方式。

第一种是和Adaboost类似的正则化项，即步长(learning rate)。定义为νν,对于前面的弱学习器的迭代

fk(x)=fk−1(x)+hk(x)f

如果我们加上了正则化项，则有

fk(x)=fk−1(x)+νhk(x)

ν的取值范围为0<ν≤1。对于同样的训练集学习效果，较小的ν意味着我们需要更多的弱学习器的迭代次数。通常我们用步长和迭代最大次数一起来决定算法的拟合效果。

第二种正则化的方式是通过子采样比例（subsample）。取值为(0,1]。注意这里的子采样和随机森林不一样，随机森林使用的是放回抽样，而这里是不放回抽样。如果取值为1，则全部样本都使用，等于没有使用子采样。如果取值小于1，则只有一部分样本会去做GBDT的决策树拟合。选择小于1的比例可以减少方差，即防止过拟合，但是会增加样本拟合的偏差，因此取值不能太低。推荐在[0.5, 0.8]之间。

使用了子采样的GBDT有时也称作随机梯度提升树(Stochastic Gradient Boosting Tree, SGBT)。由于使用了子采样，程序可以通过采样分发到不同的任务去做boosting的迭代过程，最后形成新树，从而减少弱学习器难以并行学习的弱点。

第三种是对于弱学习器即CART回归树进行正则化剪枝。在决策树原理篇里我们已经讲过，这里就不重复了。

7. GBDT小结

GBDT终于讲完了，GDBT本身并不复杂，不过要吃透的话需要对集成学习的原理，决策树原理和各种损失函树有一定的了解。由于GBDT的卓越性能，只要是研究机器学习都应该掌握这个算法，包括背后的原理和应用调参方法。目前GBDT的算法比较好的库是xgboost。当然scikit-learn也可以。

最后总结下GBDT的优缺点。

GBDT主要的优点有：

1) 可以灵活处理各种类型的数据，包括连续值和离散值。

2) 在相对少的调参时间情况下，预测的准备率也可以比较高。这个是相对SVM来说的。

3）使用一些健壮的损失函数，对异常值的鲁棒性非常强。比如 Huber损失函数和Quantile损失函数。

GBDT的主要缺点有：

1)由于弱学习器之间存在依赖关系，难以并行训练数据。不过可以通过自采样的SGBT来达到部分并行。

注意！

本站用户原创文章，版权为本站用户，转摘请注明本文地址：http://www.silva-art.net/blog/2018/03/02/ca2ef0bf7634bd1b7bff842bb9d8f58f.html。

猜您在找

线性回归原理小结线性回归原理小结 Logistic回归线性回归（转）常用的线性回归模型小结 SPSS—回归—多元线性回归（转）