逻辑回归（Logistic Regression)详解

Posted 2023-04-07 生信小兔

逻辑回归也称作logistic回归分析，是一种广义的线性回归分析模型，属于机器学习中的监督学习。其推导过程与计算方式类似于回归的过程，但实际上主要是用来解决二分类问题（也可以解决多分类问题）。通过给定的n组数据（训练集）来训练模型，并在训练结束后对给定的一组或多组数据（测试集）进行分类。其中每一组数据都是由p 个指标构成。

（1）逻辑回归所处理的数据

逻辑回归是用来进行分类的。例如，我们给出一个人的 [身高，体重] 这两个指标，然后判断这个人是属于”胖“还是”瘦“这一类。对于这个问题，我们可以先测量n个人的身高、体重以及对应的指标”胖“,"瘦”，把胖和瘦分别用0和1来表示，把这n组数据输入模型进行训练。训练之后再把待分类的一个人的身高、体重输入模型中，看这个人是属于“胖”还是“瘦”。

如果数据是有两个指标，可以用平面的点来表示数据，其中一个指标为x轴，另一个为y轴；如果数据有三个指标，可以用空间中的点表示数据；如果是p维的话(p>3)，就是p维空间中的点。

从本质上来说，逻辑回归训练后的模型是平面的一条直线（p=2),或是平面（p=3)，超平面（p>3)。并且这条线或平面把空间中的散点分成两半，属于同一类的数据大多数分布在曲线或平面的同一侧。

 如上图所示，其中点的个数是样本个数，两种颜色代表两种指标。这个直线可以看成经这些样本训练后得出的划分样本的直线。那么对于之后的样本的p1与p2的值，就可以根据这条直线来判断它属于哪一类了。

（2）算法原理

首先，我们处理二分类问题。由于分成两类，我们便让其中一类标签为0，另一类为1。我们需要一个函数，对于输入的每一组数据，都能映射成0~1之间的数。并且如果函数值大于0.5，就判定属于1，否则属于0。而且函数中需要待定参数，通过利用样本训练，使得这个参数能够对训练集中的数据有很准确的预测。

这个函数就是sigmoid函数，形式为。所以在这里我们可以设函数为

这里是测试集第i个数据，是p维列向量；w是p维列向量,为待求参数；b是一个数，也是待求参数。

我们发现，对于,其结果是。所以我们可以把

w写成,写成。就可以写成,则：

这样就可以把另一个参数b合并到w中。后面推导也方便很多。当然，我们也可以用用第一种形式来做，本质是相同的。之后就是根据训练样本来求参数w了。

（3）求解参数。

这一部分便是逻辑回归的核心问题了。兔兔在下面将给出两种方法。

（1）极大似然估计。

极大似然估计是数理统计中参数估计的一种重要方法。其思想就是一个事件发生了，那么发生这个事件的概率就是最大的。对于样本i,其类别为。对于样本i，可以把看成是一种概率。yi对应是1时，概率是h(xi)，即xi属于1的可能性;yi对应是0时，概率是1-h(xi)，即xi属于0的可能性 。那么它构造极大似然函数

.

其中i从0到k是属于类别1的个数k,i从k+1到n是属于类别0的个数n-k。由于y是标签0或1，所以上面的式子也可以写成：

这样无论y是0还是1，其中始终有一项会变成0次方，也就是1，和第一个式子是等价的。

为了方便，我们对式子取对数。因为是求式子的最大值，可以转换成式子乘以负1,之后求最小值。同时对于n个数据，累加后值会很大，之后如果用梯度下降容易导致梯度爆炸。所以可以除以样本总数n。

求最小值方法很多，机器学习中常用梯度下降系列方法。也可以采用牛顿法，或是求导数为零时w的数值等。

（2）损失函数

