逻辑斯蒂拓展以及Python实现

为什么使用sigmoid函数

在上两节中已经说到，逻辑斯蒂回归只是在线性回归的基础上套用了一个逻辑函数，并且也解释了，线性的原因。忘记的可以回顾一下。

下面我们再具体解释一下，为什么逻辑斯蒂选择使用sigmoid函数。

为什么是sigmoid函数

假设数据集有n个独立的 特征，x1到xn为样本的n个特征。常规的回归算法的目标是拟合出一个多项式函数（线性函数），使得预测值与真实值的误差最小：

但是线性回归的鲁棒性比较差，可能会因为少数的噪音，导致模型不好，主要是因为线性模型对整个实数域的敏感性是一致的。我们需要解决的是：对于大部分集中的数据是敏感的，而对于边缘的少数数据是不敏感的。另一方面，我们希望我们的模型能够具有很好的逻辑判断性质，最好是能够直接表达具有特征X的样本被分到某类的概率是多大，把概率范围控制在[0,1]。比如，假设我们的模型为f(x)，当f(x)>0.5的时候，把样本X归为正类；当f(x)<0.5的时候，把样本X归为负类。那么考虑有没有这样的一种函数，可以满足我们以上的两种需求呢？

这个时候我们的救星出现了：sigmoid函数！

首先，给书sigmoid函数的定义：

然后我们看一下sigmoid函数图像：

a.通过图像可以看出当x在0点附近的时候，函数图像比较陡峭，而在横坐标轴两端处，即x>>0或者x<<0的时候，图像比较平稳。这也就说函数对于0值附近的数据比较敏感，对于边缘数据不敏感。这满足了我们的第一个需求。

b. 通过sigmoid函数的性质，可以知道：其定义域在全体实数，值域在[0,1]之间，并且当x>0的时候f(x)>0.5(正类);当x=0的时候f(x)=0.5（中性）；当x<0的时候f(x)<0.5(负类)。这满足了我们的第二个需求。就它了！

如何将线性函数转变为sigmoid函数

设线性函数为f(x)，令P(x)=P(Y=1|X=x)：具有特征x的样本被分到类别1的概率，则，被定义为让步比(odds ratio)。则就为我们想要线性函数f(x)达到我们想要的概率结果。

例如，当P(x)>0.5的时候，P(x)/[1-P(x)]>1,则f(x)>0,反过来也就是说，只要我们的线性函数满足f（x）>0就会有P(x)>0.5,此时X属于正类。同理，当f(x)<0的时候就会有P(x)<0.5,此时x属于负类，和上一部分讲解的sigmoid函数性质完全一致。

下面 我们就进行转换：

设：

我们把P(x)解出来：

则P(x)就是我们想要的一种概率模型，也就是我们所说的逻辑斯蒂回归模型，在此也就说明了逻辑斯蒂回归模型是一种概率模型。

python实现

#算法一 调用sklearn里面的方法
# _*_ encoding:utf-8 _*_

from matplotlib import pyplot
import scipy as sp
import numpy as np
from matplotlib import pylab
from sklearn.cross_validation import train_test_split
from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.metrics import precision_recall_curve
from sklearn.metrics import classification_report
from sklearn.linear_model import LogisticRegression
import time

start_time=time.time()

#绘制R/p曲线
def plot_pr(auc_score,precision,recall,label=None):
    pylab.figure(num=None, figsize=(6, 5))
    pylab.xlim([0.0, 1.0])
    pylab.ylim([0.0, 1.0])
    pylab.xlabel('Recall')
    pylab.ylabel('Precision')
    pylab.title('P/R (AUC=%0.2f) / %s' % (auc_score, label))
    pylab.fill_between(recall, precision, alpha=0.5)
    pylab.grid(True, linestyle='-', color='0.75')
    pylab.plot(recall, precision, lw=1)
    pylab.show()

#读取
movie_data=sp.load('move_data.npy')
movie_target=sp.load('move_target.npy')
x=movie_data
y=movie_target

#向量空间模型，
count_vec = TfidfVectorizer(binary = False, decode_error = 'ignore',\
                            stop_words = 'english')
average = 0
testNum = 10
#注意，训练样本调用fit_transform接口，测试样本调用的是transform接口
for i in range(0,testNum):
    #加载数据集，切分数据集80%训练，20%测试
    x_train, x_test, y_train, y_test\
        = train_test_split(movie_data, movie_target, test_size = 0.2)
    x_train = count_vec.fit_transform(x_train)
    x_test  = count_vec.transform(x_test)
    #训练LR分类器
    clf=LogisticRegression()
    clf.fit(x_train,y_train)
    y_pred=clf.predict(x_test)
    p=np.mean(y_pred==y_test)
    print '第%d次测试的准确率为：%.5f'%(i,p)
    average+=p


#精确率与召回率
#answer = clf.predict_proba(x_test)[:,0]#属于第neg类的概率
answer = clf.predict_proba(x_test)[:,1]#属于第pos类的概率

