机器学习实战笔记2

第5章 Logistic回归

原理

【实质和感知机有点像】

———————-补充：感知机————————-

定义

给定一个数据集：$T=\left\{{(x_{1},y_{1}),(x_{2},y_{2}),…,(x_{N},y_{N})}\right\}$ , 其中$x_{i}\epsilon R^{n} , y_{i}\epsilon \left \{ +1,-1\right.\},i=1,2,…N$ ,如果存在某个超平面$S$

$$wx+b=0$$

能够将数据集的正实例和负实例点完全正确划分到超平面的两侧，即对所有$y_{i}=+1$的实例$i$ ,有$wx_{i}+b>0$ ,对所有的$y_{i}=-1$ 的实例$i$ ,有$wx_{i}+b<0$ ,则称数据集T为线性可分数据集，否则，称数据集T线性不可分

学习策略

假设训练数据集是线性可分的，这样就存在一个超平面可以将训练集的正实例点和负实例点完全分开，即存在一组w,b，问题的关键就在于如何找到w,b，而这样的超平面下是没有误分类点的，我们的初始条件$w_{0},b_{0}$ 所对应的超平面肯定是有误分类点，我们需要不断的变化$w,b$ 来使得超平面左右误分类点越来越少，这样判断的标准就是误分类点

我们选择误分类点到超平面的距离总和作为损失函数:

对于点$x_{0}$ 到超平面S的距离：

​

$$\frac{1}{||w||}|wx_{0}+b|$$

对于误分类点$(x_{i},y_{i})$ ,有

$$-y_{i}(wx_{i}+b)>0$$ ,而正确分类的点必有 $y_{i}>0,wx_{i}+b>0$ ,这样我们可以挑选出误分类点

假设超平面S的误分类点集合为M，所有误分类点到超平面S的总距离为

​

$$L(w,b)=-\frac{1}{||w||}\sum_{x_{i}\epsilon M} y_{i}(wx_{i}+b)$$

我们的目的是要使$L(w)=0$ 而且$||w||$ 添加进来会使计算复杂化，至于是否会影响查找速度，可以参见SVM，此处略过

得感知机学习的损失函数为

$$L(w,b)=-\sum_{x_{i}\epsilon M}y_{i}(wx_{i}+b)$$

我们要求解使得损失函数最小w，b

使用梯度下降法（至于为什么使用梯度下降法，而不是直接求导，解方程，下面有讲到），对损失函数求梯度

$$\frac{\partial L(w,b)}{w}=-\sum_{x_{i}\epsilon M}y_{i}x_{i}$$

$$\frac{\partial L(w,b)}{b}=-\sum_{x_{i}\epsilon M}y_{i}$$

统计学习方法中使用的是随机梯度下降法，就是遇到误分类点$(x_{i},y_{i})$，单对该点求梯度，对w，b更新 (最小值，是沿梯度下降，因此是减去梯度)

$$w_{i+1}=w_{i}+\eta y_{i}x_{i}$$

$$b_{i+1}=b_{i}+\eta y_{i}$$

这里说一下，如果是梯度下降法，需要这样更新：

$$w_{i+1}=w_{i}+\eta \sum_{x_{i}\epsilon M}y_{i}x_{i}$$

$$b_{i+1}=b_{i}+\eta \sum_{x_{i}\epsilon M}y_{i}$$

具体算法及对偶形式就不叙述了。。。

【注】感知机和逻辑斯谛回归模型的一大区别在于损失函数，在梯度下降法中感知机针对的是误分类点，逻辑斯谛回归模型是对所有的点

———————-结束————————————–

给定Sigmoid函数和输入

$$Sigmoid函数： \sigma(z)=\frac{1}{1+e^{-z}}$$

$$其中 z=w_{0}x_{0}+w_{1}x_{1}+…+w_{n}x_{n}+b ,w_{i}为特征x_{i}的权重$$

统计学习方法中讲到：

二项逻辑斯谛回归模型（二分类）是如下条件概率分布：

$$P(Y=1|x)=\frac{e^{wx+b}}{1+e^{wx+b}}$$

$$P(Y=0|x)=\frac{1}{1+e^{wx+b}}$$

