1 # -*- coding: utf-8 -*-
  2 
  3 ‘‘‘
  4 function: 实现Logistic回归，拟合直线，对数据进行分类；
  5             利用梯度上升，随机梯度上升，改进的随机梯度上升，牛顿法分别对损失函数优化；
  6             这里没有给出最后测试分类的函数；
  7 date: 2017.8.12
  8 ‘‘‘
  9 
 10 from numpy import *
 11 
 12 #从文件加载处理数据
 13 def loadDataSet():
 14     dataMat = []
 15     labelMat = []
 16     fr = open(‘testSet.txt‘)
 17     for line in fr.readlines():
 18         lineArr = line.strip().split()
 19         dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
 20         labelMat.append(int(lineArr[2]))
 21     return dataMat, labelMat
 22 
 23 #sigmoid function X: w1*x1+w2*x2+...+wn*xn
 24 def sigmoid(x):
 25     return 1 / (1 + exp(-x))
 26 
 27 #梯度上升求函数的最大值时取得的权重
 28 def gradAscent(dataMatIn, classLabels):
 29     dataMatrix = mat(dataMatIn)
 30     labelMat = mat(classLabels).transpose()   #将m维行向量转制为m维列向量
 31     m,n = shape(dataMatrix)
 32     alpha = 0.001 #设置梯度上升的步长
 33     maxCycles = 500 #最大迭代次数 
 34     weights = ones((n,1))  #weights就是theta，n维列向量,二维数组
 35     for i in range(maxCycles):
 36         h = sigmoid(dataMatrix*weights)  #计算所有数据的分类概率，h是m维向量,这里实际上进行了300次乘法运算
 37         error = labelMat - h #计算(y-h(x))：误差
 38         weights = weights + alpha * dataMatrix.transpose()*error #对所有weights同时更新
 39     print(‘shape of weights‘, shape(weights))
 40     return weights
 41 
 42 #随机上升上升
 43 def stocGradAscent0(dataMatrix, classLabels):
 44     m,n = shape(dataMatrix)
 45     alpha = 0.01
 46     weights = ones(n, float) #n 维行向量
 47     for i in range(m):
 48         h = sigmoid(sum(dataMatrix[i] * weights))
 49         error = classLabels[i] - h
 50         weights = weights + alpha * error * dataMatrix[int(i)]
 51     return weights
 52 
 53 #改进的随机上升上升
 54 def stocGradAscent1(dataMatrix, classLabels, numIter = 550):
 55     m,n = shape(dataMatrix)
 56     weights = ones(n)
 57     dataIndex = []
 58     for j in range(numIter):
 59         for k in range(m):
 60             dataIndex.append(k)
 61         for i in range(m):
 62             alpha = 4 / (i + j + 1.0) + 0.01 #aloha随着迭代次数增多不断减小但是不为0，为了解决局部波动
 63             randIndex = int(random.uniform(0, len(dataIndex))) #随机选取一个数据进行更新参数
 64             h = sigmoid(dataMatrix[randIndex] * weights)
 65             error = classLabels[randIndex] - h
 66             weights = weights + alpha * error * dataMatrix[randIndex]
 67             del(dataIndex[randIndex]) #这里会不会对本来的dataMatrix有影响,不会，因为没有操作矩阵
 68     return weights
 69 
 70 ‘‘‘
 71 functuion: 牛顿法更新theta,计算海森矩阵
 72 note: 这里一定要对 XT*X 结果求逆，不然分类效果无法直视
 73 ‘‘‘
 74 def computeHessianMatrix(dataMatrix):
 75     hessianMatrix = mat(dataMatrix).transpose().dot(mat(dataMatrix))
 76     return hessianMatrix.I #矩阵求逆
 77 
 78 ‘‘‘
 79 function: 牛顿法更新theta
 80 ‘‘‘
 81 def newtonMethod(dataMat, labelMat, numIter = 10):
 82     m,n = shape(dataMat)
 83     dataMatrix = mat(dataMat)
 84     labelMat = mat(labelMat).transpose()
 85     weights = ones((n,1)) #参数向量
 86     hessianMatrix = computeHessianMatrix(dataMatrix)
 87     print(‘shape of hessian‘, shape(hessianMatrix))
 88     for k in range(numIter):
 89         h = sigmoid(dataMatrix * weights)
 90         error = (labelMat - h)
 91         weights = weights - (dataMatrix*hessianMatrix).transpose() * error 
 92     return weights
 93 
 94 #画出决策边界
 95 def plotBestFit(weights):
 96     import matplotlib.pyplot as plt
 97     dataMat, labelMat = loadDataSet() #或取数据和标签
 98     dataArr = array(dataMat) #将二维列表转换为数组
 99     n = shape(dataArr)[0] #行数，代表数据的个数
100 
101     xcord1 = []; ycord1 = []
102     xcord2 = []; ycord2 = []
103     for i in range(n):
104         if int(labelMat[i]) == 1:
105             xcord1.append(dataArr[i,1]); ycord1.append(dataArr[i,2])
106         else:
107             xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2])
108 
109     fig = plt.figure()
110     ax = fig.add_subplot(111)
111     ax.scatter(xcord1, ycord1, s = 30, c = ‘red‘, marker = ‘s‘)
112     ax.scatter(xcord2, ycord2, s = 30, c = ‘green‘)
113     x = arange(-3.0, 3.0, 0.1)
114     print(len(x))
115     print(len(weights))
116     print(shape(weights))
117     y = (-weights[0] -  weights[1]*x) / weights[2] #把x2堪称y,x0=1,所以就是x1,x2之间的函数
118     print(len(y))
119     ax.plot(x,y) 
120     plt.xlabel(‘X1‘); plt.ylabel(‘X2‘)
121     plt.show() #显示图像
122 
123 ‘‘‘
124 function: 对测试数据进行测试
125 input: testDate 是一个向量，或者多个向量
126 ‘‘‘
127 def logisticTest(weights, testData):
128     m,n = shape(testData)
129     typeLabel = []
130     for i in range(m): 
131         result = sigmoid(sum(testData[i] * weights)) #得到一个册数数据的概率
132         if result > 0.5:
133             typeLabel.append(1)
134         else:
135             typeLabel.append(0)
136     return typeLabel
137 
138 ‘‘‘
139 function: 计算分类的正确率
140 input: calssLables 是测试数据的类别向量
141 ‘‘‘
142 def getCorrectRate(classLabels, testLabel):
143     correctRate = 0.0
144     numOfRight = 0
145     for i in range(len(testLabel)):
146         if classLabels[i] == testLabel[i]:
147             numOfRight += 1
148     correctRate = numOfRight / float(len(testLabel))
149     return correctRate
150 
151 #测试算法
152 dataMat, labelMat = loadDataSet()
153 print(‘shape of datamet, ‘, shape(dataMat))
154 weights1 = gradAscent(dataMat, labelMat)
155 weights2 = stocGradAscent0(array(dataMat), labelMat)
156 weights3 = stocGradAscent1(dataMat, labelMat)
157 
158 weights4 = newtonMethod(dataMat, labelMat, 2000) #牛顿法
159 
160 
161 #测试正确率
162 typeLabel = logisticTest(weights4, [[1,1,0],[1,1,1]])
163 print(typeLabel)
164 correctRate = getCorrectRate([1,0],typeLabel)
165 print(‘正确率： ‘ ,str(correctRate))
166 
167 #画图
168 #plotBestFit(array(weights1)) #这里因为普通梯度上升返回的是一个矩阵，所以要转换为向量
169 #plotBestFit(weights2)
170 #plotBestFit(weights3)
171 plotBestFit(array(weights4))