多变量logistic回归glm函数结果 logistic回归多分类变量

优点：计算代价小，容易实现和理解；

缺点：容易欠你和，分类精度可能不高

适用：数字型和标称型数据

在只需要两个分类的情况下，我们希望在某个临界分类线下，一侧为正类另一侧为负类，这非常像阶跃函数的性质。但是阶跃函数在非常临近分类线上的数据很容易出现误非类，所以平常往往用与其性质非常类似的sigmoid函数代替其进行分类。其公式如下：

对应图形如下：

从上图可以看出，sigmoid函数的分界线在y=0.5上，小于0.5则x<0，为负类；大于0.5则x>0，为正类。我们可以利用这个性质来对数据进行分类。设定输入为：

其中，b为常数不需要计算，从而可知，只需要计算最优的参数

，就可以对数据进行分类。

这里我们用梯度上升算法来计算最优参数。

梯度上升算法：

1）一个简单二次函数的例子

假设有二次函数如下：

我们可以根据二次函数的性质，求其导数从而获得其最大值：

令其导数等于0，可以知道，这个函数最大值为7。

同样，用梯度上升算法，我们可以设置一个初始值，然后根据学习速率，一点点的去逼近最大值，从而获取到一个较优解：

其中，α为学习速率，决定每次”爬坡“的速度。

python代码：

import matplotlib.pyplot as plt
import  numpy as np
"""
梯度上升算法
求函数f(x)=-2x^2+4x+5的最大值
"""
# 函数导数
def Derivatives(old):
    return -4*old+4

#梯度上升算法
def GradientAscent():
    xcord = []  #保存x坐标
    ycord = []  #保存y坐标
    xold = -1
    xnew = 0    #梯度算法的初始值
    alpha = 0.01 # 学习速率（一般根据经验和具体需要设置）
    maxCycle = 100 #学习次数
    for k in range(maxCycle):
        xold = xnew
        xnew = xold + alpha *Derivatives(xold) #梯度上升公式
        xcord.append(xnew)
        ycord.append(-2*np.square(xnew) + 4*xnew+5)
    return xcord,ycord #返回每一次学习完后的x,y值

#画图
def plotFig(xcord,ycord):
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xcord, ycord, s=4, c='red', marker='s', alpha=.5) #绘制每次学习后x,y值
    x = np.arange(0,2,0.01)
    y = -2*np.square(x) + 4*x+5
    ax.plot(x,y)            #绘制-2x^2+4x+5 函数
    plt.title('BestFit')  # 绘制title
    plt.xlabel('X')
    plt.ylabel('Y')  # 绘制label
    plt.show()

if __name__== '__main__':
    xcord, ycord = GradientAscent()
    plotFig(xcord,ycord)

得到下图：

可以看到经过100次训练后，获取的最优值已经非常接近真实的最大值7了。

2）sigmoid函数例子

根据sigmoid函数性质假定：

我们可以将其合并，合并的函数称为代价函数：

为了方便计算，利用取对数不改变函数单调性的性质，取代价函数的对数：

我们对其求导，以用于计算最优值：

从而我们得到梯度上升算法：

由于，sigmoid函数分类可以有很多数据，推广到多数据：

数据为testSet.txt，如下，共100条训练数据，前两列对应X1，X2的值，第三列是分类0或者1：

-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

Python代码：

import matplotlib.pyplot as plt
import numpy as np

"""
说明:利用sigmoid函数分类

"""
def loadDataSet():
    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])])        #添加数据，前面加个1是为了方便矩阵运算
        labelMat.append(int(lineArr[2]))                                #添加分类
    fr.close()                                                            
    return dataMat, labelMat

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

# 梯度上升算法
def gradAscent(dataMatIn, classLabels):
    dataMatrix = np.mat(dataMatIn)                                        #转换成numpy的mat，x0=1,x1=第一列数据，x2=第二列数据
    labelMat = np.mat(classLabels).transpose()                            #转换成numpy的mat,并进行转置
    m, n = np.shape(dataMatrix)                                            #返回dataMatrix的大小。m为行数,n为列数。
    alpha = 0.001                                                        #移动步长,也就是学习速率,控制更新的幅度。
    maxCycles = 500                                                        #最大迭代次数
    weights = np.ones((n,1))                                              #设置w^T的初始值为（1,1,1）
    for k in range(maxCycles):
        h = sigmoid(dataMatrix * weights)                                #梯度上升矢量化公式
        error = labelMat - h
        weights = weights + alpha * dataMatrix.transpose() * error
    return weights.getA()                                                #将矩阵转换为数组，返回权重数组

#画图
def plotBestFit(weights):
    dataMat, labelMat = loadDataSet()                                    #加载数据集
    dataArr = np.array(dataMat)                                            #转换成numpy的array数组
    n = np.shape(dataMat)[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])    #1为正样本
        else:
            xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2])    #0为负样本
    fig = plt.figure()
    ax = fig.add_subplot(111)                                            
    ax.scatter(xcord1, ycord1, s = 20, c = 'red', marker = 's',alpha=.5)#绘制正样本
    ax.scatter(xcord2, ycord2, s = 20, c = 'green',alpha=.5)            #绘制负样本
    x1 = np.arange(-3.0, 3.0, 0.001)                                     #生成步长为0.001的X值
    x2 = (-weights[0] - weights[1] * x1) / weights[2]                     #根据w0*x0+w1*x1+w2*x2=0获得x2的值
    ax.plot(x1, x2)                                                      #根据x1和x2画出分类直线
    plt.title('BestFit')                                                #绘制title
    plt.xlabel('X1'); plt.ylabel('X2')                                    #绘制label
    plt.show()

if __name__ == '__main__':
    dataMat, labelMat = loadDataSet()
    weights = gradAscent(dataMat, labelMat)
    plotBestFit(weights)

得到下图：

3)随机梯度上升算法

梯度上升算法，每次更新回归系数，都需要对整个数据进行计算，在数据量大的情况下效率非常低，因而出现了随机梯度上升算法。我们更新回归系数不再计算整个数据。而是利用随机性，随机的选择一组数据继续计算，并且不断的减小学习速率，已尽可能找到最优解：

Python代码（其他代码同上）： 

# 随机梯度上升算法
def stocGradAscent(dataMatrix, classLabels, numIter=150):
    m,n = np.shape(dataMatrix)
    weights = np.ones(n)   #初始值w0,w1,w2都设为1
    for j in range(numIter):
        dataIndex = list(range(m))        #根据数据大小设置随机范围最大值
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.0001    #学习速率每次都下降
            randIndex = int(np.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

if __name__ == '__main__':
    dataMat, labelMat = loadDataSet()
    weights = stocGradAscent(np.array(dataMat), labelMat)
    plotBestFit(weights)

得到下图：

4）用Logistic回归预测病马死亡率

数据准备：

现实环境中，可能会存在多个特征值，也会存在特征值丢失的情况。在这种情况下，我们就需要对特征值进行补全，常见的做法包括：

1，用特征值的均值补全；

2，用特殊值补全，如-1；

3，用相似样本特征值补全；

4，用其他机器学习方法预测特征值；

对于Logistic回归来说，我们可以采用第二种方式，把缺失的特征值置为0,。因为特征值为0，根据计算回归公式weights = weights + alpha * error*dataMatrix[randIndex]的值不会改变，即weights = weights。另外，sigmoid（0）=0.5也不会对结果产生影响。

python代码：

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

def colicTest():
    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 = stocGradAscent(np.array(trainingSet), trainingLabels, 1000) #通过随机梯度上升算法确定回归系数
    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(np.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)))