Python实现逻辑回归(Logistic Regression)

1.逻辑回归概念

逻辑分类（Logistic Classification）是一种线性模型，可以表示为$y=f(w*x+b)$，w是训练得到的权重参数（Weight）；
x是样本特征数据(逻辑回归一般要求需要对x进行归一化处理，常见的做法有最大最小值归一化：(x-min(x))/(max(x)-min(x))，0均值标准化:(x-μ)/δ)；
y是对应的分类变量（注意这里的0、1、2、3只是表示对应的标称分类，并不表示具体的数值，不能进行运算处理），二分类一般用{0,1}或{-1,1}表示，前者更常用一些，方便推导，后者也可以推导出损失函数只是表达上要更复杂一些；
b是偏置（Bias）;
f称为激活函数。就是给定一批样本数据集，和样本对象所属的分类，进行建立模型。使用模型对新来的对象分类预测。

当f函数为sigmoid函数时，就是逻辑回归。即：$y=sigmoid(w*x+b)$
在这了我们做一些定义，方便下面求导：
$h_{\theta}(x)=g(\theta^{T}x)=\frac{1}{1+e^{-x}}$

sigmoid函数求导：
$g^{'}(x)=(\frac{1}{1+e^{-x}})^{'}=\frac{e^{-x}}{(1+e^{-x})^{2}}=\frac{1+e^{-x}-1}{(1+e^{-x})^{2}}\\=\frac{1}{(1+e^{-x})}-\frac{1}{(1+e^{-x})^{2}}=\frac{1}{(1+e^{-x})}(1-\frac{1}{(1+e^{-x})})\\=g(x)(1-g(x))$
请牢记以上sigmoid函数求导结论：g(x)(1-g(x))

对于二分类，分类标号为{0,1}，可以理解为正样本的输出概率。而由于类标号只能取{0,1}，取值结果是一个二项分布，概率函数服从伯努利分布。
假设：
$P(y=1|x;\theta )=h_{\theta }(x)\\P(y=0|x;\theta )=1-h_{\theta }(x)$
将上面两个公式合成一个公式为：
$P(y|x;\theta )=(h_{\theta }(x))^{y}((1-h_{\theta }(x)))^{1-y}$
此公式没有特别复杂的地方，只是很巧妙的将y=0和y=1结合成了一个公式，当y=0时，前式为1，当y=1时后式变为了1。

预测样本x的结果为真实结果y的概率为上面的公式。那么预测一批样本，预测结果全部都和真实结果一样的概率为下面的公式。这也是最大似然的思想，用样本取估计w权重的过程，逻辑回归似然函数：

$L(\theta )=\prod_{i=1}^{m}P(y_{i}|x_{i};\theta )=\prod_{i=1}^{m}(h_{\theta }(x_{i}))^{y_{i}}((1-h_{\theta }(x_{i})))^{1-y_{i}}$

其中m表示样本数量。关于为什么最大似然函数就要取乘积，可以这样理解。最大似然就是让已经发生的事情（也就是样本）在你的模型中出现的概率尽可能大。也就是如果有一个和已知结果的样本特征一样的待测对象，那个你的预测结果也要和样本的真实结果尽可能一样。而每个样本可以理解为一个独立的事件。所有事件发生的概率尽可能大，就可以用联合概率表示所有样本都发生的概率。联合概率使用乘积的形式表示。这也就有了上式乘积的形式。
上式取对数，连乘法变加法：
$l(\theta )=logL(\theta )=\sum_{i=1}^{m}(y_{i}logh_{\theta }(x_{i})+(1-y_{i})log(1-h_{\theta }(x_{i})))$

求导：
$\frac{∂}{∂\theta_{j}}l(\theta )=(y\frac{1}{g(\theta)^{T}x}-(1-y)\frac{1}{1-g(\theta)^{T}x})\frac{∂}{∂\theta_{j}}g(\theta^{T}x)\\ =(y\frac{1}{g(\theta)^{T}x}-(1-y)\frac{1}{1-g(\theta)^{T}x})g(\theta^{T}x)(1-g(\theta^{T}x))\frac{∂}{∂\theta_{j}}\theta^{T}x\\=(y(1-g(\theta^{T}x))-(1-y)g(\theta^{T}x))x_{j}\\=(y-h_{\theta}(x))x_{j}$

上式求导中第二行直接使用了sigmoid函数求导结论：g(x)(1-g(x))，对第二行中的θ^Tx求导实际上得到的就是第j个x，记做xj。

然后我们另上式求导结果等于0，你会发现它无法解析求解，在这里我们只能借用迭代法来求解，下面我们采用梯度下降算法求解。

2.梯度下降算法

在上一章节最后求导结果中，我们稍作修改，和线性回归梯度下降算法类似，很容易即可实现基于梯度下降算法的逻辑回归。
$\theta_{j}=\theta_{j} + α(y-h_{\theta}(x))x_{j}$

基于梯度上式算法求解theta，Python源代码如下：

def sigmoid(x):
    return 1/(1+np.exp(-x))

