一、最优化算法与回归介绍

最优化算法的体现：生活中遇到过很多最优化问题，例如从怎么A走到B耗时最短，如何配比调料才能做出好菜，可见最优化的作用是很大的。
回归：假设现在有一些数据点，我们用一条直线对这些点进行拟合，这个拟合过程就称作回归。

二、logistic回归的一般过程

收集数据：采用任意方式收集数据。
准备数据：数据类型为数值型，最好将其预处理为结构化数据格式。
分析数据：采用任意方法对数据进行分析。
训练算法：大部分时间将用于训练，训练的目的是为了找到最佳的分类回归系数。
测试算法：一旦训练步骤完成，分类将会很快。
使用算法：首先，我们需要输入一些数据，并将其转换成对应的结构化数值；接着，基于训练好的回归系数就可以对这些数值进行简单的回归运算，判定他们属于哪个类别；在这之后，我们就可以在输出的类别上做一些其他分析工作。

三、线性回归与对数线性回归

线性回归模型和对数线性回归模型是为了确定基于最优化方法的最佳回归系数。

线性回归

线性模型的一般形式为：
$f ( x )= w_{1}x_{}1+w_{2}x_{2}+...+w_{d}x_{d} + b$
其中x=(x1, x2, …, xd)是由d维属性描述的样本，其中 xi是 x 在第 i 个属性上的取值。
线性回归的目的：学习一个线性模型以尽可能准确地预测实值输出标记：
$f(x)=wx_{i}+b$
使得
$f(x_{i})\simeq y_{i}$

对数线性回归

设 $y=g(x)=e^{z}$ ,z=wx+b,两边求对数，就可以得到对数线性回归模型：
$l n y = w x + b$

四、最小二乘法

如上所述，使用线性回归模型来预测样本的值，这其中必然不可能百分百与样本值重合，故此，我们引入最小二乘法来约束：回归直线应满足全部观测值与对应的回归估计值的误差平方和最小。用最小二乘法来计算，推导过程如下：
在这里插入图片描述

五、sigmoid函数

我们想要的函数应该是，能接受所有的输入然后预测出类别。例如，在两个类的情况下只有0和1.单位跃迁函数可以做到这一点，但是它也有一个问题：在跳跃点上从0瞬间跳跃到1，这个瞬间有时很难处理。而sigmoid函数则更加容易。sigmoid函数的公式如下：
$\sigma (x)=\frac{1}{1+e^{-x}}$
在这里插入图片描述
当x=0时，sigmoid函数值为0.5，随着x的增大或减小，sigmoid函数分别趋近于0和1，已经很接近跃迁函数了。因此，为了实现logistic分类器，我们可以在每个特征上都乘以一个回归系数，然后把所有的结果相加，将结果带入到sigmoid函数中，大于0.5的值分类为1，小于则分类为0.这样就将线性回归模型和sigmoid函数结合在了一起。
现在的问题变成了：最佳回归系数是多少？

六、梯度上升法

梯度上升发基于的思想是：要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。如果梯度记为 $\triangledown$ ,则函数f(x,y)的梯度由下式表示：

在这里插入图片描述
这个梯度意味着要沿x的方向移动 $\frac{\partial f(x,y)}{\partial x}$ ，沿y的方向移动 $\frac{\partial f(x,y)}{\partial y}$ 。
梯度上升算法到达每个点后都会重新估计移动的方向。从p0开始，计算完该点的梯度，函数就根据梯度移动到下一点p1。在p1点，梯度再次被重新计算，并沿新的梯度方向移动到p2.如此循环迭代，直到满足停止条件。迭代过程中，梯度算子总是保证我们能选取到最佳的移动方向。
在这里插入图片描述
梯度上升算法的迭代公式如下：
$w:=w+\alpha \bigtriangledown w f(w)$

七、训练算法

def loadDataSet():      # 打开testSet文件并读取，前两列为特征x1和x2，最后一列为类别标签
    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])])
        labelMat.append(int(lineArr[2]))
    return dataMat, labelMat


def sigmoid(inX):       # 分类函数
    return 1.0 / (1 + exp(-inX))


