logistics回归分析+代码详解

1.2 Logistic 回归的引入

1.3 分类预测函数问题的转化成求θ

1.4 梯度下降法求解θ

2.3 利用随机梯度上升法训练样本

1.2 Logistic 回归的引入

Logistic 回归是概率非线性模型，用于处理二元分类结果，名为回归实为分类。下面给出一个二元分类的例子，如下图所示：

图中数据分成两类，打叉的一类用 y = 1 表示，另一类为圆圈用 y= 0 表示。现在我们要对数据进行分类，要找到一个分类函数(边界线)，我们很容易得出一个线性函数对其分类：z = θ0 + θ1x1 + θ2x2 。但我们想要的函数应该是，能接受所有的输入然后预测类别。例如：在两个分类的情况下，函数输出 0 或 1。因此，我们就需要引入Sigmoid 函数，这个函数的性质可以满足要求。Sigmoid 函数：

Sigmoid 函数的值域为（0,1），且具有良好的从0 到 1 的跳跃性，如在两个不同的坐标尺度下的函数图形：

所以，我们把线性方程和Sigmoid 函数结合起来就能解决问题。即：分类预测函数 hθ (x) = g( θ0 + θ1x1 + θ2x2) .我们就可以对样本数据进行分类，如下图所示：

其中，θT 是向量 θ (θ0, θ1,... ,θn) 的转置，向量 x是 ( x0 , x1 ,... , xn)，x0 =1，这是便于计算。

1.3 分类预测函数问题的转化成求θ

通过上面的分析，我们得出了分类预测函数 hθ(x) , 但其中向量 x 是已知的，向量 θ 未知，即我们把求分类函数问题转化成求向量 θ 。
因为Sigmoid 函数的取值区间（0,1），那我们可以看做概率 P（y = 1 | xi ; θ）= hθ(x) , 表示在 xi 确定的情况下，类别 y = 1 的概率。
由此我们也可以得出在 xi 确定的情况下，类别 y = 0 的概率 P（y = 0 | xi ; θ）= 1 - P（y = 1 | xi ; θ）= 1 - hθ(x) . 即：

我们可以将这两个式子合并得：(其中的 y = 0 或 1 )：类别为0或1 的概率函数

这时候我们可以利用最大似然函数对向量 θ 求值，可以理解为选取的样本数据的概率是最大的，

①那么样本数为 m 的似然函数为：（概率分布函数的联合分布密度）

②通过对数的形式对似然函数进行变化，对数似然函数：

此时我们不求解似然方程，而是利用梯度下降法。
这里的最大似然函数的值就是我们所要求的向量 θ , 而求解的方法利用梯度下降法。

1.4 梯度下降法求解θ

