ML刻意练习第二周之逻辑斯特回归

一、逻辑斯特回归和Sigmoid函数的分类概述
我们想要的函数可以接受所有的输入然后预测出类别。之前接触过的阶跃函数便有这种性质，但是阶跃函数在x=0除是不可导的，数学上的计算很难处理，因此我们采用也有类似性质的sigmoid函数，其计算公式如下：

其图形如下所示：

二、基于最优化方法的最佳回归系数确定
sigmoid函数的输入记为z，由下面公式给出：

我们需要得到的结果就是n个系数wi，i = 1,2,…,n。

1.梯度上升法
因为梯度是函数增长最快的方向，因此可以利用梯度来快速得到函数的最值。
其中梯度上升法用来求函数的最大值，迭代公式如下：

而之前学过的梯度下降法是来求函数的最小值，迭代公式只需将上面的加法变成减法即可。

2.训练算法：使用梯度上升找到最佳参数
首先，加载文件中的数据集特征和标签。

def loadDataSet():
    dataMat = []
    labelMat = []
    fr = open('E:\学习资料\机器学习算法刻意练习\机器学习实战书电子版\machinelearninginaction\Ch05\\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]))
    fr.close()
    return dataMat,labelMat

注意：程序中strip函数的用法如下：
s为字符串，rm为要删除的字符序列。
s.strip(rm) 删除s字符串中开头、结尾处，位于 rm删除序列的字符
s.lstrip(rm) 删除s字符串中开头处，位于 rm删除序列的字符
s.rstrip(rm) 删除s字符串中结尾处，位于 rm删除序列的字符
当rm为空时，默认删除空白符（包括’\n’, ‘\r’, ‘\t’, ’ ')

之后定义sigmoid函数，因为此时是矩阵运算，需要用numpy模块下的exp指数函数：

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

接下来实现梯度上升法：

def gradAscent(dataMatIn, classLabels):        #梯度上升算法（只规定了迭代次数，没有规定误差小到什么程度退出）
    dataMatrix = np.mat(dataMatIn)             #转换为numpy下的矩阵,m*n
    labelMat = np.mat(classLabels).transpose() #转换为numpy下的矩阵并取转置变为列向量，m*1
    m,n = np.shape(dataMatrix)                 #m行数（数据个数）,n列数（特征个数+1）
    alpha = 0.001                              #步长取0.001
    maxCycles = 500                            #最多迭代500次
    weights = np.ones((n,1))                   #存储最后的系数Wi，列向量n*1
    for k in range(maxCycles):                 #迭代过程
        h = sigmoid(dataMatrix*weights)        #m*1阶矩阵
        error = (labelMat - h)                 #与最终标签的差值，m*1阶
        weights = weights + alpha * dataMatrix.transpose()* error
    return weights

注意：该过程中要特别注意矩阵的阶数，即是否取转置。此处必须使用mat函数来生成矩阵，因为mat生成的矩阵中星号*才表示矩阵的乘法，array生成的矩阵中星号*表示对应元素的乘积。

运行结果如下：

[[ 4.12414349]
 [ 0.48007329]
 [-0.6168482 ]]

分别对应三个W的值。

3.分析数据：画出决策边界
本章节将画出所有的数据与最后得出的分隔线，从而让优化过程更容易理解。

def plotBestFit(weights):
    dataMat,labelMat = loadDataSet()
    dataArr = np.array(dataMat)#这步是为了什么，为什么特意将mat型矩阵转化为array呢
    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')
    x = np.arange(-3.0, 3.0, 0.1)#-3.0到3.0，步长为0.1的数字形成的列表，左闭右开
    y = (-np.array(weights[0]) - weights[1]*x)/weights[2]#因为weights[0] + weights[1]*x + y *weights[2] = 0
    #weight[0]，weight[1]，weight[2]都是1*1矩阵， 所以weight[1]*x是 1*60的矩阵
    #weight[0]在作减法时会自动填充为1*60的矩阵，weight[2]也是
    ax.plot(x, y)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

其图形如下所示：

由图形可见，把两类别的数据分别标成了红色和绿色，蓝色的直线为通过逻辑斯特回归得到的最终分界线，可以大致将两类数据分开。

注意：
①.其参数weights需要用array做一下转换，否则weights[0]是只有一个元素的1*1矩阵，如果array做过转换后weights[0]就是一维列表了。
②.

 x = np.arange(-3.0, 3.0, 0.1)#-3.0到3.0，步长为0.1的数字形成的列表，左闭右开

以上函数是取一个列表，列表内是从-3.0到3.0，步长为0.1的数字，左闭右开。

③.列表不能直接取负，转换成矩阵就可以了，例如weights[0]是列表，不能用-weights[0]取负，需要转换成列表后取-np.array(weights[0])。

4.训练算法：随机梯度上升
因为梯度上升算法在每次更新回归系数的时候都需要遍历数据集，因此当数据集很大时其复杂度会很高，所以我们对其做以下改进：每次只用一个样本点来更新回归系数，成为随机梯度上升。
如果把梯度上升比喻成数据的批处理，那么随机梯度上升就好比是数据的在线处理。
改进后的随机梯度上升算法如下所示：

def stocGradAscent0(dataMatrix, classLabels):
    m,n = np.shape(dataMatrix)#m*n
    alpha = 0.01
    weights = np.ones(n)   #1*3
    for i in range(m):
        h = sigmoid(sum(dataMatrix[i]*weights))
        error = classLabels[i] - h
        weights = weights + alpha * np.array(dataMatrix)[i]* error
    return weights

注意运算过程中类型的转换，特别是列表和矩阵，因为好多可以在矩阵中运行的操作不能作用于列表，比如取负和乘小数。

改进后的完成代码如下所示：

import numpy as np
import matplotlib.pyplot as plt

def loadDataSet():
    dataMat = []
    labelMat = []
    fr = open('E:\学习资料\机器学习算法刻意练习\机器学习实战书电子版\machinelearninginaction\Ch05\\testSet.txt')
    for line in fr.readlines():
        lineArr = line.strip().split()         #根据空格切分原来的长字符串
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])#最开始的1.0是最后的常数项
        labelMat.append(int(lineArr[2]))
    fr.close()
    return dataMat,labelMat

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

