参考：http://blog.csdn.net/dongtingzhizi/article/details/15962797
1.简述
在线性回归中，h函数的输出值为连续值，当需要进行归类时，输出的应该是离散值，如何将连续值转换成离散值？
如果分类结果只有两个，用1,0表示。我们希望有：函数1/(1+e^(-z))，这样就可以将函数结果限定在0~1之间。
Logistic Regression 有三个主要组成部分：回归、线性回归、Logsitic方程。
1）回归其实就是对已知公式的未知参数进行估计。
比如已知公式是y = a*x + b，未知参数是a和b。我们现在有很多真实的(x,y)数据（训练样本）， 回归就是利用这些数据对a和b的取值去自动估计。估计的方法大家可以简单的理解为，在给定训练样本点和已知的公式后，对于一个或多个未知参数，机器会自动 枚举参数的所有可能取值（对于多个参数要枚举它们的不同组合），直到找到那个最符合样本点分布的参数（或参数组合）。
注意，回归的前提是公式已知，否则回归无法进行。而现实生活中哪里有已知的公式啊（G=m*g 也是牛顿被苹果砸 了脑袋之后碰巧想出来的不是？哈哈），因此回归中的公式基本都是数据分析人员通过看大量数据后猜测的（其实大多数是拍脑袋想出来的，嗯...）。根据这些 公式的不同，回归分为线性回归和非线性回归。线性回归中公式都是“一次”的（一元一次方程，二元一次方程...），而非线性则可以有各种形式（N元N次方 程，log方程 等等）。具体的例子在线性回归中介绍吧。
2）线性回归
直接来一个最简单的一元变量 的例子：假设要找一个y和x之间的规律，其中x是鞋子价钱，y是鞋子的销售量。（为什么要找这个规律呢？这样的话可以帮助定价来赚更多的钱嘛，小学的应用 题经常做的呵呵）。已知一些往年的销售数据（x0,y0), (x1, y1), ... (xn, yn)做样本集, 并假设它们满足线性关 系：y = a*x + b （其中a,b的具体取值还不确定），线性回归即根据往年数据找出最佳的a, b取值，使 y = a * x + b 在所 有样本集上误差最小。
需要注意的是，这里线性回归能过获得好效果 的前提是y = a*x + b 至少从总体上是有道理的（因为我们认为鞋子越贵，卖的数量越少，越便宜卖的越多。另外鞋子质量、广告投入、客流量等都有 类似规律）；但并不是所有类型的变量都适合用线性回归，比如说x不是鞋子的价格，而是鞋子的尺码），那么无论回归出什么样的（a,b），错误率都会极高 （因为事实上尺码太大或尺码太小都会减少销量）。总之：如果我们的公式假设是错的，任何回归都得不到好结果。
3）Logistic方程
上面我们的sell是一个具体的实数值，然而很多情况下，我们需要回归产生一个类似概率值的0~1之间的数值（比如某一双鞋子今天能否卖出去？或者某一个广 告能否被用户点击? 我们希望得到这个数值来帮助决策鞋子上不上架，以及广告展不展示）。这个数值必须是0~1之间，但sell显然不满足这个区间要求。 于是引入了Logistic方程，来做归一化。这里再次说明，该数值并不是数学中定义的概率值。那么既然得到的并不是概率值，为什么我们还要费这个劲把数值归一化为0~1之间呢？归一化的好处在于数值具备可比性和收敛的边界，这样当你在其上继续运算时（比如你不仅仅是关心鞋子的销量，而是要对鞋子卖出的可 能、当地治安情况、当地运输成本 等多个要素之间加权求和，用综合的加和结果决策是否在此地开鞋店时），归一化能够保证此次得到的结果不会因为边界 太大 /太小导致 覆盖其他feature 或 被其他feature覆盖。（举个极端的例子，如果鞋子销量最低为100，但最好时能卖无限多个，而当地治安 状况是用0~1之间的数值表述的，如果两者直接求和治安状况就完全被忽略了）这是用logistic回归而非直接线性回归的主要原因。到了这里，也许你已 经开始意识到，没错，Logistic Regression 就是一个被logistic方程归一化后的线性回归，仅此而已。
2. 具体过程
2.1 预测函数
（对应是一个取值在0和1之间的S型曲线）
接下来需要确定数据划分的边界类型，对于图2和图3中的两种数据分布，显然图2需要一个线性的边界，而图3需要一个非线性的边界。接下来我们只讨论线性边界的情况
图2：线性边界的形式：
预测函数：
hθ(x)函数的值有特殊的含义，它表示结果取1的概率，因此对于输入x分类结果为类别1和类别0的概率分别为：
2.2 构造一个Cost函数（损失函数），该函数表示预测的输出（h）与训练数据类别（y）之间的偏差，可以是二者之间的差（h-y）或者是其他的形式。综合考虑所有训练数据的
“损失”，将Cost求和或者求平均，记为J(θ)函数，表示所有训练数据预测值与实际类别的偏差。
2.3 显然，J(θ)函数的值越小表示预测函数越准确（即h函数越准确），所以这一步需要做的是找到J(θ)函数的最小值。
找函数的最小值有不同的方法，Logistic Regression实现时用的是梯度下降法（Gradient Descent）
θ的更新过程：
将J(θ)函数待入推导可得：
因为式中α本来为一常量，所以1/m一般将省略，所以最终的θ更新过程为：
θ更新过程的vectorization： 1）求A=x.θ
2）求E=g(A)-y；
3）求θ:=θ-α.x‘.E,x‘表示矩阵x的转置。
3. python 实现
3.1 加载数据
#-*-coding:utf-8-*-‘
from numpy import *
def loadDataSet():
dataMat = []; labelMat = []
fr = open(‘D://machine_learn//Ch05//testSet.txt‘)
for line in fr.readlines():
lineArr = line.strip().split()
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
labelMat.append(int(lineArr[2]))
return dataMat,labelMat
data
3.2 sigmod函数
1 def sigmoid(inX): 2 return 1.0/(1+exp(-inX))
sigmoid
3.3 梯度上升法寻找最优参数
1 def gradAscent(dataMatIn, classLabels): 2 dataMatrix = mat(dataMatIn) #convert to NumPy matrix 3 labelMat = mat(classLabels).transpose() #convert to NumPy matrix 4 m,n = shape(dataMatrix) 5 alpha = 0.001 6 maxCycles = 500 7 weights = ones((n,1)) 8 for k in range(maxCycles): #heavy on matrix operations 9 h = sigmoid(dataMatrix*weights) #1) 10 error = (labelMat - h) #2) 11 weights = weights + alpha * dataMatrix.transpose()* error #3) 12 return weights
gradAscent
3.4 画图
1 def plotBestFit(weights): 2 import matplotlib.pyplot as plt 3 dataMat,labelMat=loadDataSet() 4 dataArr = array(dataMat) 5 n = shape(dataArr)[0] 6 xcord1 = []; ycord1 = [] 7 xcord2 = []; ycord2 = [] 8 for i in range(n): 9 if int(labelMat[i])== 1: 10 xcord1.append(dataArr[i,1]); ycord1.append(dataArr[i,2]) 11 else: 12 xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2]) 13 fig = plt.figure() 14 ax = fig.add_subplot(111) 15 ax.scatter(xcord1, ycord1, s=30, c=‘red‘, marker=‘s‘) 16 ax.scatter(xcord2, ycord2, s=30, c=‘green‘) 17 x = arange(-3.0, 3.0, 0.1) 18 y = (-weights[0]-weights[1]*x)/weights[2] 19 ax.plot(x, y) 20 plt.xlabel(‘X1‘); plt.ylabel(‘X2‘); 21 plt.show() 22 dataMat,labelMat=loadDataSet() 23 weights=gradAscent(dataMat, labelMat) 24 plotBestFit(weights.getA())
plot
注意：print weights >>>[[ 4.12414349]
[ 0.48007329]
[-0.6168482 ]]
print weights[0] #[[ 4.12414349]]
print weights.getA()[0] #[ 4.12414349]
3.4 随机梯度上升法：改进后的梯度上升法
梯度下降算法在每次更新回归系数的时候都需要遍历整个数据集（计算整个数据集的回归误差），该方法对小数据集尚可。但当遇到有数十亿样本和成千上万的特征 时，就有点力不从心了，它的计算复杂度太高。改进的方法是一次仅用一个样本点（的回归误差）来更新回归系数。这个方法叫随机梯度下降算法。由于可以在新的 样本到来的时候对分类器进行增量的更新（假设我们已经在数据库A上训练好一个分类器h了，那新来一个样本x。对非增量学习算法来说，我们需要把x和数据库 A混在一起，组成新的数据库B，再重新训练新的分类器。但对增量学习算法，我们只需要用新样本x来更新已有分类器h的参数即可），所以它属于在线学习算 法。与在线学习相对应，一次处理整个数据集的叫“批处理”。
随机梯度下降算法的伪代码如下：
初始化回归系数为1
重复下面步骤直到收敛{
对数据集中每个样本
计算该样本的梯度
使用alpha xgradient来更新回归系数
}
返回回归系数值
1 def stocGradAscent0(dataMatrix, classLabels): 2 m,n = shape(dataMatrix) 3 alpha = 0.01 4 weights = ones(n) #initialize to all ones 5 for i in range(m): 6 h = sigmoid(sum(dataMatrix[i]*weights)) 7 error = classLabels[i] - h 8 weights = weights + alpha * error * dataMatrix[i] 9 return weights 10 #dataMat,labelMat=loadDataSet() 11 #weights=stocGradAscent0(array(dataMat), labelMat) 12 #plotBestFit(weights)
stocGradAscent0
3.5 随机梯度上升法的改进
1 def stocGradAscent1(dataMatrix, classLabels, numIter=150): 2 m,n = shape(dataMatrix) 3 weights = ones(n) #initialize to all ones 4 for j in range(numIter): 5 dataIndex = range(m) 6 for i in range(m): 7 alpha = 4/(1.0+j+i)+0.0001 #apha decreases with iteration, does not 8 randIndex = int(random.uniform(0,len(dataIndex)))#go to 0 because of the constant 9 h = sigmoid(sum(dataMatrix[randIndex]*weights)) 10 error = classLabels[randIndex] - h 11 weights = weights + alpha * error * dataMatrix[randIndex] 12 del(dataIndex[randIndex]) 13 return weights 14 #dataMat,labelMat=loadDataSet() 15 #weights=stocGradAscent1(array(dataMat), labelMat) 16 #plotBestFit(weights)
stocGradAscent1
对随机梯度下降算法，我们做两处改进来避免上述的波动问题：
1） 在每次迭代时，调整更新步长alpha的值。随着迭代的进行，alpha越来越小，这会缓解系数的高频波动（也就是每次迭代系数改变得太大，跳的跨度太 大）。当然了，为了避免alpha随着迭代不断减小到接近于0（这时候，系数几乎没有调整，那么迭代也没有意义了），我们约束alpha一定大于一个稍微 大点的常数项，具体见代码。
2）每次迭代，改变样本的优化顺序。也就是随机选择样本来更新回归系数。这样做可以减少周期性的波动，因为样本顺序的改变，使得每次迭代不再形成周期性。
改进的随机梯度下降算法的伪代码如下：
初始化回归系数为1
重复下面步骤直到收敛{
对随机遍历的数据集中的每个样本
随着迭代的逐渐进行，减小alpha的值
计算该样本的梯度
使用alpha x gradient来更新回归系数
}
返回回归系数值
4. 示例：从疝气病症预测病马的死亡率
4.1 用logistic回归进行分类的时候，处理缺失值的方法：选择实数0来代替。
处理无标签的数据：丢弃
4.2 用logistic回归进行分类：1）把测试集上每个特征向量乘以最优方法得来的回归系数
2）再将该乘积结果求和
3）最后输入到sigmoid函数即可，如果函数值大于0.5就预测类别为1，反之为0.
1 def classifyVector(inX, weights): 2 prob = sigmoid(sum(inX*weights)) 3 if prob > 0.5: return 1.0 4 else: return 0.0 5 6 def colicTest(): 7 frTrain = open(‘horseColicTraining.txt‘); frTest = open(‘horseColicTest.txt‘) 8 trainingSet = []; trainingLabels = [] 9 for line in frTrain.readlines(): 10 currLine = line.strip().split(‘\t‘) 11 lineArr =[] 12 for i in range(21): 13 lineArr.append(float(currLine[i])) 14 trainingSet.append(lineArr) 15 trainingLabels.append(float(currLine[21])) 16 trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 1000) 17 errorCount = 0; numTestVec = 0.0 18 for line in frTest.readlines(): 19 numTestVec += 1.0 20 currLine = line.strip().split(‘\t‘) 21 lineArr =[] 22 for i in range(21): 23 lineArr.append(float(currLine[i])) 24 if int(classifyVector(array(lineArr), trainWeights))!= int(currLine[21]): 25 errorCount += 1 26 errorRate = (float(errorCount)/numTestVec) 27 print "the error rate of this test is: %f" % errorRate 28 return errorRate 29 30 def multiTest(): 31 numTests = 10; errorSum=0.0 32 for k in range(numTests): 33 errorSum += colicTest() 34 print "after %d iterations the average error rate is: %f" % (numTests, errorSum/float(numTests))
example