逻辑回归中常用交叉熵损失函数，交叉熵损失函数和上面极大似然法得到的损失函数是相同的。这里不再赘述。另一种也可以采用平方损失函数(均方误差），即

这个是比较直观的。就是让这个预测函数h(xi)与实际的分类1或0越接近越好。也就是损失函数越小越好。求最小值还是用到上面说到的方法。

目前我们这里已经得到了这两种函数。我们以梯度下降为例，即求损失函数的导数。

对于损失函数（1），导数求解过程如下（需要用到矩阵求导）。

对于损失函数（2），推导过程如下：

 

（3）算法实现。

兔兔在这里以Dry_Bean_Dataset文件为例。同学们可以在 www.kaggle.com 中下载该数据集，应该比较方便的。

import pandas as pd
data=pd.read_csv('Dry_Bean_Dataset.csv')
df=pd.DataFrame(data)
print(df.columns，df.shape)

我们先看一下数据集的情况，里面有很多个指标，包括菜豆的区域，周长，长轴长、短轴长等。class代表对应菜豆的种类。我们这里为了直观方便（便于在二维平面表示），就选MajorAxisLength,MinorAxisLength这两个指标。由于是二分类，我们就选SEKER和BARBUNYA这两个种类，即前3349组数据。绘制散点图观察。

import pandas as pd
import matplotlib.pyplot as plt
data=pd.read_csv('Dry_Bean_Dataset.csv')
df=pd.DataFrame(data)
color=[]
for i in df['Class'][0:3349]:
    if i=='SEKER':
        color.append('red')
    else:
        color.append('blue')
plt.scatter(df['MajorAxisLength'][0:3349],df['MinorAxisLength'][0:3349],color=color)
plt.xlabel('MajorAxisLength')
plt.ylabel('MinorAxisLength')
plt.show()

散点图如下：

 可见，该数据是可以用一条直线分成两个部分的。

之后，我们需要处理一下数据。把每组数据两个指标及其对应的分类绑在一起（如果用随机梯度下降需要这样处理，批量梯度下降可以不用这样处理），并把每组数据转成列向量。

import numpy as np
import pandas as pd
data=pd.read_csv('Dry_Bean_Dataset.csv')
df=pd.DataFrame(data)
label=[]
for i in df['Class'][0:3349]:
    if i=='SEKER':
        label.append(0)
    else:
        label.append(1)
x1=df['MajorAxisLength'][0:3349]
x2=df['MinorAxisLength'][0:3349]
train_data=list(zip(x1,x2,label))
class Logistic_Regression:
    def __init__(self,traindata,alpha=0.001,circle=1000,batchlength=40):
        self.traindata=traindata #训练数据集
        self.alpha=alpha #学习率
        self.circle=circle #学习次数
        self.batchlength=batchlength #把3349个数据分成多个部分，每个部分有batchlength个数据
        self.w=np.random.normal(size=(3,1)) #随机初始化参数w
    def data_process(self):
        '''做随机梯度下降，打乱数据顺序，并把所有数据分成若干个batch'''
        np.random.shuffle(self.traindata)
        data=[self.traindata[i:i+self.batchlength]
              for i in range(0,len(self.traindata),self.batchlength)]
        return data
    def train1(self):
        '''根据损失函数（1）来进行梯度下降，这里采用随机梯度下降'''
        for i in range(self.circle):
            batches=self.data_process()
            print('the  epoch'.format(i)) #程序运行时显示执行次数
            for batch in batches:
                d_w=np.zeros(shape=(3,1)) #用来累计w导数值
                for j in batch: #取batch中每一组数据
                    x0=np.r_[j[0:2],1] #把数据中指标取出，后面补1
                    x=np.mat(x0).T #转化成列向量
                    y=j[2] #标签
                    dw=(self.sigmoid(self.w.T*x)-y)[0,0]*x
                    d_w+=dw
                self.w-=self.alpha*d_w/self.batchlength

    def train2(self):
        '''用均方损失函数来进行梯度下降求解'''
        for i in range(self.circle):
            batches=self.data_process()
            print('the  epoch'.format(i)) #程序运行时显示执行次数
            for batch in batches:
                d_w=np.zeros(shape=(3,1)) #用来累计w导数值
                for j in batch: #取batch中每一组数据
                    x0=np.r_[j[0:2],1] #把数据中指标取出，后面补1
                    x=np.mat(x0).T #转化成列向量
                    y=j[2] #标签
                    dw=((self.sigmoid(self.w.T*x)-y)*self.sigmoid(self.w.T*x)*(1-self.sigmoid(self.w.T*x)))[0,0]*x
                    d_w+=dw
                self.w-=self.alpha*d_w/self.batchlength
   
    def sigmoid(self,x):
        return 1/(1+np.exp(-x))

    def predict(self,x):
        '''测试新数据属于哪一类，x是2维列向量'''
        s=self.sigmoid(self.w.T*x)
        if s>=0.5:
            return 1
        elif s<0.5:
            return 0
if __name__=='__main__':
    regr=Logistic_Regression(traindata=train_data)
    regr.train1() #采用1的方式进行训练

处理过程需要注意数据的类型。例如当行向量与列向量相乘为一个数时，我们计算是一个数，但是numpy 返回是一个1x1的矩阵，这样再与后面向量相乘就会出错，所以需要注意把矩阵转换成数。

那么，对于这个模型如何比较直观地观察效果呢？我们知道，最终x代入h(x)中，值越接近于1说明越可能属于1类，越接近0越属于0类；那么当h(x)=0.5时，这个时候x就在分界线上面了。所以由，最终得出，即。这个便是分界直线方程。

w=regr.w
w1=w[0,0]
w2=w[1,0]
w3=w[2,0]
x=np.arange(190,500)
y=-w1*x/w2-w3/w2
plt.plot(x,y)
plt.show()

为了动态观察直线变换，也是可以把这个部分放在train()函数当中的。

对于上面的运行结果，我们可以发现效果不是很好。一方面运行速度很慢，而且alpha设置稍稍大一点很容易造成sigmoid函数下方溢出现象，另一方面，我们发现直线并未把两组数据完全分割开，甚至是穿过了两个区域。

关于这个现象，兔兔认为：数据集数量较大，每处理一遍数据都要很大的计算量，导致运行速度慢，解决方法是：可以每次circle只从中随机选取部分数据进行训练；对于函数溢出的现象，在于运行过程中变成很小的负数，导致指数运算过大，所以可以把学习率变小，防止一次迭代后w就变成使得很小的数值；直线没有完全分割两个部分，我们可以从数据特征方面来考虑，由于我们人为地随便选两个指标，然后根据这两个指标训练并判断分类标准，这样做是很容易出现问题的，我们不知道两个指标的数据特征与关系，也不清楚其它的指标是否对分类起决定性作用，而且两类数据有很多公共的部分重叠，也同样导致分类出现问题。当我们训练结束后，发现直线虽然穿过两个区域，是斜率为正的直线，但实际上它也把两个区域的核心部分（两类当中数据最集中部分）分成了两个部分，说明这一部分数据起了主要作用，而且逻辑回归训练最终目的也是让损失函数达到最小，说明最终的曲线也是符合要求的。

 为了直观体现逻辑回归的训练特点，兔兔选了iris数据集一部分进行训练，效果如下图所示。

 像上面的情况基本上循环几十次就稳定在最优解附近。如果学习率等参数没有弄好，可能会出现下面情况：

 （4）非线性逻辑回归

与线性的逻辑回归相比，非线性的逻辑回归应用应该更加普遍。例如，当两组数据无法用一条直线或一个平面分割，而是需要曲线或曲面才能分割好。这个时候就可以用非线性的逻辑回归了。比如用一个圆、椭圆、曲线等把两组数据分开。

非线性回归的训练和推导过程和前面的是一致的。只是把x1,x2两个指标做一下处理。这个过程和兔兔之前讲到的非线性回归是一致的（详见：《线性回归（线性拟合）与非线性回归（非线性拟合）原理、推导与算法实现（一）》）。

前面我们发现，我们最后训练得到的曲线方程是。那么，如果我们令,把每组输入的向量x处理成：,也就是类似多项式回归。这样训练之后就可以用形式为这样的曲线方程来分割区域了。对于三维、p维也是如此。我们甚至也可以根据需要调整多项式次数，或者函数形式等，从而达到理想的效果。但是这时要注意过拟合、欠拟合的情况发生，并且需要正则化等处理方法。

（5）逻辑回归的多分类问题。

前面所讲的都是逻辑回归的二分类问题，那么逻辑回归是否可以处理多分类呢？答案是肯定的。这时我们就不再使用sigmoid 函数，而是另一个叫做softmax函数。函数形式如下：

那么h(x)函数就是

这里同样，我们把k个类用数字1，2.....k来表示，在sigmoid 函数当中，函数值表示概率。在这里也是如此，x经过h函数处理后，得到向量里面对应位置（分类）的数值就是取对应位置（分类）的概率。例如，对于三分类问题，如果向量中p(y=1)=0.7,p(y=2)=0.2,p(y=3)=0.1,那么x属于1类的概率最大，故判别为1类。

它的推导过程于前面类似，同样构造损失函数，求损失函数对w的导数，做梯度下降处理计算。

 

这里用的是均方损失，推导过程如上所示，对各个w求偏导，之后相应的w做梯度下降即可。

（6）总结。

逻辑回归可以分为线性与非线性，也可以根据类的个数分为二分类与多分类问题，使用时需要灵活应用，能够构造损失函数并求梯度，同时能够用算法实现并进行训练预测。

事实上，细心的同学会发现，在逻辑回归中，我们发现是多个输入（即p个指标），最终输出一个结果（0或1），处理过程是输入乘上权重w加偏置b（本文权重w与偏置b都合并到w中了），再对结果用sigmoid 函数处理，这个过程其实很接近于神经网络了，而逻辑回归的模型更接近于感知机。对于神经网络，它不只有输入和输出两层，而且增加了更多的隐藏层，每一层的处理结果都作为下一层的输入，那么它的损失函数与梯度的求解也将更加复杂，模型也复杂许多。

逻辑回归（logistic regression）

1、简介

回归是一种极易理解的模型，就相当于y=f(x)，表明自变量x与因变量y的关系。最简单的回归是线性回归，然而线性回归的鲁棒性很差，使回归模型在训练集上表现都很差。这主要是由于线性回归在整个实数域内敏感度一致，而分类范围，需要在[0,1]。逻辑回归就是一种减小预测范围，将预测值限定为[0,1]间的一种回归模型。

2、模型

逻辑回归其实仅为在线性回归的基础上，套用了一个逻辑函数，首先介绍逻辑斯谛分布，该分布的定义是：设X是连续随机变量,X服从逻辑斯谛分布是指X服从如下分布函数和密度函数:

其中，为位置参数，> 0 为形状参数。

可以通过其图像观察：

右边的逻辑斯蒂分布函数以点中心对称，即满足：

形状参数越小，曲线在中心的增长速度越快。

二项逻辑斯蒂回归模型

这是一种由条件概率表示的模型，其条件概率模型如下：

其中，exp为以e为底的指数函数，x∈Rn是输入，y∈{0,1}输出，w，b是模型参数——w是权值向量，b称作偏置，w·x是向量内积。

有了后验概率，逻辑斯蒂回归模型选择二分类中较大的那一个完成分类。

另外，逻辑斯特回归模型还有一个方便的形式，如果将权值向量w和输入向量x拓充为w=(w(1),w(2),…w(n),b)T,x=(x(1),…x(n),1)T，此时逻辑斯谛模型可以表示为：

为什么要重新提一个形式出来呢？这是因为，这个形式跟几率的表达式很像。

定义事件的几率:发生概率与不发生概率的比值——。

定义对数几率：

 

将逻辑斯蒂模型的便捷形式做一个变换恰好可以得到：

 

这也就是说，在逻辑斯蒂回归模型中，输出Y=1的对数几率是输入x的线性函数。或者说输出Y=1的对数几率是由输入x的线性函数表示的模型，即逻辑斯蒂回归模型。反过来讲，如果知道权值向量，给定输入x，就能求出Y=1的概率：

 

线性函数w·x的值越接近正无穷，概率值越接近1；反之，越接近负无穷，概率值越接近0——这就是逻辑斯谛回归模型。

模型参数估计

在模型学习的时候，对于给定训练集T = {(x1,y1)…(xN,yN)},x∈Rn,y∈{0,1}

设

定义似然函数

则有对数似然函数

这个好说，把后面括号里的负π提到前面去就行了。

对L(w)求极大值就可以得出权值向量w的估计值。

解决以L(w)为目标函数的最优化问题的一般方法是梯度下降法及拟牛顿法。

梯度下降算法

求函数最小值的时候用的是梯度下降算法，而此处求的是对数似然函数的最大值，所以应该称为梯度上升算法。函数的梯度由其偏导数构成：

梯度是函数增长最快的方向，记移动补偿为α，则梯度算法的迭代公式为：

假定权值向量w有了，怎么计算模型输出呢？

特征向量乘以权值向量得出一个实数z

希望通过该实数输出一个0或1的类别，这时候就需要利用Sigmoid函数了：

其图像如下：

将该实数代入Sigmoid函数后，得到一个0~1之间的数，大于0.5归入1，小于0.5归入0即可。

利用Sigmoid函数，梯度上升算法的代码如下：

from numpy import *
 
def sigmoid(inX):
    return 1.0/(1+exp(-inX))
 
def gradAscent(dataMatIn, classLabels):
    """
    逻辑斯谛回归梯度上升优化算法
    :param dataMatIn:输入X矩阵（100*3的矩阵，每一行代表一个实例，每列分别是X0 X1 X2）
    :param classLabels: 输出Y矩阵（类别标签组成的向量）
    :return:权值向量
    """
    dataMatrix = mat(dataMatIn)             #转换为 NumPy 矩阵数据类型
    labelMat = mat(classLabels).transpose() #转换为 NumPy 矩阵数据类型
    m,n = shape(dataMatrix)                 #矩阵大小
    alpha = 0.001                           #步长
    maxCycles = 500
    weights = ones((n,1))
    for k in range(maxCycles):              #最大迭代次数
        h = sigmoid(dataMatrix*weights)     #矩阵内积
        error = (labelMat - h)              #向量减法
        weights += alpha * dataMatrix.transpose() * error  #矩阵内积
    return weights

 

随机梯度上升算法

梯度下降算法在每次更新权值向量的时候都需要遍历整个数据集，该方法对小数据集尚可。但如果有数十亿样本和成千上万的特征时，它的计算复杂度就太高了。一种改进的方法是一次仅用一个样本点的回归误差来更新权值向量，这个方法叫随机梯度下降算法。由于可以在遇到新样本的时候再对分类器进行增量式更新，所以随机梯度上升算法是一个在线学习算法；与此对应，一次处理完所有数据的算法（如梯度上升算法）被称作“批处理”。

随机梯度上升算法的代码：

def stocGradAscent0(dataMatrix, classLabels, history_weight):
    """
    随机梯度上升算法
    :param dataMatIn:输入X矩阵（100*3的矩阵，每一行代表一个实例，每列分别是X0 X1 X2）
    :param classLabels: 输出Y矩阵（类别标签组成的向量）
    :return:权值向量
    """
    dataMatrix = array(dataMatrix)
    m,n = shape(dataMatrix)
    alpha = 0.01
    weights = ones(n)                               #初始化为单位矩阵
    for i in range(m):
        h = sigmoid(sum(dataMatrix[i]*weights))     #挑选（伪随机）第i个实例来更新权值向量
        error = classLabels[i] - h
        weights = weights + dataMatrix[i] * alpha * error
        history_weight.append(copy(weights))
    return weights

 

改进的随机梯度上升算法

既然随机梯度上升算法最终给出的参数不好，那是否仅仅是因为参数没有足够收敛，而算法本质是优秀的呢？对此，可以逐步减小步长，避免参数周期性的抖动。代码：

def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    """
    改进的随机梯度上升算法
    :param dataMatIn:输入X矩阵（100*3的矩阵，每一行代表一个实例，每列分别是X0 X1 X2）
    :param classLabels: 输出Y矩阵（类别标签组成的向量）
    :param numIter: 迭代次数
    :return:
    """
    dataMatrix = array(dataMatrix)
    m,n = shape(dataMatrix)
    weights = ones(n)                                           #初始化为单位矩阵
    for j in range(numIter):
        dataIndex = range(m)
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.0001                          #步长递减，但是由于常数存在，所以不会变成0
            randIndex = int(random.uniform(0,len(dataIndex)))   #总算是随机了
            h = sigmoid(sum(dataMatrix[randIndex]*weights))
            error = classLabels[randIndex] - h
            weights = weights + alpha * error * dataMatrix[randIndex]
            del(dataIndex[randIndex])                           #删除这个样本，以后就不会选到了
    return weights

以上三种算法的效果：

此外，还有多元逻辑回归，假设离散型随机变量Y的取值集合是{1,2,…K}，那么多项逻辑斯蒂回归模型是： 
 
对于参数估计方法仍然可以采用最大似然估计。

3、总结

个人感觉，只要特征找的准，数据量足够大，逻辑斯蒂回归将会非常好用。另外，还要注意避免过拟合。

特征选择的话，由于逻辑斯蒂回归的优点，开始的时候不用考虑各个特征之间是否有相关性，直接把能用的特征全部线性加权起来就好。经过初步训练，观察各个特征的权值，如果权值接近为0，那么就可以将这个特征看做是不相关的可以去除的特征。总结起来就是：先做加法再做减法。

解决过拟合的方法不过两种，一种是减少特征的个数；另一种是模型选择的正则化方法。正则化的话，可以参考岭回归方法。逻辑回归的优缺点总结如下：优点：计算代价不高，易于理解和实现，且若采用随机梯度上升法可以在线学习； 缺点：可能容易欠拟合，分类精度不高，这个可能是因为我们无法找到足够的特征。

4、实现

依据以上分析，得到如下代码：

‘‘‘
Created on Oct 27, 2010
Logistic Regression Working Module
@author: Peter
‘‘‘
from numpy import *

def loadDataSet():
    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]))
    return dataMat,labelMat