def gradAscent(dataMatIn, classLabels):     # 梯度上升寻找最佳参数w
    dataMatrix = mat(dataMatIn)
    labelMat = mat(classLabels).transpose()     # 类别标签矩阵转置
    m, n = shape(dataMatrix)
    alpha = 0.001   # 步长
    maxCycles = 500     # 迭代次数
    weights = ones((n, 1))
    for k in range(maxCycles):
        h = sigmoid(dataMatrix * weights)       # sigmoid函数计算出预测值
        error = (labelMat - h)      # 计算真是类别和预测类别的差值
        weights = weights + alpha * dataMatrix.transpose() * error      #梯度上升公式，将100*3的dataMatrix转置后与100*1的error相乘
    return weights      # 返回回归系数

运行结果为系数矩阵：

在这里插入图片描述

八、画出决策边界

def plotBestFit(weights):       # 画出数据集和logistic回归最佳拟合直线

    dataMat, labelMat =loadDataSet()
    dataArr=array(dataMat)
    n=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=arange(-3.0,3.0,0.1)      # 创建等差数组
    y=(-weights[0]-weights[1]*x)/weights[2]
    ax.plot(x,y)
    plt.xlabel("X1")
    plt.ylabel("X2")
    plt.show()

在这里插入图片描述

九、随机梯度上升算法及改进

def stocGradAscent(dataMatrix,classLabels):     # 随机梯度上升算法
    m,n=shape(dataMatrix)
    alpha=0.01
    weights=ones(n)
    for i in range(m):
        h=sigmoid(sum(dataMatrix[i]*weights))
        error=classLabels[i]-h
        weights=weights+alpha*error*dataMatrix[i]
    return weights

在这里插入图片描述

可以发现gradAscent函数和stocGradAscent函数的梯度上升算法很相似，但有区别。前者的error计算出的是向量，后者是数值，并且后者没有矩阵转换的过程。
拟合出来的直线也并不像梯度上升的效果那样好，这归咎于迭代次数的差距。

改进的随机梯度上升算法

def stocGradAscent1(dataMatrix,classLabels,numIter=1500):   # 改进的随机梯度上升算法
    m,n=shape(dataMatrix)
    weights=ones(n)
    for j in range(numIter):
        dataIndex=list(range(m))
        for i in range(m):
            alpha=4/(1.0+j+i)+0.01
            randIndex=int(random.uniform(0, len(dataIndex)))
            h=sigmoid(sum(dataMatrix[randIndex]*weights))
            error=classLabels[randIndex]-h
            weights=weights+alpha*error*dataMatrix[randIndex]
            del(dataIndex[randIndex])
    return weights

在这里插入图片描述

改进后的算法将每次迭代的步长设置为随着迭代次数不断减小。
这里通过随机选取样本来更新回归系数。这种方法将减少周期性的波动。

十、logistic回归预测病马死亡率

训练数据：在这里插入图片描述
测试数据：

分类代码：

def classifyVector(inX,weights):    # 分类0和1
    prob=sigmoid(sum(inX*weights))
    if prob >0.5:
        return 1.0
    else:
        return 0.0

def colicTest():        # logistic分类测试函数
    frTrain=open("horseColicTraining.txt")
    frTest=open("horseColicTest.txt")
    trainingSet=[]
    trainingLabels=[]
    for line in frTrain.readlines():
        currLine=line.strip().split("\t")
        lineArr=[]
        for i in range(21):
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr)
        trainingLabels.append(float(currLine[21]))
    trainWeights=stocGradAscent1(array(trainingSet),trainingLabels,500)
    errorCount=0
    numTestVec=0.0
    for line in frTest.readlines():
        numTestVec+=1.0
        currLine=line.strip().split("\t")
        lineArr=[]
        for i in range(21):
            lineArr.append(float(currLine[i]))
        if int(classifyVector(array(lineArr),trainWeights)!=int(currLine[21])):
            errorCount+=1
    errorRate=float(errorCount)/numTestVec
    print("the error rate of this test is: %f" %errorRate)
    return errorRate

def multiTest():
    numTests=10
    errorSum=0.0
    for k in range(numTests):
        errorSum+=colicTest()
    print("after %d iterations the average error rate is: %f" %(numTests,errorSum/float(numTests)))