逻辑斯谛回归及其物理含义

逻辑斯谛回归

先来看看书上逻辑斯谛回归模型的公式吧，再了解公式的情况下，我们再来谈谈它实际的物理含义。吼吼，它貌似蛮有内涵的，也是从生物学上挖过来的一条曲线吧。

定义：
设X是连续随机变量，X服从逻辑斯蒂回归分布是指X具有下列分布函数和密度函数：

$F (x) = P (X \leq x) = 1 1 + e - ( x - μ ) / γ$ $F(x) = P(X \le x) = \frac{1}{1+ e^{-(x-\mu)/\gamma}}$
$f (x) = F' (x) = e - ( x - μ ) / γ γ ( 1 + e - ( x - μ ) / γ ) 2$ $f(x) = F'(x) = \frac{e^{-(x-\mu)/\gamma}}{\gamma(1+e^{-(x-\mu)/\gamma})^2}$
式中， $\mu$ 为位置参数， $\gamma \gt 0$ 为形状参数。

附上逻辑斯蒂分布的密度函数与分布函数。如下：
alt text

一张密度函数，一张分布函数，它们有什么物理内涵么？为什么有些分类方法可以用这种逻辑斯蒂回归模型去做分类？而不是其它函数？不急，咱们留个疑问，先来看看接下来这模型具体怎么工作，起码先知道该怎么用吧。

二项逻辑斯蒂回归模型

二项逻辑斯蒂回归模型是一种分类模型，由条件概率分布 $P(Y | X)$ 表示，形式为参数化的逻辑斯蒂分布。这里，随机变量X取值为实数，随机变量Y取值为1或0。我们通过监督学习的方法来估计模型参数。

定义：
二项逻辑斯蒂回归模型是如下的条件概率分布：

$P (Y = 1 | x) = e x p ( w \cdot x + b ) 1 + e x p ( w \cdot x + b )$ $P(Y = 1 | x ) = \frac{exp(w\cdot x+ b)}{1 + exp(w\cdot x + b)}$
$P (Y = 0 | x) = 1 1 + e x p ( w \cdot x + b )$ $P(Y = 0 | x ) = \frac{1}{1 + exp(w\cdot x + b)}$
这里， $x\in R^n$ 是输入， $Y \in {0,1}$ 是输出， $w \in R^n$ 是参数， $w$ 称为权值向量， $b$ 称为偏置， $w\cdot x$ 为 $w和x$ 的内积。

二项逻辑斯蒂回归模型的定义还是很有意思的，该模型最终形式表现为一种概率模型，这是和几何分类模型最大的区别，如SVM，KNN等。所以个人认为，概率模型的物理含义要更加清晰，但同时导致的结果是应用范围将变得狭窄很多，因为此处的随机变量 $x$ 是重点考虑的对象。是任何特征都能加权充当 $w\cdot x + b$ 的一个影响因子的么？

ok,既然是概率模型，问题定当迎刃而解了，概率有了，就差模型了。模型是什么？很简单，就是每个数据样本出现的概率的乘机，为该分类模型。《统计学习方法》中也给出了该模型，如下。

定义：
逻辑斯蒂回归模型学习时，对于给定的训练数据集 $T = {(x_1,y_1),(x_2,y_2),...,(x_n,y_n)}$ ，其中， $x_i \in R^n, y_i\in {0,1}$ ，可以应用极大似然估计法估计模型参数，从而得到逻辑斯蒂回归模型。
设：

$P (Y = 1 | x) = π (x), P (Y = 0 | x) = 1 - π (x)$ $P(Y = 1| x ) = \pi(x), P(Y = 0 | x) = 1 - \pi (x)$
似然函数为：
$\prod i = 1 N [π (x i)] y i [1 - π (x i)] 1 - y i$ $\prod_{i=1}^N[\pi(x_i)]^{y_i}[1-\pi(x_i)]^{1-y_i}$
对数似然函数为：
$L (w) = \sum i = 1 N [y i log π (x i) + (1 - y i) log (1 - π (x i))] = \sum i = 1 N [y i log (π ( x i ) 1 - π ( x i )) + log (1 - π (x i))] = \sum i = 1 N [y i (w \cdot x i) - log (1 + e x p (w \cdot x i))]$ $\begin{align*} L(w) & = \sum_{i = 1}^N[y_i\log\pi(x_i)+ (1-y_i)\log(1-\pi(x_i))]\\ &= \sum_{i=1}^N[y_i\log(\frac{\pi(x_i)}{1-\pi(x_i)})+\log(1-\pi(x_i))]\\ &=\sum_{i=1}^N[y_i(w\cdot x_i) - \log (1+ exp(w\cdot x_i))] \end{align*}$