#这个地方可以查看官方文档，
precision, recall, thresholds = precision_recall_curve(y_test, answer)
#thresholds对应着一个向量，个数为answer中不重复的值的个数，每个位置上的值对应着一个threshold,并且是越来越大
#precision 和recall也是一个向量，个数比thresholds多一个，precision最后一个为1,recall最后一个值为0，确保图像从x轴=0开始
#precision每个值对应着，当阈值为thresholds对应位置值的时候的准确度
#例如：当取值取thresholds[0]的时候，precision[0]表示所有的训练集所对应的准确率recall[0]表示召回率

report = answer > 0.5#这里的0.5可以换为其他值
#我们通常使用0.5来划分两类数据，但是我们可以根据P/R图分析，选择一个合适的优秀的阈值。
#print report
print(classification_report(y_test, report, target_names = ['neg', 'pos']))
print("平均精确度为:", average/testNum)
print("花费的时间为:", time.time() - start_time)

plot_pr(0.6, precision, recall, "pos")

实验结果：

第0次测试的准确率为：0.81071
第1次测试的准确率为：0.79643
第2次测试的准确率为：0.79643
第3次测试的准确率为：0.80000
第4次测试的准确率为：0.80714
第5次测试的准确率为：0.81071
第6次测试的准确率为：0.75714
第7次测试的准确率为：0.78214
第8次测试的准确率为：0.78214
第9次测试的准确率为：0.81429

             precision    recall  f1-score   support

        neg       0.84      0.79      0.81       145
        pos       0.79      0.84      0.81       135

avg / total       0.82      0.81      0.81       280

平均精度率为: 0.795714285714
花费时间为: 12.1490001678

通过上图可以看出，如果选择的阈值过低，那么更多的测试样本都将分为1类，因此召回率得到提升，但是要牺牲相应的准确率。
注意precision_recall_curve（）方法中的thresholds中的阈值是逐渐增大的，对应到图像是就会，x轴从左到右对应的阈值是逐渐减小的。

# -*- encoding:utf-8 -*-
import matplotlib.pyplot as plt
from matplotlib import animation
import numpy as np


#加载数据集
def loadDataSet():
    '''
    :return: 输入向量矩阵和输出向量
    '''
    dataMat=[]
    labelMat=[]
    fr=open('test.txt','r')#一共有三列，每一列为x1,x2,Y
    for line in fr.readlines():
        lineArr=line.strip().split('##')
        #注意在这里将字符串转换成float型的时候总是遇到问题
        #感觉每一行的开头存在其他字符，所以我们多添加了一列“##”分隔符，是的第一个有效数字从第二个开始
        dataMat.append([1.0,float(lineArr[1]),float(lineArr[2])])#设x0为1，构成拓展之后的输入向量
        labelMat.append(int(lineArr[3]))
    return dataMat,labelMat

#可视化数据，画出数据集合逻辑斯蒂最佳回归直线
def plotBestFit(weights):
    dataMat,labelMat=loadDataSet()
    dataArr=np.array(dataMat)
    n=dataArr.shape[0]
    x1=[]
    y1=[]
    x2=[]
    y2=[]
    for i in range(n):
        if int(labelMat[i])==1:
            x1.append(dataArr[i,1])
            y1.append(dataArr[i,2])
        else:
            x2.append(dataArr[i,1])
            y2.append(dataArr[i,2])

    fig=plt.figure()
    ax=fig.add_subplot(111)
    ax.scatter(x1,y1,s=30,c='red',marker='s')
    ax.scatter(x2,y2,s=30,c='green',)
    if weights is not None:
        x=np.arange(-3.0,3.0,0.1)
        y=(-weights[0]-weights[1]*x)/weights[2]# #令w0*x0 + w1*x1 + w2*x2 = 0，其中x0=1，解出x1和x2的关系
        ax.plot(x,y)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

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


#逻辑斯提回归梯度上升批量算法,静态展示
def gradAscent_static(dataMatIn,classLabels):
    '''
    :param dataMatIn: 输入X矩阵，每一行代表一个实例，每一列分别是x0,x1,x2
    :param classLabels: 类别标签组成的向量
    :return: 权值向量
    '''
    dataMatrix=np.mat(dataMatIn)#转换为numpy矩阵数据类型,(100,3)
    labelMat=np.mat(classLabels).transpose()#转换为numpy矩阵数据类型,(100,1)
    m,n=dataMatrix.shape
    alpha=0.001
    maxCycles=500
    weights=np.ones((n,1))#(3,1)
    for k in range(maxCycles):
        h=sigmoid(dataMatrix*weights)
        error=(labelMat-h)#向量减法
        weights+=alpha*dataMatrix.transpose()*error#矩阵内积

    return weights

def draw_line(weights,lines):
    x = np.arange(-5.0, 5.0, 0.1)
    y = (-weights[0]-weights[1]*x)/weights[2]   #令w0*x0 + w1*x1 + w2*x2 = 0，其中x0=1，解出x1和x2的关系
    lines.set_data(x, y)
    return lines


def main1():#实现批量梯度上升，展示静态图像
    datas,labels=loadDataSet()
    weights=gradAscent_static(datas,labels)
    plotBestFit(weights)#静态图
if __name__ == '__main__':
    main1()

实验结果：

最终得到的权重向量为： [[ 4.12414349]
                    [ 0.48007329]
                    [-0.6168482 ]]

所谓的不平凡就是平凡的N次幂
                     ---------By Ada