Python logRegres

from numpy import *#数据处理函数，读取文件testSet.txt的内容#前两列分别为x1和x2值，第3列为数据的类别def loadDataSet():    #存放训练数据    dataMat = [];    #每条数据的类别    labelMat = []    #打开文件读取训练样本    fr = open('testSet.txt')    #按行读取文件内容    for line in fr.readlines():        lineArr = line.strip().split()                #为计算方便，将x0设置为1        #线性回归为h(x)=W0*1+W1*X1+W2*X2        #(W0,W1,W2)*(1,X1,X2),(W0,W1,W2)为所求回归系数                #以下为(1,X1,X2)的值        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])                #每条数据类别存入列别向量labelMat中        labelMat.append(int(lineArr[2]))        #返回训练数据与类别标签    return dataMat,labelMat#sigmoid函数，输入值为inXdef sigmoid(inX):    return 1.0/(1+exp(-inX))#输入为训练数据，类别标签def gradAscent(dataMatIn, classLabels):    #转为NumPy可识别的矩阵    dataMatrix = mat(dataMatIn)                 #为了便于计算，classLabels为行向量转为列向量    labelMat = mat(classLabels).transpose()     #获取输入数据的条数m，特征数n    m,n = shape(dataMatrix)    #设定迭代的步长alpha    alpha = 0.001    #设置循环次数500次，即训练次数，人为给定    maxCycles = 500    #权值初始化为1，后面根据样本数据调整    #训练结束得到最优权值    #weights为n行，1维。为列向量。    weights = ones((n,1))    #循环maxCycles次，    #每次根据模型输出结果与真实值的误差，调整权值.    for k in range(maxCycles):        #dataMatrix*weights矩阵的乘法。        #事实上包含300次的乘积        #h为模型给出的一个预测值        h = sigmoid(dataMatrix*weights)        #计算误差，每条记录真实值与预测值之差        error = (labelMat - h)        #权值调整(未知参数调整)        weights = weights + alpha * dataMatrix.transpose() *  error         #循环次数结束，返回回归系数    return weightsdef plotBestFit(weights):    import matplotlib.pyplot as plt    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()        #随机梯度上升法#输入：训练数据矩阵，类别标签#输出：权值wdef stocGradAscent0(dataMatrix, classLabels):    #获取训练数据条数m，维度n    m,n = shape(dataMatrix)    #步长    alpha = 0.01    #权值初始化为1，后面根据样本数据调整    weights = ones(n)       #遍历每一条数据    for i in range(m):        #h为当前样本的预测值，批处理梯度上升算法的h为所有样本的模型输出        #此处h为一个值，即一次只有一个样本更新        #dataMatrix[i]*weights也为当前样本行乘以权值weights        h = sigmoid(sum(dataMatrix[i]*weights))        #误差，此处error为一个值        error = classLabels[i] - h        #只选择当前样本进行权值更新        weights = weights + alpha * error * dataMatrix[i]        #返回权值    return weights    #输入：训练数据、类别标签、默认迭代次数#输出：权值def stocGradAscent1(dataMatrix, classLabels, numIter=150):    #获取训练数据条数m，列数n    m,n = shape(dataMatrix)    #权值初始化为1，后面根据样本数据调整    weights = ones(n)    #循环numIter次    for j in range(numIter):        #dataIndex用于随机选择样本来更新系数        dataIndex = range(m)        #遍历每一条数据        for i in range(m):            #步长不在为定值            #每次迭代时都调整alpha值，使其减小            alpha = 4/(1.0+j+i)+0.0001                #随机选取样本对系数(权值)进行更新            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#定义分类函数#输入：需要进行类别判断的数据inX，训练好的权值weights#输出类别0和1def classifyVector(inX, weights):    #使用sigmoid函数计算    prob = sigmoid(sum(inX*weights))        #计算结果大于0.5为第1类，小于0.5为第0类    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():        #原始数据按tab键分开，读取数据        currLine = line.strip().split('\t')        lineArr =[]        #样本数据共21个属性特征        #遍历每列数据，将改行数据转换为行向量，便于整体操作        for i in range(21):            lineArr.append(float(currLine[i]))        #每条数据处理完后加入训练数据集中        trainingSet.append(lineArr)        #最后一列为每条数据类别，加入类别标签矩阵        trainingLabels.append(float(currLine[21]))        #使用改进随机梯度上升(下降)法训练1000次        trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 1000)    #记录错误次数    errorCount = 0;    #记录测试次数    numTestVec = 0.0    #按行读取测试集数据    for line in frTest.readlines():        #每读一行测试次数加1        numTestVec += 1.0        #制表符分割出每行数据        currLine = line.strip().split('\t')        lineArr =[]        #遍历每条数据的每列的值        for i in range(21):            #转为行向量便于操作            lineArr.append(float(currLine[i]))            #预测的类别与真实类别不一致，则错误次数加1        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#一次性测试结果错误率并不准确，#可以进行多次测试求平均错误率。#multiTest()函数用于调用colicTest()函数#即调用10次，每次得到一个错误率，求10次的平均错误率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))if __name__=="__main__":    dataMatIn,classLabels=loadDataSet()        """批量梯度下降法"""    #weights=gradAscent(dataMatIn, classLabels)    #print weights    #plotBestFit(weights.getA())        """随机梯度下降法"""    #weights=stocGradAscent0(array(dataMatIn), classLabels)    #print weights    #plotBestFit(weights)            """改进的随机梯度下降法 """    weights= stocGradAscent1(array(dataMatIn), classLabels)    print weights    plotBestFit(weights)                  """实例 疝气病症预测病马是否存活"""    #multiTest()