def LogRegGD(x,y):
    # 梯度上升算法
    alpha=0.05  # 步长
    max_loop=500 # 循环次数
    #x = StandardScaler().fit_transform(x) # 数据进行标准化
    x=np.mat(x)
    y=np.mat(y).transpose()
    m,n=np.shape(x)
    weights=np.ones((n,1))
    
    for k in range(max_loop):
        # k=2
        h=sigmoid(x*weights)
        error=y-h
        weights=weights+alpha*x.transpose()*error # 梯度下降算法公式
    print('k:',k,'weights:',weights.T)
    return weights.getA()

3.随机梯度下降算法

随机梯度下降就是每一个样本更新一次权值，超过样本范围就循环从第一个样本再循环，同等数据量情况下（数据量不是特别巨大时）训练速度与批量梯度下降要慢一些，但是适合在线学习，来一条数据迭代训练一次。
在上节代码中增加按一条数据循环一次，即可实现：

def LogRegSGD(x,y):
    # 随机梯度上升算法
    alpha=0.05
    max_loop=500
    #x = StandardScaler().fit_transform(x) # 数据进行标准化
    x=np.mat(x)
    y=np.mat(y).T
    m,n=np.shape(x)
    weights=np.ones((n,1))
    
    for k in range(max_loop):
        for i  in range(0,len(x)):
            # k=0;i=0
            h=sigmoid(x[i]*weights)
            error=y[i]-h[0]
            weights=weights+alpha*x[i].T*error[0]  # 随机梯度下降算法公式
            print('k:',k,'i:',i,'weights:',weights.T)
    return weights.getA()

以上代码汇总如下：

# -*- coding: utf-8 -*-
"""
 @Time    : 2018/10/26 13:13
 @Author  : hanzi5
 @Email   : **@163.com
 @File    : linear_regression_bgd.py
 @Software: PyCharm
"""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn import preprocessing
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_iris
from sklearn.pipeline import Pipeline
from sklearn.linear_model import LogisticRegression

def sigmoid(x):
    return 1/(1+np.exp(-x))

def LogRegGD(x,y):
    # 梯度上升算法
    alpha=0.05  # 步长
    max_loop=500 # 循环次数
    #x = StandardScaler().fit_transform(x) # 数据进行标准化
    x=np.mat(x)
    y=np.mat(y).transpose()
    m,n=np.shape(x)
    weights=np.ones((n,1))
    
    for k in range(max_loop):
        # k=2
        h=sigmoid(x*weights)
        error=y-h
        weights=weights+alpha*x.transpose()*error # 梯度下降算法公式
    #print('k:',k,'weights:',weights.T)
    return weights.getA()

def LogRegSGD(x,y):
    # 随机梯度上升算法
    alpha=0.05
    max_loop=500
    #x = StandardScaler().fit_transform(x) # 数据进行标准化
    x=np.mat(x)
    y=np.mat(y).T
    m,n=np.shape(x)
    weights=np.ones((n,1))
    
    for k in range(max_loop):
        for i  in range(0,len(x)):
            # k=0;i=0
            h=sigmoid(x[i]*weights)
            error=y[i]-h[0]
            weights=weights+alpha*x[i].T*error[0]  # 随机梯度下降算法公式
            #print('k:',k,'i:',i,'weights:',weights.T)
    return weights.getA()

# 对新对象进行预测
def predict1(weights,testdata):
    testdata.insert(0, 1.0)       #现在首部添加1代表b偏量
    testMat = np.mat([testdata])
    z = testMat * np.mat(weights)
    y=sigmoid(z)
    if y>0.5:
    # if z>0:    #h>0.5的判断等价于  z>0
        return 1,y
    else:
        return 0,y
    

if __name__ == "__main__":
    # 读入100条数据，数据文件原名“testSet.txt”
    data=[
[-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]
]
    data_array=np.array(data) # 转换为numpy array类型
    x = data_array[:, :-1]    # 除最后一列，其他列作为x
    a= np.ones(len(x))        # 新增一列，默认值设为1，作为截距b的初始化值
    x=np.insert(x, 0, values=a, axis=1)
    
    y = data_array[:,-1] # 分类变量y
    
    # 使用sklearn包实现逻辑回归
    # x = StandardScaler().fit_transform(x)
    # lr = LogisticRegression()   # Logistic回归模型
    # lr.fit(x, y.ravel())        # 根据数据[x,y]，计算回归参数
    # 等价形式
    lr = Pipeline([('sc', StandardScaler()),('clf', LogisticRegression()) ])
    lr.fit(x, y.ravel())
    
    # 1、梯度下降算法实现逻辑回归
    weights1=LogRegGD(x,y) 
    type,y_hat1 = predict1(weights1,[0.317029,14.739025])
    print('梯度下降算法实现逻辑回归',type,y_hat1)
    
    # 2、随机梯度下降算法实现逻辑回归
    weights2=LogRegSGD(x,y)
    type,y_hat2 = predict1(weights2,[0.317029,14.739025])
    print('随机梯度下降算法实现逻辑回归',type,y_hat2)

参考资料：
1、Machine-Learning-With-Python
2、《机器学习实战》Peter Harrington著
3、《机器学习》西瓜书，周志华著
4、 斯坦福大学公开课 ：机器学习课程
5、机器学习视频，邹博
6、python机器学习案例系列教程——逻辑分类/逻辑回归LR/一般线性回归（softmax回归）
7、逻辑回归----Python实现
8、逻辑回归python实现
9、机器学习算法与Python实践之（七）逻辑回归（Logistic Regression）