def sigmoid(inX):
    return 1.0/(1+exp(-inX))

def gradAscent(dataMatIn, classLabels):
    dataMatrix = mat(dataMatIn)             #convert to NumPy matrix
    labelMat = mat(classLabels).transpose() #convert to NumPy matrix
    m,n = shape(dataMatrix)
    alpha = 0.001
    maxCycles = 500
    weights = ones((n,1))
    for k in range(maxCycles):              #heavy on matrix operations
        h = sigmoid(dataMatrix*weights)     #matrix mult
        error = (labelMat - h)              #vector subtraction
        weights = weights + alpha * dataMatrix.transpose()* error #matrix mult
    return weights

def plotBestFit(weights):
    import matplotlib.pyplot as plt
    dataMat,labelMat=loadDataSet()
    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(-3.0, 3.0, 0.1)
    y = (-weights[0]-weights[1]*x)/weights[2]
    ax.plot(x, y)
    plt.xlabel(‘X1‘); plt.ylabel(‘X2‘);
    plt.show()

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

def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    m,n = shape(dataMatrix)
    weights = ones(n)   #initialize to all ones
    for j in range(numIter):
        dataIndex = range(m)
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.0001    #apha decreases with iteration, does not 
            randIndex = int(random.uniform(0,len(dataIndex)))#go to 0 because of the constant
            h = sigmoid(sum(dataMatrix[randIndex]*weights))
            error = classLabels[randIndex] - h
            weights = weights + alpha * error * dataMatrix[randIndex]
            del(dataIndex[randIndex])
    return weights

def classifyVector(inX, weights):
    prob = sigmoid(sum(inX*weights))
    if prob > 0.5: return 1.0
    else: return 0.0

def colicTest():
    frTrain = open(‘horseColicTraining.txt‘); frTest = open(‘horseColicTest.txt‘)
    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 = stocGradAscent1(array(trainingSet), trainingLabels, 1000)
    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

def multiTest():
    numTests = 10; errorSum=0.0
    for k in range(numTests):
        errorSum += colicTest()
    print "after %d iterations the average error rate is: %f" % (numTests, errorSum/float(numTests))
        

相关博客：

http://www.hankcs.com/ml/the-logistic-regression-and-the-maximum-entropy-model.html

https://www.douban.com/note/323644915/

http://www.cnblogs.com/zhizhan/p/5007540.html