logistic回归和Python实现

一、Logistic Regression的基本内容

通过学习了台湾的林教授和Stanford的课程后发现，他们两个人的基本思路虽然一致，但是具体做法有所差异，下面简单介绍一下两种实现方式。
1、台湾的林教授方法 
     使用Logistic回归的思路是，我们不希望只给输出一个确定的{-1，+1}的结果，而是想给出一个概率值，这时就可以用到Logistic回归，得到概率的输出值在0~1之间。 
   但是实际情况是，我们不能知道准确的概率值，只能知道输出的结果。比如以心脏病复发为例。绿色和红色分别表示理想情况下和实际的结果。 

 可以将实际的情况看成是含有噪声的理想训练数据。和以前一样，同样用输入加权的形式得到一个分数，并将此分数通过称为Logistic function或Sigmoid function的函数转换为一个0~1之间的值，在这里用θ（s）表示。s越大患病风险越高。假设函数h为。 

     和平时一样，方法是最小化Logistic回归的错误Ein。而Logistic回归的目标函数为

     假设存在一个数据集 D：，则产生这个数据集的概率为各个输入样本对应概率的连乘。通过把目标函数带入可以得到带有f（x）的形式，而得到f（x）也就是根本目的。

    但实际上f（x）是未知的，已知的只有假设函数h。是否可以用假设函数代替f呢？意味着假设函数h产生同样的数据集样本D的可能性有多大。这在数学上称为似然（likelihood）。（注意：这里和Stanford后面讲的就不是太一样了，个人感觉这种思路比较好理解。） 替代公式变成了

      如果假设函数h和未知函数f很接近（err很小），那么h产生数据样本D的可能性与f产生同样数据D的可能性很接近。由于f（x）产生了样本D，那么可以认为f产生该数据样本的几率很大。因此推断出最好的假设函数g，就是似然性最大的假设函数h。即。

         通过h的性质，可以将似然函数写为：，P（x）对似然的大小没有影响，只与函数h对每个样本的连乘有关。并将连乘通过取ln变成加和形式，为了得到Ein通过加负号，变成求最小值。

         这里，定义了交叉熵误差（cross entropy error）来衡量训练误差

       我们得到了训练误差的具体形式，由于这个训练误差函数是可微并且是凸函数，所以依照之前的思路，对这个函数求梯度。对每一个权值向量wi求偏微分，最后合并。

        这时分为两种情况，一种是线性可分时权值θ（.）为0，即y与wTX同号。排除这种特殊情况，当加权求和为0时，不能用类似解线性回归时的闭合求解方式，这时的计算就要用迭代求解的方式了。类似梯度上升、梯度下降和牛顿方法都是用于这种加权和的形式。

        按照类似的思路，一步一步的去修正w，使得Ein的结果越来越小。在logistic回归中，Ein的式子是平滑的。现在我们可以将这个Ein想象成一个山谷，我们要到达山谷的最低点，就是要沿着当前的梯度最大的方向每次迈出一小步，直到到达谷底，使得Ein最小，得到最佳解w。迭代函数为（怎么选择的？）。增加一个η表示步伐，v表示方向，则迭代函数可以变为。

         最终目标是最小化Ein（Wt+η'v）。要确定一个w使得Ein最小，不可能一步到位，指导思路还是由繁化简，用线性近似(linear approximation)的方式来解决问题，我们用多维度的泰勒展开公式来近似Ein，只要给定一个小的η，就可以近似这个Ein。

这里Ein(wt)和η都是已知的，Ein的梯度表示了下降的方向，也可以求出来，唯一要考虑的就是最好的向量v该如何选择。于是，对于v和Ein的梯度这两个向量，使得其值是最小的方法就是让v向量的方向和Ein的梯度的方向相反，这样使得两个向量的内积是最小的，另外由于v的模是1，还需要有个归一化的步骤。

为了找到一个适合的η，防止η太小时下降慢，太大时不稳定。让η随着梯度的变化变化（与梯度大小成正比）。 所以公式变为这里的η为一个新的值，被称为固定的学习速率。

总结：logistic回归求解最优化解的步骤为：

2、Stanford机器学习方法

       和上面的方法主要的区别从likelihood开始。上面方法主要是想让h能够最好的拟合f（x），使产生相同数据集的几率最大。Stanford的选择为对应参数θ（上面的方法参数名为w）的似然函数。但实际上是一样的，因为上面的方法在求解最大值得过程中略去了固定值P（x），得到的似然函数其实和下面的一样。但是上面的方法用的方式不是像下面这种方法将似然函数分开写。

下面是具体方法：

接着为了求解方便，取了对数。而目的就是最大化似然函数。

最大化似然函数有两种方法，第一种是梯度上升法，一种是牛顿方法。

    （1）梯度上升法

              ，不断更新参数。而梯度求解方法如下：

这里面用了一个公式，对于sigmoid function g(x),

所以，最终的梯度上升规则为

然后repeat所有的j。 但是其实不需要for所有的j，可以用矩阵来解决。（下面的做个参考就好，实际代码中一般用矩阵虽然全面，但是相应的速度也会很慢，代码中会用随机梯度上升法来解决这一问题）。

      （2）其他方法

              优化算法除了梯度下降算法外，还包括：

• Conjugate gradient method(共轭梯度法)
• Quasi-Newton method(拟牛顿法)
• BFGS method
• L-BFGS(Limited-memory BFGS)

             这些方法的有点是不用计算步伐就可以。

        下面简单介绍一下牛顿方法。牛顿法的目的是找到f（x）= 0的点，方法用的是泰勒展开不断的找到△x不断的迭代更新。公式为：。下图为不断迭代的过程。