在用梯度下降法时，我们将会利用Sigmoid 函数的一个性质： g, (z) = g(z）[ 1- g(z) ] 。

构造一个Cost函数（损失函数），该函数表示预测的输出（h）与训练数据类别（y）之间的偏差。综合考虑所有训练数据的“损失”，将Cost求和或者求平均，记为J(θ)函数（代价函数），表示所有训练数据预测值与实际类别的偏差。

损失函数：

J(θ)代价函数：

其中，x(i) 每个样本数据点在某一个特征上的值，即特征向量x的某个值，y(i) 是类别号，m 是样本对象个数。

梯度下降法含义：

梯度下降法，就是利用负梯度方向来决定每次迭代的新的搜索方向，使得每次迭代能使待优化的目标函数逐步减小。梯度其实就是函数的偏导数。

这里对用梯度下降法对 J (θ) 求最小值，与求似然函数的最大值是一样的。则 J(θ) 最小值的求解过程：

其中 α 是步长。
具体步骤为：

即

这就是一开始，我对代码中公式困惑的地方。在这里我在补充一点，以上的梯度下降法可以认为是批量梯度下降法（Batch Gradient Descent），由于我们有m个样本，这里求梯度的时候就用了所有m个样本的梯度数据。下面介绍随机梯度下降法（Stochastic Gradient Descent）：

随机梯度下降法，其实和批量梯度下降法原理类似，区别在与求梯度时没有用所有的m个样本的数据，而是仅仅选取一个样本j来求梯度。随机梯度下降法，和批量梯度下降法是两个极端，一个采用所有数据来梯度下降，一个用一个样本来梯度下降。自然各自的优缺点都非常突出。
对于训练速度来说，随机梯度下降法由于每次仅仅采用一个样本来迭代，训练速度很快，而批量梯度下降法在样本量很大的时候，训练速度不能让人满意。
对于准确度来说，随机梯度下降法用于仅仅用一个样本决定梯度方向，导致解很有可能不是最优。
对于收敛速度来说，由于随机梯度下降法一次迭代一个样本，导致迭代方向变化很大，不能很快的收敛到局部最优解。公式为：

到这里已经把Logistics 回归的原理推导完成。这里对以上推导做一次总结：

二分类边界线 —要利用—> Sigmoid 函数 ——>要求向量 θ —在 xi 确定的情况下求类别 y = 0、1 的概率—> P(y | x ; θ) —— > 似然函数 —得到—> 损失函数——>代价函数 J (θ) ——> 梯度下降法求解向量 θ ——> 最终公式 θj

步骤简略公式为：

2 使用python 实现Logistic 回归

from numpy import *
import matplotlib.pyplot as plt######################################################读取数据############################################################
def loadDataSet():'''数据集的前两个值分别为X1和X2,第三个值是数据对应的类别标签，为了方便计算，把X0的值设置成了1.0'''dataMat = []labelMat = []fr = open('testSet.txt')for line in fr.readlines():lineArr = line.strip().split()  #按行读取并按空格拆分，每一行为一个列表， ['-0.017612', '14.053064', '0']# 为了方便计算，我们将 X0 的值设为 1.0 ，也就是在每一行的开头添加一个 1.0 作为 X0# 因为线性回归化式为 H(x) = W0 + W1*X1 + W2*X2，即为 (W0, W1, W2)*(1.0, X1, X2),其中 (W0, W1, W2) 即为所求回归系数 W，所以模拟了一个数据列全为1.0dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])#[[1.0, -0.017612, 14.053064], [1.0, -1.395634, 4.662541],........]含有100个元素的列表，每一个元素都为一个列表，包含3个元素labelMat.append(int(lineArr[2])) #100个元素（类别标签）的列表return dataMat, labelMat################################################分类器的分类（转换）函数#####################################################
def sigmoid(inX):return 1.0 / (1 + exp(-inX))################################### 梯度上升算法，用来计算出最佳回归系数w(weights)#############################################
#输入数据特征与数据的类别标签，dataMatIn存放的是100*3的矩阵"""第一个参数(dataMatIn)是3维数组，每列代表每个不同特征，每行代表每个训练样本第二个参数(labelMat)是类别标签，1*100的行向量，为便于计算，将其转换为列向量，即进行转置，并赋值给labelMat转化矩阵[[0,1,0,1,0,1.....]]为[[0],[1],[0].....]"""
def gradAscent(dataMatIn, classLabels):dataMatrix = mat(dataMatIn) #转为100*3的矩阵 [ 1.0000000e+00 -1.7612000e-02  1.4053064e+01]labelMat = mat(classLabels).transpose()  #转为100*1的矩阵，transpose() 行列转置函数,将行向量转化为列向量   =>  矩阵的转置# m->数据量、样本数100    n->特征数3m, n = shape(dataMatrix)# shape()函数是numpy.core.fromnumeric中的函数，它的功能是查看矩阵或者数组的维数alpha = 0.001  # 步长，向函数增长最快的方向的移动量，即学习率maxCycles = 500# 迭代次数weights = ones((n, 1))  # 转为3*1的矩阵，元素为1的矩阵赋给weihts，即回归系数初始化为1# 循环 maxCycles次, 每次都沿梯度向真实值 labelMat 靠拢for k in range(maxCycles):# 求当前sigmoid函数的值h = sigmoid(dataMatrix * weights)  # 矩阵相乘 包含了300次的乘积 [100*3] * [3*1] = [100*1]error = (labelMat - h)  # 向量减法，计算真实类别与预测类别的差值，h是一个列向量，列向量的元素个数等于样本数，即为100weights = weights + alpha * dataMatrix.transpose() * error  # 矩阵相乘，dataMatrix.transpose()* error 就是梯度f(w)，按照该差值的方向调整回归系数return weights###################################################分析数据：画出决策边界####################################################
# 画出数据集和Logistic回归最佳拟合直线
def plotBestFit(weights):dataMat, labelMat = loadDataSet()dataArr = array(dataMat)n = shape(dataArr)[0]   # n->数据量、样本数# xcord1,ycord1代表正例特征1,正例特征2# xcord2,ycord2代表负例特征1,负例特征2xcord1 = []; 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()# “111”表示“1×1网格，第一子图”，“234”表示“2×3网格，第四子图”。ax = fig.add_subplot(111)ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')ax.scatter(xcord2, ycord2, s=30, c='green')# 画线# 设定边界直线x和y的值,并且以0.1为步长从-3.0到3.0切分。x = arange(-3.0, 3.0, 0.1)#[-3.000000000 -2.90000000 -2.80000000 -2.70000000"""y的由来？首先理论上是这个样子的:dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])w0*x0+w1*x1+w2*x2=f(x)x0最开始就设置为1， x2就是我们画图的y值。0是两个类别的分界处所以： w0+w1*x+w2*y=0 => y = (-w0-w1*x)/w2   """y = (-weights[0] - weights[1] * x) / weights[2]  #最佳拟合直线ax.plot(x, y)plt.xlabel('X1');plt.ylabel('X2');plt.show()# 输出运用梯度上升优化算法后得到的最理想的回归系数的值
if __name__ =='__main__':dataArr,labelMat=loadDataSet()weights=gradAscent(dataArr,labelMat)plotBestFit(weights.getA())

数据集为

-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

2.3 利用随机梯度上升法训练样本

from numpy import *
import matplotlib.pyplot as plt
import numpy as np######################################################读取数据############################################################
def loadDataSet():'''数据集的前两个值分别为X1和X2,第三个值是数据对应的类别标签，为了方便计算，把X0的值设置成了1.0'''dataMat = []labelMat = []fr = open('testSet.txt')for line in fr.readlines():lineArr = line.strip().split()  #按行读取并按空格拆分，每一行为一个列表， ['-0.017612', '14.053064', '0']# 为了方便计算，我们将 X0 的值设为 1.0 ，也就是在每一行的开头添加一个 1.0 作为 X0# 因为线性回归化式为 H(x) = W0 + W1*X1 + W2*X2，即为 (W0, W1, W2)*(1.0, X1, X2),其中 (W0, W1, W2) 即为所求回归系数 W，所以模拟了一个数据列全为1.0dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])#[[1.0, -0.017612, 14.053064], [1.0, -1.395634, 4.662541],........]含有100个元素的列表，每一个元素都为一个列表，包含3个元素labelMat.append(int(lineArr[2])) #100个元素（类别标签）的列表return dataMat, labelMat################################################分类器的分类（转换）函数#####################################################
def sigmoid(inX):return 1.0 / (1 + exp(-inX))def stocGradAscent1(dataMatrix, classLabels, numIter=150):dataMatrix = np.array(dataMatrix)m, n = np.shape(dataMatrix)  # 返回dataMatrix的大小。m为行数,n为列数。weights = np.ones(n)  # 参数初始化for j in range(numIter):dataIndex = list(range(m))for i in range(m):alpha = 4 / (1.0 + j + i) + 0.01  # 降低alpha的大小，每次减小1/(j+i)。randIndex = int(random.uniform(0, len(dataIndex)))  # 随机选取样本h = sigmoid(sum(dataMatrix[randIndex] * weights))  # 选择随机选取的一个样本，计算herror = classLabels[randIndex] - h  # 计算误差weights = weights + alpha * error * dataMatrix[randIndex]  # 更新回归系数del (dataIndex[randIndex])  # 删除已经使用的样本return weights  # 返回###################################################分析数据：画出决策边界####################################################
# 画出数据集和Logistic回归最佳拟合直线
def plotBestFit(weights):dataMat, labelMat = loadDataSet()dataArr = array(dataMat)n = shape(dataArr)[0]   # n->数据量、样本数# xcord1,ycord1代表正例特征1,正例特征2# xcord2,ycord2代表负例特征1,负例特征2xcord1 = []; 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()# “111”表示“1×1网格，第一子图”，“234”表示“2×3网格，第四子图”。ax = fig.add_subplot(111)ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')ax.scatter(xcord2, ycord2, s=30, c='green')# 画线# 设定边界直线x和y的值,并且以0.1为步长从-3.0到3.0切分。x = arange(-3.0, 3.0, 0.1)#[-3.000000000 -2.90000000 -2.80000000 -2.70000000"""y的由来？首先理论上是这个样子的:dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])w0*x0+w1*x1+w2*x2=f(x)x0最开始就设置为1， x2就是我们画图的y值。0是两个类别的分界处所以： w0+w1*x+w2*y=0 => y = (-w0-w1*x)/w2   """y = (-weights[0] - weights[1] * x) / weights[2]  #最佳拟合直线ax.plot(x, y)plt.xlabel('X1');plt.ylabel('X2');plt.show()# 输出运用梯度上升优化算法后得到的最理想的回归系数的值
if __name__ =='__main__':dataArr,labelMat=loadDataSet()weights=stocGradAscent1(dataArr,labelMat)plotBestFit(weights)