逻辑斯谛回归
关于逻辑斯谛回归,这篇文章http://blog.csdn.net/zouxy09/article/details/20319673 讲的很好;Andrew Ng的机器学习公开课也很不错(中文笔记也很好http://blog.csdn.net/stdcoutzyx );还有《机器学习实战》,都是不错资料。
在逻辑斯谛回归中,因为使用梯度上升(gradient ascent)收敛较慢,固本文采用牛顿法(Newton’s Method)进行参数求解,试验发现通常迭代10次左右就可达到收敛,而梯度上升法则需要迭代上百甚至上千次,当然实际的迭代次数也要视实际数据而定。
牛顿法
牛顿法与梯度下降法的功能一样,都是最优化的常用方法。
对于一个函数,如果要求函数值为0时的值,如图所示:
先随机选一个点,然后求出该点的切线,即导数,延长切线与横轴相交,以相交时的的值作为下一次迭代的值,更新规则如下
对于逻辑斯谛回归,需要求的是似然函数L(θ)的最大值,当L(θ)的导数L’(θ)为0时即为L(θ)的最大值,即求L’(θ)=0的参数,则可使用牛顿法进行求解,此时参数更新规则为
使用牛顿法的另一个好处是不需要像梯度法一样指定学习率(即步长)。但是牛顿法需要对二阶导(Hessian矩阵)进行求逆,不过随着拟牛顿法(BFGS)以及限域拟牛顿法(LBFGS)的提出,大大减少了求逆的计算量,不过在本文还是使用牛顿法进行参数求解。
牛顿法求解逻辑斯谛回归参数
迭代中需要进行的主要步骤包括如下:
(1) 初始化参数θ
(2) 获取数据x
(3) 对数据进行预测h
(4) 得到对数似然函数L(θ)
(5) 根据L(θ)计算梯度g
(6) 根据L(θ)计算Hessian矩阵H
(7) 更新参数θ
具体计算
(1) 初始化θ=(b, θ(1), θ(2), …, θ(n))T,初始时θ=(0, 0, 0, …, 0)T
(2) 获取x=(1, x(1), x(2), …, x(n))T,其中x(1)到x(n)为数据的n个特征的值
(3) 使用逻辑斯谛函数计算预测值h
(4) 得到似然函数
对数似然函数为
(5) 根据L(θ)计算梯度,一个维度的梯度计算公式为
n+1维的梯度向量为
(6) 根据L(θ)计算Hessian矩阵H,其中j行k列的值为
(n+1)*(n+1)的Hessian矩阵为
(7) 更新参数θ,θ = θ – H-1 g
本文使用的迭代规则是求出m个样本的梯度g以及Hessian矩阵H,对m个g和H求平均,然后更新θ,伪代码:
while(iterNum) //迭代次数 { <span style="white-space:pre"> </span>for i range(1,m) //m个样本 <span style="white-space:pre"> </span>{
<span style="white-space:pre"> </span>(2) (3) (4) (5) (6) //具体计算的步骤标号 } g = (1/m)g H = (1/m)H θ = θ – H-1g iterNum-- }
测试
(1)使用Andrew Ng教授的成绩数据进行测试,对原始数据进行了一些预处理http://openclassroom.stanford.edu/MainFolder/courses/MachineLearning/exercises/ex4materials/ex4Data.zip
测试效果如下:
(2)使用马疝气病数据http://archive.ics.uci.edu/ml/datasets/Horse+Colic 进行了测试,迭代7次后收敛,分类错误率0.253731,使用的数据是《机器学习实战》中预处理后的数据。
tie dai ma
依赖包括numpy库和matplotlib库。
代码比较乱,因为开始时并没有使用numpy库,所以向量的点积、相加等计算都是自己实现的。后来使用牛顿法需要求矩阵逆,自己偷懒不想写,所以用了numpy。中间还有一些计算本来可以使用numpy中的方法直接进行计算,但可能是对矩阵的初始化方法不对,始终无法用numpy计算矩阵乘法,无奈只能自己用笨方法进行计算。代码较乱,做好瞎眼准备,ready! go!
# -*- coding: gb18030 -*- __author__ = 'jinyu' from numpy import * import math ## # 读取训练数据 # # def loadDataSet(dataFileName): dataMat = [] labelMat = [] fr = open(dataFileName) for line in fr: lineArr = line.strip().split() dataMat.append([1.0, 10.0*float(lineArr[0]), 10.0*float(lineArr[1])]) labelMat.append(int(lineArr[2])) return dataMat, labelMat ## # sigmod函数 # # def sigmoid(x): return 1.0 / (1+math.exp(-x)) ## # 梯度上升法 # # def gradientAscent(dataMat, labelMat, alpha): m = len(dataMat) #训练集个数 n = len(dataMat[0]) #数据特征纬度 theta = [0] * n iter = 1000 while(iter): for i in range(m): hypothesis = sigmoid(computeDotProduct(dataMat[i], theta)) error = labelMat[i] - hypothesis gradient = computeTimesVect(dataMat[i], error) theta = computeVectPlus(theta, computeTimesVect(gradient, alpha)) iter -= 1 return theta ## # 牛顿法 # # def newtonMethod(dataMat, labelMat, iterNum=10): m = len(dataMat) #训练集个数 n = len(dataMat[0]) #数据特征纬度 theta = [0.0] * n while(iterNum): gradientSum = [0.0] * n hessianMatSum = [[0.0] * n]*n for i in range(m): try: hypothesis = sigmoid(computeDotProduct(dataMat[i], theta)) except: continue error = labelMat[i] - hypothesis gradient = computeTimesVect(dataMat[i], error/m) gradientSum = computeVectPlus(gradientSum, gradient) hessian = computeHessianMatrix(dataMat[i], hypothesis/m) for j in range(n): hessianMatSum[j] = computeVectPlus(hessianMatSum[j], hessian[j]) #计算hessian矩阵的逆矩阵有可能异常,如果捕获异常则忽略此轮迭代 try: hessianMatInv = mat(hessianMatSum).I.tolist() except: continue for k in range(n): theta[k] -= computeDotProduct(hessianMatInv[k], gradientSum) iterNum -= 1 return theta ## # 计算hessian矩阵 # # def computeHessianMatrix(data, hypothesis): hessianMatrix = [] n = len(data) for i in range(n): row = [] for j in range(n): row.append(-data[i]*data[j]*(1-hypothesis)*hypothesis) hessianMatrix.append(row) return hessianMatrix ## # 计算两个向量的点积 # # def computeDotProduct(a, b): if len(a) != len(b): return False n = len(a) dotProduct = 0 for i in range(n): dotProduct += a[i] * b[i] return dotProduct ## # 计算两个向量的和 # # def computeVectPlus(a, b): if len(a) != len(b): return False n = len(a) sum = [] for i in range(n): sum.append(a[i]+b[i]) return sum ## # 计算某个向量的n倍 # # def computeTimesVect(vect, n): nTimesVect = [] for i in range(len(vect)): nTimesVect.append(n * vect[i]) return nTimesVect def plotBestFit(dataMat, labelMat, weights): import matplotlib.pyplot as plt 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(0.0, 70.0, 0.1) y = (-weights[0]-weights[1]*x)/weights[2] ax.plot(x, y) plt.xlabel('X1'); plt.ylabel('X2'); plt.show() def classifyVector(inX, weights): prob = sigmoid(sum(inX*weights)) if prob > 0.5: return 1.0 else: return 0.0 def colicTest(tainFileName, testFileName): frTrain = open(tainFileName); frTest = open(testFileName) 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 = newtonMethod(trainingSet, trainingLabels, 10) 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 dataMat, labelMat = loadDataSet("ex4x.dat") theta = newtonMethod(dataMat, labelMat, 10) print theta plotBestFit(dataMat, labelMat, array(theta)) #colicTest('horseColicTraining.txt', 'horseColicTest.txt')
代码和数据
代码和数据可以在https://github.com/0zone/LogisticRegression 下载