对应于我们的问题，为了求解最大化的似然函数，我们可以令，得到迭代函数：

但在我们的logistic回归中，θ是向量，要把公式对应到高维的情况下，这时我们得到了：

其中H叫做Hessian。缺点是当特征数量很多（几千）计算的代价是很大的。

3、总结

       两种方法有各自的优缺点，第一种方法计算相对麻烦一点，但是得到的迭代函数的学习率却是随着梯度自适应变化的。第二种的计算方法比较直观，但是学习率不好处理，只能靠经验或测试观察。

      其实还有一种方法是直接用代价函数（cost function）做的。

还有更多的内容，如防止过拟合和欠拟合，多分类等。

4、python实现

      代码主要参考了《机器学习实战》这本书，还有http://blog.csdn.net/zouxy09/article/details/20319673所写的程序，他们的程序多少都有一点错误，这些我做了适当的修改，而且即使没有修改的地方也做了备注。程序主要用的是梯度上升法。

# -*- coding: utf-8 -*-

import random
from numpy import *
from math import *
import matplotlib.pyplot as plt

# 逻辑函数
def sigmoid(inX):
	return 1.0/(1+exp(-inX))
#加载数据
def loadDataSet(filename):
	numFeat = len(open(filename).readline().split('\t'))-1
	dataMat = []; labelMat = []
	fr = open(filename)
	for line in fr.readlines():
		lineArr = []
		curLine = line.strip().split('\t')
		for i in range(numFeat):
			lineArr.append(float(curLine[i]))
		dataMat.append([1]+lineArr)  #这里是为了使 x0 等于 1
		labelMat.append(float(curLine[-1]))
	return mat(dataMat),mat(labelMat).T

#Logistic 回归
def logisticRegression(dataMat,labelMat,opts):
	numSamples,numFeature = shape(dataMat)
	alpha = opts['alpha'] ; maxIter = opts['maxIter']
	weight = ones((numFeature,1))
	for j in range(maxIter):
		#选择方法共三种：梯度下降，随机梯度下降，优化的随机梯度下降
		if opts['optimizeType'] == 'gradDescent':# gradient descent algorilthm
			h = sigmoid(dataMat * weight)  #注意这里会报错，因为exp（）的自变量只能是数，不能是矩阵或数组
			error = labelMat - h
			weight = weight + alpha * dataMat.T * error
		elif opts['optimizeType'] == 'stocGradDescent': # stochastic gradient descent
			for i in range(numSamples):
				h = sigmoid(sum(dataMat[i]*weight))
				error = labelMat[i] - h
				weight = weight + alpha * dataMat[i].T* error
		elif opts['optimizeType'] == 'smoothStocGradDescent':# smooth stochastic gradient descent
			dataIndex = range(numSamples)
			for i in range(numSamples):
				randIndex = int(random.uniform(0,len(dataIndex)))#产生0到dataIndex之间的随机数
				alpha = 4.0/(1.0+j+i)+0.01 #j是迭代次数，i是迭代时，第i个选出的样本
				h = sigmoid(sum(dataMat[randIndex]*weight))
				error = labelMat[randIndex] - h
				weight = weight + alpha * dataMat[randIndex].T* error
				del(dataIndex[randIndex])
		else:
			raise NameError('Not support optimize method type!')
	return weight

#准确度计算
def testLogRegres(weights, test_x, test_y):
	numSamples, numFeatures = shape(test_x)
	matchCount = 0
	for i in xrange(numSamples):
		predict = sigmoid(test_x[i, :] * weights)[0, 0] > 0.5
		if predict == bool(test_y[i, 0]):
			matchCount += 1
	accuracy = float(matchCount) / numSamples
	return accuracy


#画出决策边界
def showLogRegres(weight,dataMat,labelMat):
	#输入的x和y为矩阵
	numSamples, numFeatures = shape(dataMat)
	weights = weight.getA() #将矩阵转为数组，经测试应该是对应于[[0],[1],[2]]这种的才有效
	if numFeatures != 3:
		print "Sorry! I can not draw because the dimension of your data is not 2!"
		return 1
		#xrange 与range相比返回的是一个xrange对象，而不是一个数列。xrange的对象是单个元素，可以节省空间。本质为一个生成器
		#在python3中，将xrange命名为range，而取消了range
	for i in xrange(numSamples):
		if int(labelMat[i]) == 0:
			plt.plot(dataMat[i,1],dataMat[i,2],'or')
		elif int(labelMat[i]) == 1:
			plt.plot(dataMat[i,1],dataMat[i,2],'ob')

	#画分类线
	min_x = min(dataMat[:,1])[0,0] # 只有用[0,0]才能得到一个数值
	max_x = max(dataMat[:, 1])[0,0]
	y_min_x = float(-weights[0] - weights[1] * min_x) / weights[2] # 这是令sigmoid函数的输入为零，得到分界处的值
	y_max_x = float(-weights[0] - weights[1] * max_x) / weights[2]
	plt.plot([min_x, max_x], [y_min_x, y_max_x], '-g')
	plt.xlabel = ('X1'); plt.ylabel = ('X2')
	plt.show()

其他的方法会在后面介绍。

参考资料:http://blog.csdn.net/abcjennifer/article/details/7716281（讲的代价函数的方法）

               http://blog.chinaunix.net/uid-20761674-id-4336428.html

                http://www.open-open.com/lib/view/open1421575796468.html