# def gradAscent(dataMatIn, classLabels):        #梯度上升算法（只规定了迭代次数，没有规定误差小到什么程度退出）
#     dataMatrix = np.mat(dataMatIn)             #转换为numpy下的矩阵,m*n
#     labelMat = np.mat(classLabels).transpose() #转换为numpy下的矩阵并取转置变为列向量，m*1
#     m,n = np.shape(dataMatrix)                 #m行数（数据个数）,n列数（特征个数+1）
#     alpha = 0.001                              #步长取0.001
#     maxCycles = 500                            #最多迭代500次
#     weights = np.ones((n,1))                   #存储最后的系数Wi，列向量n*1
#     for k in range(maxCycles):                 #迭代过程
#         h = sigmoid(dataMatrix*weights)        #m*1阶矩阵
#         error = (labelMat - h)                 #与最终标签的差值，m*1阶
#         weights = weights + alpha * dataMatrix.transpose()* error
#     return weights

def stocGradAscent0(dataMatrix, classLabels):
    m,n = np.shape(dataMatrix)#m*n
    alpha = 0.01
    weights = np.ones(n)   #1*3
    for i in range(m):
        h = sigmoid(sum(dataMatrix[i]*weights))
        error = classLabels[i] - h
        weights = weights + alpha * np.array(dataMatrix)[i]* error
    return weights

def plotBestFit(weights):
    dataMat,labelMat = loadDataSet()
    dataArr = np.array(dataMat)#这步是为了什么，为什么特意将mat型矩阵转化为array呢
    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')
    x = np.arange(-3.0, 3.0, 0.1)#-3.0到3.0，步长为0.1的数字形成的列表，左闭右开
    y = (-np.array(weights[0]) - weights[1]*x)/weights[2]#因为weights[0] + weights[1]*x + y *weights[2] = 0
    #weight[0]，weight[1]，weight[2]都是1*1矩阵， 所以weight[1]*x是 1*60的矩阵
    #weight[0]在作减法时会自动填充为1*60的矩阵，weight[2]也是
    ax.plot(x, y)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

dataMat ,labelMat = loadDataSet()
# stocGradAscent0(dataMat, labelMat)
weights  = stocGradAscent0(dataMat, labelMat)
plotBestFit(np.array(weights))

运行结果如下图所示：

由此可见该算法在数据不够大时可能出现欠拟合现象。