对 $L(w)$ 求极大值，得到 $w$ 的估计值。这样，问题就变成了以对数似然函数为目标函数的最优化问题。逻辑斯蒂回归学习中通常采用的方法是梯度下降法及拟牛顿法。

物理含义

这就是逻辑斯蒂回归的主要过程了，接下来我们来研究下它实际的物理意义。也就是说该模型为什么选择逻辑斯蒂回归的密度函数和分布函数？且是个概率模型。

先来回顾下逻辑斯蒂曲线长什么样，如下图：
alt text
它是一个分布函数，函数如下

P (Y = 1 | x) = e x p ( w \cdot x + b ) 1 + e x p ( w \cdot x + b )

$P(Y = 1 | x) = \frac{exp(w\cdot x+ b)}{1 + exp(w\cdot x + b)}$
我们再令

X=w⋅x+b $X = w\cdot x + b$ ，即

X $X$ 可以表示为复合的随机变量。可为什么就说

X $X$ 的分布就符合该函数分布呢？而且我们从密度函数去分析的话，

P(Y=1|x) $P(Y = 1 | x)$ 的概率实际是

P(Y=1|−∞≤x≤X) $P(Y = 1 | -\infty \le x \le X )$ 的概率，它是对密度函数从负无穷到当前

X $X$ 的累加。为什么要这样做？

OK，我更愿把随机变量 $X$ 描述为条件。也就是要分类的一个样例达到当前条件( $X$ )下的概率是多少！我们再来看看高中生物学过的一个类似的曲线。如下图：
alt text

刚开始，种群的数量非常少，繁殖的速度会比较慢。随着数量的增加，繁殖速度越来越快，然后，食物不足，有天敌出现等原因，增速开始下降，最后稳定在一个范围内。Logistic曲线非常好的描述了这个变化规律。如果用函数去表达的话，它这公式其实就是我们的分布函数。无非现在的变量换成了t，对啊，是时间！！！我起初一直没有明白什么样的随机变量需要天然的对之前所有的状态进行累加，而时间就是这种天然的随机变量。在15分钟后，我们可以明确得知，草履虫种群的密度将处于平衡状态，且再今后也很少会发生较大幅度的变化。所以，我们就认为该种群出现的条件成熟了。而所谓的条件成熟便是我们的时间，和分类问题一样，随机变量 $X$ 无非是当前条件成熟的一个标志，是对所有过去条件的累加，逐步形成的。

