机器学习实战读书笔记(4)--logistic回归

Logistic回归

假设我们有一些数据点,用一条直线对这些点进行拟合,这个拟合的过程成为回归.在计量经济学中我们大量的使用过线性回归,线性回归的模型试图得到一个通过属性的线性组合来进行预测的函数,即

$$f(x)=W^TX+b$$
式中的W和b通过学习获得,模型就被标定了.而logistic回归的思想是,回归的结果不是一个实值,而是一个二分类的概率.这个变量到函数值的映射关系应满足这样一些关系:（1）我们需要Y代表的是概率即Y∈（0,1）。（2）我们需要X各维度叠加和在0附近变化幅度比较大，并且是非线性的变化。而在很大或很小的时候，几乎不变化，这是基于概率的一种认识与需要。感性的一个例子，想想你学习努力的程度与从60分提高到80分和80提高到100分并不是线性的。（3）这个关系的公式要在之后形成的cost function是凸函数。

为了从属性映射到概率,我们首先想到的是一个分段函数

$$y(x)= \begin{cases} 1 & &X<0 \\ 0 & &X>0 \end{cases}$$
但是这个函数并不连续,不能满足第三个条件,这时出现了sigmoid函数,其形式如下:
$$y=\frac{1}{1+e^{(z)}}$$
这个函数的图像如下:

x = arange(-50,50,1)
y = 1.0 / (1 + exp(-x))
plt.figure()
plt.plot(x,y)
ax = plt.gca()
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.xaxis.set_ticks_position('bottom')
# set_position中的参数元组的第二个值可取-1，0，1分别代表相对‘data’的不同的位置
ax.spines['bottom'].set_position(('data', 0.5))
ax.yaxis.set_ticks_position('left')
ax.spines['left'].set_position(('data', 0))
plt.show()

<matplotlib.figure.Figure at 0x7fd6ce41ada0>






[<matplotlib.lines.Line2D at 0x7fd6ce202550>]

这个函数满足我们上文中三个条件的要求,并且有较好的数值表现,效果接近理想的分段函数,但是又是连续的.所以,我们可以使用这个函数将特征映射到二分类的概率.在线性回归中,我们使用最小二乘法求解参数,在logistic函数的参数标定的目标函数是什么呢?

构造最优化问题的目标函数

考虑单调可谓函数g(x),令

$$y=g(w^T+b)$$
这样的模型成为广义线性模型,其中函数g(x)称为联系函数.在logistic回归中,这个函数就是sigmoid函数.对其线性边界情况,边界形式就是 $W^TX$.

构造预测函数为:

$$h_{\theta}(X)=g(W^TX)=\frac{1}{1+e^{(W^TX+b)}}$$
函数 $h_{\theta}(X)$有特殊含义:它表示结果取1的概率,因此对于输入X分类结果为类别1和类别0的概率分别为:
$$P(y=1|x;\theta)=h_{\theta}(X)=(h_{\theta}(X))^y\\ P(y=0|x;\theta)=1-h_{\theta}(X)=(1-h_{\theta}(X))^(1-y)$$
上面两个式子是互斥的所以可以统一成一个式子,即:
$$P(y|x;\theta)=(h_{\theta}(X))^y (1-h_{\theta}(X))^(1-y)$$
对于所有的X,目标是使P(X)尽可能大,对所有数据取似然函数
$$L(\theta)=\prod_{i=1}^m(h_{\theta}(X))^y_i (1-h_{\theta}(X))^(1-y_i)$$
对两边取对数化简,有
$$l(\theta)=\sum_{i=1}^m(y_i log h_{\theta}(X))+ (1-y_i)log (1-h_{\theta}(X))$$
这就是我们的第一版目标函数,对数似然函数最大,考虑规划问题的标准形式和平均概率的问题,对其乘以一个系数-1/m,有
$$min J(\theta)=-\frac{1}{m}l(\theta)$$
求得最佳参数 $\theta$,完成回归函数的标定.

目标函数的一种求解方法–梯度下降法

有了目标函数,如何求解,在非线性规划问题中一个最基础的求解办法就是梯度下降法,尽管梯度下降法存在一些慢收敛的问题,但是还是可以作为讨论的第一个方法,学到一些有趣的东西.

首先复习一下梯度的概念.对一个函数$f(x,y)$求梯度$\nabla f(x,y)$,即是对x,y分别求偏导,构成一个向量.

$$\nabla f(x,y)= \begin{bmatrix} \frac{\partial f(x,y)}{ \partial x} \\ \frac{\partial f(x,y)}{ \partial y} \end{bmatrix}$$
那就是对 $log h_{\theta}(X))+ (1-y_i)log (1-h_{\theta}(X))$求导嘛,这里面有一个较重要的基础知识点是h(x)函数是sigmoid函数,sigmoid函数求导很有意思,也很重要(之后神经网络中常用的激活函数),该函数的求导h’(x)=h(x)[1-h(x)],具体自己去推导,直接给结论,目标函数的梯度中第i项为 $\sum_{i=1}^m (h(\theta)-y_i)x_i$,如果令 $Error=h(\theta)-y_i$那么梯度的矩阵表示就变成了 $ErrorX^T$,代码中直接使用了结论,但是为了不做调包侠,这些要懂,这部分有点难理解,一定要自己推导,但是反过来发现,这最终的结果的形式和线性回归的目标函数的梯度的形式是一样的啊!

def loadDataSet():
    '''这是个准备数据集的小函数'''
    dataMat = []; labelMat = []
    fr = open('testSet.txt')
    for line in fr.readlines():
        lineArr = line.strip().split('\t')
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelMat.append(int(lineArr[2]))
    return dataMat, labelMat

def sigmoid(inX):
    '''构造sigmoid函数'''
    return 1.0/(1 + exp(-inX))


def gradAscent(dataMatIn, classLabels):
    '''梯度下降法计算weight, 没有推导过程只有结果
    主要数据都是mat格式, 方便进行矩阵乘法'''
    dataMatrix = mat(dataMatIn)
    labelMatrix = mat(classLabels).transpose()
    m,n = shape(dataMatrix)
    weights = ones((n,1))
    MaxCycle = 500
    alpha = 0.001
    for i in range(MaxCycle):
        # result = dataMatrix * weights
        result = sigmoid(dataMatrix * weights)
        # 注意两个问题,一个是sigmoid函数f(z)的求导是f(z)[1-f(z)]
        # 回归的目标函数不是简单的方差最小,由于这是概率问题,需要采用最大似然函数作为目标函数
        Error = labelMatrix - result
        weights = weights + alpha * dataMatrix.T * Error
    return weights