对于上式的w，x，b加以扩充修改

$$w=(w^{(1)},w^{(2)},…,w^{(n)},b) , x=(x^{(1)},x^{(2)},…,x^{(n)},1)$$

条件概率分布变为

$$P(Y=1|x)=\frac{e^{wx}}{1+e^{wx}}$$

$$P(Y=0|x)=\frac{1}{1+e^{wx}}$$

继续分析影响上述条件变量的因素，一个事件的几率是指该事件发生的概率与该事件不发生的概率的比值。如果事件发生的概率为p，那么该事件发生的几率是$\frac{p}{1-p}$ ，对逻辑斯谛回归而言，其对数几率为：

$$log\frac{P(Y=1|x)}{1-P(Y=1|x)}=wx$$

即输出Y=1的对数几率是输入x的线性函数，该线性函数的值越接近正无穷，概率值就越接近1；该线性函数的值越接近负无穷，概率值就越接近0，这样的模型就是逻辑斯谛回归模型

关于逻辑斯谛回归模型的参数估计

对于给定的训练数据集$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)=\pi(x)$ , $P(Y=0|x)=1-\pi(x)$

其似然函数为

$$\prod_{i=1}^{N}[\pi(x_{i})]^{y_{i}}[1-\pi(x_{i})]^{1-y_{i}}$$

对数似然函数为

​

$$L(w)=\sum_{i=1}^{N}\bigg(y_{i}log\pi(x_{i})+(1-y_{i})log(1-\pi(x_{i}))\bigg)$$

​

$$=\sum_{i=1}^{N}\bigg(y_{i}log\frac{\pi(x_{i})}{1-\pi(x_{i})}+log(1-\pi(x_{i}))\bigg)$$

这里将

$$\pi(x_{i})=\frac{e^{wx_{i}}}{1+e^{wx_{i}}}$$ 代入上式

$$L(w)=\sum_{i=1}^{N}\bigg(y_{i}wx_{i}-log(1+e^{wx_{i}})\bigg)$$

似然估计法在找到似然函数$L(w)$后，对其求极大值，得到w的估计值

【至于极大似然估计法为什么是求解最大值，百科上讲到：极大似然估计是建立在这样的思想上：已知某个参数能使这个样本出现的概率最大，我们当然不会再去选择其他小概率的样本，所以干脆就把这个参数作为估计的真实值。】

大概意思就是说按照极大似然估计的方法对似然函数求解最大值，所对应的参数是使误分类点最小的参数。。。

这种多维下求解目标函数的最优化问题不同于一二维下直接求导，令导函数等于0，然后解出参数，多维下很难求解，于是有了梯度下降法等方法

我们通过参数的梯度来知道目标函数下降(上升)最快的方向，一点点的更新参数，使函数值向着极大值迈进，直到变化量小于某个阈值，停止

对对数似然函数求导：

$$\frac{\partial L}{\partial w}=\sum_{i=1}^{N}\bigg(y_{i}x_{i}-\frac{x_{i}e^{wx_{i}}}{1+e^{wx_{i}}}\bigg)$$

$$w_{i+1}=w_{i}+\alpha \frac{\partial L}{\partial w}$$

不断的迭代w，直至$L(w_{i+1})-L(w_{i})$的差值在某一个阈值内，停止迭代，所求的w为极优参数

【这里需要注意一点，机器学习实战中使用的函数是

$$P(Y=1|x)=\frac{1}{1+e^{-wx}}$$ ,实质上是没有区别的,将统计学习方法中的函数分子分母同除 $e^{wx}$ ，便可得到此式，对于机器学习实战中对应代码的公式: 即上面的梯度括号第二项分子分母同除 $e^{wx_{i}}$ 后的结果】

实例1

现有很多二维坐标系下的点，每个点上有一个标签，代表分类。如下所示：

-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

下面使用Logistis回归来进行分类预测

首先是数据加载，训练集如上所示，前两列为特征$x^{(1)},x^{(2)}$，最后一列为类别$y$，$x$ 扩充为$(1,x^{(1)},x^{(2)})$,代码如下：

def loadDataSet():
    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])])#注意这里1的添加，这里代表的是b
        labelMat.append(int(lineArr[2]))
    return dataMat,labelMat