因此，我们可以看出该模型下，复合随机变量 $X$ 对特征的选取有指导意义，时间是很好的一个标准，如对于某类分类问题，如质量特征m(t)，若随时间变化趋于稳定，我们便可以用该特征来做逻辑斯蒂回归，而有些变量，若跟动态变化无关，那么选取的意义不大，逻辑斯蒂回归模型的预测能力也不会太好。（日后在做实际的预测时，将做进一步验证，此处仅为自己的思考，不完全正确。）
Code Time
以下测试数据来自《机器学习实战》以及博文逻辑斯谛回归与最大熵模型- 码农场
测试数据
《机器学习实战》中给出了一个testSet.txt：
-0.017612 14.053064 0 -1.395634 4.662541 1 -0.752157 6.538620 0 -1.322371 7.152853 0 0.423363 11.054677 0 0.406704 7.067335 1 0.667394 12.741452 0 -2.460150 6.866805 1 0.569411 9.548755 0 -0.026632 10.427743 0 0.850433 6.920334 1 1.347183 13.175500 0 1.176813 3.167020 1 -1.781871 9.097953 0 -0.566606 5.749003 1 0.931635 1.589505 1 -0.024205 6.151823 1 -0.036453 2.690988 1 -0.196949 0.444165 1 1.014459 5.754399 1 1.985298 3.230619 1 -1.693453 -0.557540 1 -0.576525 11.778922 0 -0.346811 -1.678730 1 -2.124484 2.672471 1 1.217916 9.597015 0 -0.733928 9.098687 0 -3.642001 -1.618087 1 0.315985 3.523953 1 1.416614 9.619232 0 -0.386323 3.989286 1 0.556921 8.294984 1 1.224863 11.587360 0 -1.347803 -2.406051 1 1.196604 4.951851 1 0.275221 9.543647 0 0.470575 9.332488 0 -1.889567 9.542662 0 -1.527893 12.150579 0 -1.185247 11.309318 0 -0.445678 3.297303 1 1.042222 6.105155 1 -0.618787 10.320986 0 1.152083 0.548467 1 0.828534 2.676045 1 -1.237728 10.549033 0 -0.683565 -2.166125 1 0.229456 5.921938 1 -0.959885 11.555336 0 0.492911 10.993324 0 0.184992 8.721488 0 -0.355715 10.325976 0 -0.397822 8.058397 0 0.824839 13.730343 0 1.507278 5.027866 1 0.099671 6.835839 1 -0.344008 10.717485 0 1.785928 7.718645 1 -0.918801 11.560217 0 -0.364009 4.747300 1 -0.841722 4.119083 1 0.490426 1.960539 1 -0.007194 9.075792 0 0.356107 12.447863 0 0.342578 12.281162 0 -0.810823 -1.466018 1 2.530777 6.476801 1 1.296683 11.607559 0 0.475487 12.040035 0 -0.783277 11.009725 0 0.074798 11.023650 0 -1.337472 0.468339 1 -0.102781 13.763651 0 -0.147324 2.874846 1 0.518389 9.887035 0 1.015399 7.571882 0 -1.658086 -0.027255 1 1.319944 2.171228 1 2.056216 5.019981 1 -0.851633 4.375691 1 -1.510047 6.061992 0 -1.076637 -3.181888 1 1.821096 10.283990 0 3.010150 8.401766 1 -1.099458 1.688274 1 -0.834872 -1.733869 1 -0.846637 3.849075 1 1.400102 12.628781 0 1.752842 5.468166 1 0.078557 0.059736 1 0.089392 -0.715300 1 1.825662 12.693808 0 0.197445 9.744638 0 0.126117 0.922311 1 -0.679797 1.220530 1 0.677983 2.556666 1 0.761349 10.693862 0 -2.168791 0.143632 1 1.388610 9.341997 0 0.317029 14.739025 0
加载数据与数据可视化
# 记载数据 def loadDateSet(): """ 加载数据集 :return: 输入向量矩阵和输出向量 """ dataMat = [];labelMat = [] fr = open('testSet.txt') for line in fr.readlines(): lineArr = line.strip().split() dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])]) labelMat.append(int(lineArr[2])) return dataMat,labelMat # 可视化数据 import matplotlib.pyplot as plt import numpy as np def plotBestFit(weights): """ 画出数据集和逻辑斯蒂回归直线 :param weights: """ dataMat,labelMat = loadDateSet() dataArr = np.array(dataMat) n = np.shape(dataArr)[0] xcord1 = [];ycord1 = [] xcord2 = [];ycord2 = [] for i in range(n): if int(labelMat[i]) == 1: xcord1.append(dataArr[i,1]) ycord1.append(dataArr[i,2]) else: xcord2.append(dataArr[i,1]) ycord2.append(dataArr[i,2]) fig = plt.figure() ax = fig.add_subplot(111) ax.scatter(xcord1,ycord1,s=30,c='red',marker='s') ax.scatter(xcord2,ycord2,s=30,c='green') if weights is not None: x = range(-3.0,3.0,0.1) y = (-weights[0] - weights[1] * x) / weights[2] ax.plot(x,y) plt.xlable('x1');plt.ylabel('x2') plt.show()
数据呈现如下图所示：
梯度下降算法
在《机器学习实战》中，称求函数最小值的时候用的是梯度下降算法，而此处求的是对数似然函数的最大值，所以应该称为梯度上升算法。其实梯度下降算法在我的博文【算法细节系列（3）：梯度下降法，牛顿法，拟牛顿法】中也详细阐述了，这里重新复习下。
函数的梯度由其偏导数构成：

$\nabla f (x, y) = ⎛ ⎝ ⎜ ⎜ ⎜ \partial f ( x , y ) \partial x \partial f ( x , y ) \partial y ⎞ ⎠ ⎟ ⎟ ⎟$ $\nabla f(x,y) = \begin{pmatrix} \frac{\partial f(x,y)}{\partial x}\\ \frac{\partial f(x,y)}{\partial y}\\ \end{pmatrix}$
梯度是函数增长最快的方向，记移动补偿为 $\alpha$ ，则梯度算法的迭代公式为：

$w : = w + α \nabla w f (w)$ $w : = w + \alpha \nabla_wf(w)$
假定权值向量 $w$ 有了，怎么计算模型输出呢？
特征向量乘以权值向量得出一个实数z：

