Logistic回归

和很多其他机器学习算法一样，逻辑回归也是从统计学中借鉴来的，尽管名字里有回归俩字儿，但它不是一个需要预测连续结果的回归算法。
与之相反，Logistic 回归是二分类任务的首选方法。它输出一个 0 到 1 之间的离散二值结果。简单来说，它的结果不是 1 就是 0。
癌症检测算法可看做是 Logistic 回归问题的一个简单例子，这种算法输入病理图片并且应该辨别患者是患有癌症（1）或没有癌症（0）。

• Logistic是如何工作的:
  Logistic 回归通过使用其固有的 logistic 函数估计概率，来衡量因变量（我们想要预测的标签）与一个或多个自变量（特征）之间的关系。
  然后这些概率必须二值化才能进行预测。这就是 logistic 函数的任务，也称为 sigmoid 函数。Sigmoid 函数是一个 S 形曲线，它可以将任意实数值映射到介于 0 和 1 之间的值，但并不会取到 0/1。然后使用阈值分类器将 0 和 1 之间的值转换为 0 或 1。
  我们希望随机数据点被正确分类的概率最大化，这就是最大似然估计。最大似然估计是统计模型中估计参数的通用方法。
  还可以使用不同的方法（如优化算法）来最大化概率。牛顿法也是其中一种，可用于查找许多不同函数的最大值（或最小值），包括似然函数。也可以用梯度下降法代替牛顿法。
• Logistic回归的公式
  一般来说，回归模型一般不用在分类问题上，因为回归是连续型模型，而且受噪声的因素很大，但是，若需要选择，可以选择使用logisti 回归。
  对数回归本质上是线性回归，只是在特征到结果的映射里加入了一层函数映射，选择g(z)=1/(1+exp（-z）)作为sigmoid函数进行映射，可以将连续值映射到0-1之间。
  其中g(z)函数的图像如下：可以看到，函数的取值始终在0-1之间。
  对于分类问题，我们可以建立假设：
  if（z >= 0.5） g(z)=1; if（z < 0.5） g(z)=0;

算法思想：
对于输出值为y={0,1}的两类分类问题，我们做出一个假设

函数g(z)即为上文提到的sigmoid函数，其导数形式为：

根据这个函数，我们可以得到对于一个样本的概率分布为：

综合起来就是：

将问题转化为求Logistic回归的最佳系数，因为logistic回归可以被看做一种概率模型，且y发生的概率与回归参数Θ有关，因此我们可以对Θ进行最大似然估计，使得y发生的概率最大，此时的Θ 便是最优回归系数：

对数据集求似然函数，并取对数计算：

要使得最大化，则运用梯度上升法，求出最高点：

计算结果为：

此公式便是梯度上升算法的更新规则，α是学习率，决定了梯度上升的快慢。可以看到与线性回归类似，只是增加了特征到结果的映射函数。
sigmoid函数中Z的计算公式为：
是一个矩阵，θ是参数列向量(要求解的)，x是样本列向量(给定的数据集)。θ^T表示θ的转置。g(z)函数实现了任意实数到[0,1]的映射，这样我们的数据集([x0,x1,…,xn])，不管是大于1或者小于0，都可以映射到[0,1]区间进行分类。hθ(x)给出了输出为1的概率。比如当hθ(x)=0.7，那么说明有70%的概率输出为1。输出为0的概率是输出为1的补集，也就是30%。
如果我们有合适的参数列向量θ([θ0,θ1,…θn]^T)，以及样本列向量x([x0,x1,…,xn])，那么我们对样本x分类就可以通过上述公式计算出一个概率，如果这个概率大于0.5，我们就可以说样本是正样本，否则样本是负样本。
一、加载数据