求解梯度，迭代数据，更新参数

def sigmoid(inX):
    return 1.0/(1+exp(-inX))

#dataMat为list，需要转化为ieNumPy矩阵
def gradAscent(dataMat,labelMat):
    dataMatrix=mat(dataMat)
    labelMat=mat(labelMat).transpose()#转置
    m,n=shape(dataMatrix)
    alpha=0.001
    maxCycles=500
    weights=ones((n,1))#生成n行1列的值为1的矩阵
    for k in range(maxCycles):
        h=sigmoid(dataMatrix*weights)
        error=(labelMat-h)
        weights=weights+alpha*dataMatrix.transpose()*error
    return weights

上面代码中的公式在前面是有推导过的，此处代码的亮点在于使用矩阵简化代码，如果用传统的计算不知多么复杂。下面详细讲解：

$$\begin{bmatrix} 1& x_{1}^{(1)} & x_{1}^{(2)} \\ 1& x_{2}^{(1)} & x_{2}^{(2)} \\ .& .& . \\ .& .& .\\1& x_{N}^{(1)} &x_{N}^{(2)} \end{bmatrix}\begin{bmatrix}w_{1}\\w_{2}\\w_{3}\end{bmatrix}=\begin{bmatrix}z_{1}\\z_{2}\\.\\.\\z_{N}\end{bmatrix}\overset{Sigmoid(y)}{\rightarrow}=\begin{bmatrix}h_{1}\\h_{2}\\.\\.\\h_{N}\end{bmatrix}$$

关于代码中的参数weights的更新：

借用前面的公式如下

$$\frac{\partial L}{\partial w}=\sum_{i=1}^{N}\bigg(y_{i}x_{i}-\frac{x_{i}e^{wx_{i}}}{1+e^{wx_{i}}}\bigg)$$

对括号内的第二项分子分母同除$e^{wx_{i}}$

$$\frac{\partial L}{\partial w}=\sum_{i=1}^{N}\bigg(y_{i}x_{i}-\frac{x_{i}}{1+e^{-wx_{i}}}\bigg)$$

$$=\sum_{i=1}^{N}(y_{i}-\frac{1}{1+e^{-wx_{i}}})x_{i}$$

$$=\sum_{i=1}^{N}(y_{i}-h_{i})\begin{bmatrix}1\\x_{i}^{(1)}\\x_{i}^{(2)}\end{bmatrix}$$

即：

$$\begin{bmatrix}w_{1}\\w_{2}\\w_{3}\end{bmatrix}_{i+1}=\begin{bmatrix}w_{1}\\w_{2}\\w_{3}\end{bmatrix}_{i}+\alpha \begin{bmatrix}1&1&.&.&1\\x_{1}^{(1)}&x_{2}^{(1)}&.&.&x_{N}^{(1)}\\x_{1}^{(2)}&x_{2}^{(2)}&.&.&x_{N}^{(2)}\end{bmatrix}\begin{bmatrix}y_{1}-h_{1}\\y_{2}-h_{2}\\.\\.\\y_{N}-h_{N}\end{bmatrix}$$

画决策边界

迭代500次获得参数weights，超平面即为weights[0]+weights[1:]*x=0

代码如下：

#分析数据，画出决策边界
def plotBestFit(dataMat,labelMat,weights):
    data0=[]
    data1=[]
    #分离数据为label为0的和为1的
    for i in range(len(labelMat)):
        if labelMat[i]==0:
            data0.append(dataMat[i])
        else:
            data1.append(dataMat[i])
    #这一块可以参考knn中的绘图部分
    plt.scatter([x[1] for x in data0], [x[2] for x in data0], c='r', label='0')  # red
    plt.scatter([x[1] for x in data1], [x[2] for x in data1], c='g', label='1')  # green
    plt.legend()

    x1=arange(-4,4,0.1)#生成-4，到3得列表，并且步长为0.1
    #w1*x1+w2*x2+b=0,解出x2
    x2=-(float(weights[0])+float(weights[1])*x1)/float(weights[2])#注意weight的格式，需要将其转为数值
    plt.plot(x1,x2)
    plt.show()

如下图所示：