$z = w 0 x 0 + w 1 x 1 + \dots + w n x n$ $z = w_0x_0 + w_1x_1+\cdots+w_nx_n$
希望通过实数输出一个0或1的类别，这时候就需要利用Logistic函数了：
$δ (z) = 1 1 + e - z$ $\delta(z) = \frac{1}{1+e^{-z}}$
其图像如下：
参考《机器学习实战》，加了一些注释的代码：
from numpy import * import logistic def sigmoid(inX): return 1.0/(1+exp(-inX)) def gradAscent(dataMatIn,classLabels): """ 逻辑斯蒂回归梯度上升优化算法 :param dataMatIn: 输入X矩阵 :param classLabels: 输出Y矩阵 :return: 权值向量 """ dataMatrix = mat(dataMatIn) labelMat = mat(classLabels).transpose() m,n = shape(dataMatrix) alpha = 0.001 maxCycles = 500 weights = ones((n,1)) for k in range(maxCycles): h = sigmoid(dataMatrix*weights) error = (labelMat - h) weights += alpha * dataMatrix.transpose() * error return weights dataArr, labelMat = logistic.loadDataSet() weights = gradAscent(dataArr, labelMat) logistic.plotBestFit(weights)
分类效果如下图所示：
随机梯度上升算法
梯度下降算法在每次更新权值向量的时候都需要遍历整个数据集，该方法对小数据集尚可。但如果有数十亿样本和成千上万的特征时，它的计算复杂度就太高了。一种改进的方法是一次仅用一个样本点的回归误差来更新权值向量，这个方法叫随机梯度下降算法。由于可以在遇到新样本的时候再对分类器进行增量式更新，所以随机梯度上升算法是一个在线学习算法；与此对应，一次处理完所有数据的算法（如梯度上升算法）被称作“批处理”。
代码实现：
def stocGradAscent0(dataMatrix, classLabels, history_weight): """ 随机梯度上升算法 :param dataMatIn:输入X矩阵（100*3的矩阵，每一行代表一个实例，每列分别是X0 X1 X2） :param classLabels: 输出Y矩阵（类别标签组成的向量） :return:权值向量 """ dataMatrix = array(dataMatrix) m,n = shape(dataMatrix) alpha = 0.01 weights = ones(n) #初始化为单位矩阵 for i in range(m): h = sigmoid(sum(dataMatrix[i]*weights)) #挑选（伪随机）第i个实例来更新权值向量 error = classLabels[i] - h weights = weights + dataMatrix[i] * alpha * error history_weight.append(copy(weights)) return weights
可见随机梯度上升算法中h和error都是数值，没有进行复杂的矩阵运算。
可视化：
将原程序中的gradAscent换成stocGradAscent0，得出如下结果——
可以看到，最终拟合出来的直线效果并不如梯度上升算法，大约错了1/3的样本。
不过这种比较并不公平，毕竟随机梯度上升算法每次迭代的复杂度小得多，而且也只迭代了样本个数（200）次。
改进的随机梯度上升算法
既然随机梯度上升算法最终给出的参数不好，那是否仅仅是因为参数没有足够收敛，而算法本质是优秀的呢？对此，可以逐步减小步长，避免参数周期性的抖动。
代码实现：
def stocGradAscent1(dataMatrix, classLabels, numIter=150): """ 改进的随机梯度上升算法 :param dataMatIn:输入X矩阵（100*3的矩阵，每一行代表一个实例，每列分别是X0 X1 X2） :param classLabels: 输出Y矩阵（类别标签组成的向量） :param numIter: 迭代次数 :return: """ dataMatrix = array(dataMatrix) m,n = shape(dataMatrix) weights = ones(n) #初始化为单位矩阵 for j in range(numIter): dataIndex = range(m) for i in range(m): alpha = 4/(1.0+j+i)+0.0001 #步长递减，但是由于常数存在，所以不会变成0 randIndex = int(random.uniform(0,len(dataIndex))) #总算是随机了 h = sigmoid(sum(dataMatrix[randIndex]*weights))
可视化：
三者的收敛速度如图：
以上内容摘自博文【逻辑斯谛回归与最大熵模型- 码农场】
参考文献
机器学习中 Sigmoid 函数的物理含义
逻辑斯谛回归与最大熵模型- 码农场
算法细节系列（3）：梯度下降法，牛顿法，拟牛顿法
李航. 统计学习方法[M]. 北京：清华大学出版社，2012
Peter Harrington. Machine Learning in Action[M]. 北京：人民邮电出版社，2013

原文链接：https://blog.csdn.net/u014688145/article/details/53608203