def loadDataSet():
    """
    加载数据
    :return:
      dataMat - 数据列表
    labelMat - 标签列表
    """
    dataMat = []  # 创建数据列表
    labelMat = [] #创建标签列表
    fr = open('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])) #添加标签

    fr.close() #关闭文件
    return dataMat, labelMat #返回

为了方便计算，我们将 X0 的值设为 1.0 ，也就是在每一行的开头添加一个 1.0 作为 X0。因为假设Sigmoid函数的输入记为z，那么z=w0x0 + w1x1 + w2x2，即可将数据分割开。其中，x0为全是1的向量，x1为数据集的第一列数据，x2为数据集的第二列数据。另z=0，则0=w0 + w1x1 + w2x2。横坐标为x1，纵坐标为x2。这个方程未知的参数为w0，w1，w2，也就是我们需要求的回归系数(最优参数)。
二、 训练算法
通过上面的公式推导我们可以得到下梯度上升迭代公式：
将上述公式矢量化：
根据矢量化的公式，编写训练代码如下：

#定义sigmoid函数
def sigmoid(inX):
    return 1.0/(1+exp(-inX)) #输入数据特征与数据的类别标签
    # #返回最佳回归系数(weights)
def gradAscent(dataMatIn, classLabels):
    """
    梯度上升算法
    :param dataMatIn: 数据集
    :param classLabels: 数据标签
    :return:
    weights.getA() - 求得的权重数组(最优参数)
    """

    dataMatrix = mat(dataMatIn) #转换为numpy型
    labelMat = mat(classLabels).transpose()#  转化为矩阵[[0,1,0,1,0,1.....]]，并将行向量转化为列向量[[0],[1],[0].....]  
    m,n = shape(dataMatrix)# m->数据量，样本数 n->特征数
    alpha = 0.001 #步长
    maxCycles = 500 #迭代次数 
    weights = ones((n,1)) #初始化权值向量，每个维度均为1.0
    for k in range(maxCycles): #求当前sigmoid函数的值
        h = sigmoid(dataMatrix * weights) 
        error = (labelMat - h)
        weights = weights + alpha * dataMatrix.transpose() * error
    return array(weights)

三、随机梯度上升
每次更新回归系数都要遍历整个数据集，如果数据量庞大的话，该方法的计算复杂度就很高。一种改进的方法是一次仅使用一个样本点来更新回归系数，我们叫做随机梯度上升法，并且它是一种在线学习算法。前面的梯度上升法我们叫做(Batch learning)批处理

def stoGradAscent0(dataMatrix, classLabels):
	 m, n = shape(dataMatrix) 
	 alpha = 0.01 
	 weights = ones(n)
	 for i in range(m):
	 	h = sigmoid(sum(dataMatrix[i] * weights))
	 	error = classLabels[i] - h 
	 	weights = weights + alpha * error * dataMatrix[i]
  return weights

两者在代码上的区别：

1.后者的变量h和误差error都是向量

2.前者没有矩阵的转换过程，所有变量的数据类型都是Numpy数组

可视化：

从结果上来看，该方法不如直接使用梯度上升法那样完美。

直接这样比较是不公平的，我们知道在梯度上升法中迭代了500次。

一个判断优化算法优劣的可靠方法是看它是否收敛，也就是说参数是否达到了稳定值，是否还会不断的变化？通过分析，我们使用改进的随机梯度上升算法

def stocGradAscent1(dataMatrix, classLabels, numIter=150): 
	m,n = np.shape(dataMatrix) #返回dataMatrix的大小。m为行数,n为列数。 
	weights = np.ones(n) #参数初始化
	 for j in range(numIter): 
	 	dataIndex = list(range(m)) 
	 	for i in range(m): 
	 		alpha = 4/(1.0+j+i)+0.01 #降低alpha的大小，每次减小1/(j+i)。
	 		 randIndex = int(random.uniform(0,len(dataIndex))) #随机选取样本
	 		 h = sigmoid(sum(dataMatrix[randIndex]*weights)) #选择随机选取的一个样本，计算h 
	 		 error = classLabels[randIndex] - h #计算误差
	 		 weights = weights + alpha * error * dataMatrix[randIndex] #更新回归系数 
	 		 del(dataIndex[randIndex]) #删除已经使用的样本
  return weights #返回

该算法第一个改进之处在于，alpha在每次迭代的时候都会调整，并且，虽然alpha会随着迭代次数不断减小，但永远不会减小到0，因为这里还存在一个常数项。必须这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。如果需要处理的问题是动态变化的，那么可以适当加大上述常数项，来确保新的值获得更大的回归系数。另一点值得注意的是，在降低alpha的函数中，alpha每次减少1/(j+i)，其中j是迭代次数，i是样本点的下标。第二个改进的地方在于跟新回归系数(最优参数)时，只使用一个样本点，并且选择的样本点是随机的，每次迭代不使用已经用过的样本点。这样的方法，就有效地减少了计算量，并保证了回归效果。
五、示例 从疝气病症预测病马的死亡率
由于数据中会有缺失值，因此我们要处理数据中的缺失值
这里只是列了几个方法，其实数据清洗以及特征提取是一项很庞大的工程。

1.使用可用特征的均值来填补缺失值
2.使用特殊值来填补缺失值
3.忽略有缺失值的样本
4.使用相似样本的均值填补缺失值
5.使用另外的机器学习算法来预测缺失值

  def classifyVector(inX, weights):
    prob = sigmoid(sum(inX * weights)) #大于0.5 返回 1；否则返回0
    if prob > 0.5:
        return 1.0
    else:
      return 0.0
def colicTest():
    frTrain = open('HorseColicTraining.txt')
    frTest = open('HorseColicTest.txt')
    trainingSet = []
    trainingLabels = [] # 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])) #使用改进后的随机梯度下降算法得到回归系数
    trainingWeights =stocGradAscent1(array(trainingSet), trainingLabels, 500) #统计测试集预测错误样本数量和样本总数
    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])) # 利用分类预测函数对该样本进行预测，并与样本标签进行比较
         if int(classifyVector(array(lineArr), trainingWeights)) != int(currLine[21]): #如果预测错误，错误数加1
            errorCount += 1 #计算错误率
    errorRate = (float(errorCount) / numTestVec)
    print('the error rate of this test is : %f' % errorRate)
    return errorRate
def multiTest():
    numTests = 10
    errorSum = 0.0
    for k in range(numTests):
        errorSum += colicTest()
    print('after %d iterations the average error rete is : %f ' %(numTests,errorSum / float(numTests))) 

此外，我们通过调整步长与迭代次数该结果肯定也会有改变。