使用logistic regression对数据进行分类

来自集智百科
跳转到: 导航搜索

目录

logistic回归的原理

我们先制造一份示例数据来看看logistic 回归是如何工作的。

首先调用模块

    import numpy as np
    from sklearn import linear_model, datasets
    import matplotlib.pyplot as plt
    from scipy.stats import norm

然后生成模拟数据并制图

 
    np.random.seed(3)
    n = 40
    X = np.hstack((norm.rvs(loc=2, size=n, scale=2), norm.rvs(loc=8, size=n, scale=3)))
    y = np.hstack((np.zeros(n),np.ones(n)))
    plt.figure(figsize=(10, 4))
    plt.xlim((-5, 20))
    plt.scatter(X, y, c=y)
    plt.xlabel("feature value")
    plt.ylabel("class")
    plt.grid(True, linestyle='-', color='0.75')
    plt.savefig("E:/wulingfei/logistic_classification/logistic_classify1.png", bbox_inches="tight")

我们就看到这样一张图

Logistic classify1.png

上图有80个数据点,在模拟X的时候,我们使用了两个正态分布,分别制定各自的均值,方差,生成40个点。模拟Y的时候我们直接生成40个0和40个1。两个正态分布均值虽然不同,但在方差作用下,在X上的分布有重叠。我们假设这是数据中某个feature的表现,要求使用一个模型,根据这个feature来判别所属的Y类别(class)。

这种情况比较适合使用logistic regression,为什么呢?

我们先来看一下线性回归的公式


y_i = \beta_0 +\beta_1 x_i

如果我们现在手头需要预测的因变量不是连续的,而是一个0/1变量,那怎么办呢?一种思路就是把公式进行变形,使得我们原来积累的线性回归的公式继续能用。

怎么做到这一点呢?我们可以考虑如下三个变形步骤:

1. 从离散变成连续:不是直接判定一个case是0还是1,而是判定它属于0或者1的概率;

2. 从0\sim 1之间的连续变成0\sim \infty的连续:不考虑0或者1的绝对概率,而是考虑取1的概率和取0的概率的比值odds ratio;

3. 从0\sim \infty的连续变成-\infty\sim \infty的连续:把odds ratio取对数。


    k = np.arange(0, 1, 0.001)
    plt.figure(figsize=(10, 5))
    plt.subplot(1, 2, 1)
    plt.xlim((0, 1)); plt.ylim((0, 10))
    plt.plot(k, k / (1 - k))
    plt.xlabel("P"); plt.ylabel("odds = P / (1-P)")
    plt.grid(True, linestyle='-', color='0.75')
    plt.subplot(1, 2, 2)
    plt.xlim((0, 1))
    plt.plot(k, np.log(k / (1 - k)))
    plt.xlabel("P"); plt.ylabel("log(odds) = log(P / (1-P))")
    plt.grid(True, linestyle='-', color='0.75')
    plt.tight_layout(pad=0.4, w_pad=0, h_pad=1.0)     
    plt.savefig("E:/wulingfei/logistic_classification/logistic_classify2.png", bbox_inches="tight")

通过上述代码,我们看到经过这三步变形,我们得到了一个理论值在负无穷到正无穷的因变量,如下图所示,这样,我们之前的线性回归的知识就可以用上了。

Logistic classify2.png

那么,原来线性回归,现在就变成了logistic回归:


log(\frac{p_i}{1-p_i}) = \beta_0 +\beta_1 x_i

或者也可以写成


p_i = \frac{1}{1+ e^{-(\beta_0 +\beta_1 x_i)}}

现在,我们手头已经有了模型,剩下的任务就是参数拟合。也就是说,我们要找到这样的的参数\beta_0\beta_1,使得所有的{x_i,p_i} pair根据这个模型拟合之后,总的误差最小。实际的参数拟合中,因为我们实际上只有一个p_i,所以往往采取最大似然估计的方法。

使用sklearn包自带的拟合函数,我们可以对示例数据做线性和logistic拟合,作为对比

 
    #---linear regression----------
    xs=np.linspace(-5,15,10)
    from sklearn.linear_model import LinearRegression
    clf = LinearRegression()
    clf.fit(X.reshape(n * 2, 1), y)
    def lin_model(clf, X):
        return clf.intercept_ + clf.coef_ * X
    #---logistic regression--------
    from sklearn.linear_model import LogisticRegression
    logclf = LogisticRegression()
    logclf.fit(X.reshape(n * 2, 1), y)
    def lr_model(clf, X):
        return 1.0 / (1.0 + np.exp(-(clf.intercept_ + clf.coef_ * X)))
    #----plot---------------------------    
    plt.figure(figsize=(10, 5))
    plt.subplot(1, 2, 1)
    plt.scatter(X, y, c=y)
    plt.plot(X, lin_model(clf, X),"o",color="orange")
    plt.plot(xs, lin_model(clf, xs),"-",color="green")
    plt.xlabel("feature value")
    plt.ylabel("class")
    plt.title("linear fit")
    plt.grid(True, linestyle='-', color='0.75')
    plt.subplot(1, 2, 2)
    plt.scatter(X, y, c=y)
    plt.plot(X, lr_model(logclf, X).ravel(),"o",color="c")
    plt.plot(xs, lr_model(logclf, xs).ravel(),"-",color="green")
    plt.xlabel("feature value")
    plt.ylabel("class")
    plt.title("logistic fit")
    plt.grid(True, linestyle='-', color='0.75')
    plt.tight_layout(pad=0.4, w_pad=0, h_pad=1.0)     
    plt.savefig("E:/wulingfei/logistic_classification/logistic_classify3.png", bbox_inches="tight")

Logistic classify3.png

其中连续曲线展示了我们的模型,黄色和绿色的数据点分别对应着线性模型和Logistic模型对原始数据点的预测。logistic的表现比较好,例如对于x略小于5左右的数据点,logistic预测属于class 1概率已经大于0.6,目测来说也是比较合理的。而我们使用线性回归判别的话(我们很粗暴地假设线性预测值>0.5则为class 1),此时仍然是被判别为class 0的。

使用logistic regression处理iris数据

Iris example.jpg


接下来,我们使用logistic来处理一个真实的数据,iris花的数据。关于这个数据的基本描述可以看这里

这张图展示了4维特征空间中的数据所有的二维投影:

Iris clustering 2.png

这里的情况比之前的模拟数据复杂,应为我们现在的自变量有两个而不是一个(我们从四个维度中选取两个),因变量有三个而不是两个。不过幸好sklearn的fit函数和predict函数让我们处理起这些问题来变得很容易。

使用logistic regression处理音乐数据

音乐数据训练样本的获得和使用快速傅里叶变换(FFT)预处理的方法见之前在这里已经介绍过。我们现在的任务是

1. 把训练集扩大到每类100个首歌而不是之前的10首歌,类别仍然是六类:jazz,classical,country, pop, rock, metal。

2. 同时使用logistic回归和KNN作为分类器。

3. 引入一些评价的标准来比较Logistic和KNN在测试集上的表现。

1.扩大样本规模

因为我们每次都要对音频进行FFT处理,这个速度是比较慢的,所以为了加快训练的速度,我们可以把处理完的数据先存出来,在训练分类器时直接读取这个保存的数值型数据。为了实现这一点,我们操作如下:

    def create_fft(g,n):
        rad="E:/wulingfei/music_classification/genres/"+g+"/converted/"+g+"."+str(n).zfill(5)+".au.wav"
        sample_rate, X = wavfile.read(rad)
        fft_features = abs(scipy.fft(X)[:1000])
        sad="E:/wulingfei/music_classification/trainset/"+g+"."+str(n).zfill(5)+ ".fft"
        np.save(sad, fft_features)
 
    #-------creat fft--------------
 
    genre_list = ["classical", "jazz", "country", "pop", "rock", "metal"]
    for g in genre_list:
        for n in range(100):
            create_fft(g,n)
 
    #-------read fft-------------- 
 
    X=[]
    Y=[]
    for g in genre_list:
        for n in range(100):
            rad="E:/wulingfei/music_classification/trainset/"+g+"."+str(n).zfill(5)+ ".fft"+".npy"
            fft_features = np.load(rad)
            X.append(fft_features)
            Y.append(genre_list.index(g))
 
    X=np.array(X)
    Y=np.array(Y)

在这个过程中,1.35G的原始音乐数据提出来的特征数据被存储成4.6M的矩阵数据,数据被压缩了3万倍,自然处理时间大大加快。


2.同时使用Logistic和KNN作为分类器

首先我们要将原始数据分为训练集和测试集,这里是随机抽样80%做测试集,剩下20%做训练集

    import random
    randomIndex=random.sample(range(len(Y)),int(len(Y)*8/10))
    trainX=[];trainY=[];testX=[];testY=[]
    for i in range(len(Y)):
        if i in randomIndex:
            trainX.append(X[i])
            trainY.append(Y[i])
        else:
            testX.append(X[i])
            testY.append(Y[i])

接下来,我们使用sklearn,来构造和训练我们的两种分类器

    #------train logistic classifier-------------- 
    from sklearn.linear_model import LogisticRegression
    logclf = LogisticRegression()
    logclf.fit(trainX, trainY)
    predictYlogistic=map(lambda x:logclf.predict(x)[0],testX)
 
    #----train knn classifier-----------------------
    from sklearn.neighbors import NearestNeighbors
    neigh = NearestNeighbors(n_neighbors=1)
    neigh.fit(trainX) 
    predictYknn=map(lambda x:trainY[neigh.kneighbors(x,return_distance=False)[0][0]],testX)

将predictYlogistic以及predictYknn与testY对比,我们就可以知道两者的判定正确率

    a = array(predictYlogistic)-array(testY)
    accuracyLogistic = 1-np.count_nonzero(a)/len(a)
    b = array(predictYknn)-array(testY)
    accuracyKNN = 1-np.count_nonzero(b)/len(b)

logistic的正确率是46%,KNN的正确率是63%。后者是更好的分类器。当然,更严格的检验要使用cross-validation,如下图所示

3.Logistic和KNN作为分类器的效果比较

Iris clustering 4.png

    def accuracyRate(classifier,trainX,trainY,testX,testY):
        if classifier=="logistic":
            logclf = LogisticRegression()
            logclf.fit(trainX, trainY)
            predictY=map(lambda x:logclf.predict(x)[0],testX)
        if classifier=="knn":
            neigh = NearestNeighbors(n_neighbors=1)
            neigh.fit(trainX) 
            predictY=map(lambda x:trainY[neigh.kneighbors(x,return_distance=False)[0][0]],testX)
        a = array(predictY)-array(testY)
        accuracy = 1-np.count_nonzero(a)/len(a)
        return accuracy
 
    index=range(len(Y))
    random.shuffle(index)
    ARS=[]
    for i in range(5):
        l1=range(i*100,(i+1)*100)
        l2=map(lambda x:index[x], l1)
        l3=np.setdiff1d(index,l2)
        testX=map(lambda k:X[k], l2)
        testY=map(lambda k:Y[k], l2)
        trainX=map(lambda k:X[k], l3)
        trainY=map(lambda k:Y[k], l3)
        r1=accuracyRate("logistic",trainX,trainY,testX,testY)
        r2=accuracyRate("knn",trainX,trainY,testX,testY)
        ARS.append([r1,r2])
 
    ARS=np.array(ARS)

得到的ARS就包含了两种分类器5次测试的结果。我们也可以把这个结果用柱形图画出来

    ind = np.arange(len(ARS))  # the x locations for the groups
    width = 0.3 
    fig = plt.figure()
    ax = fig.add_subplot(111)
    g1=ax.bar(ind, ARS[:,0], width, color='c',alpha=0.7)
    g2=ax.bar(ind+width, ARS[:,1], width, color='orange',alpha=0.7)
    ax.set_ylabel('Accuracy rate')
    ax.set_title('Accuracy rate of two classifiers')
    ax.set_xticks(ind+width)
    ax.set_xticklabels(('test1', 'test2', 'test3', 'test4', 'test5'))
    legend( (g1[0], g2[0]), ('Logistic', 'KNN'),loc="upper center" )
    plt.plot()
    plt.savefig("E:/wulingfei/logistic_classification/logistic_classify6.png", bbox_inches="tight")

Logistic classify6.png

计算得到,logistci分类器的平均准确率是41%,KNN分类器的平均正确率是60%。还是KNN表现比较好

但是,之前的这种cross-validation只能用于测量分类器的整体表现,无法具体分析分类器在某一个具体类别上的表现,不利于我们诊断bad case和提高feature engieering的水平,为了解决这个问题,我们引入confusion matrix来诊断分类器的问题。


 
    from sklearn.metrics import confusion_matrix
    cmlogistic = confusion_matrix(testY, predictYlogistic)
    cmknn = confusion_matrix(testY, predictYknn)
 
    def plotCM(cm,title,colorbarOn,givenAX):
        ncm=cm/cm.max()
        plt.matshow(ncm, fignum=False, cmap='Blues', vmin=0, vmax=1.0)
        if givenAX=="":
            ax=plt.axes()
        else:
            ax = givenAX
        ax.set_xticks(range(len(genre_list)))
        ax.set_xticklabels(genre_list)
        ax.xaxis.set_ticks_position("bottom")
        ax.set_yticks(range(len(genre_list)))
        ax.set_yticklabels(genre_list)
        plt.title(title,size=12)
        if colorbarOn=="on":
            plt.colorbar()
        plt.xlabel('Predicted class')
        plt.ylabel('True class')
        for i in range(cm.shape[0]):
            for j in range(cm.shape[1]):
                text(i,j,cm[i,j],size=15)
 
    plt.figure(figsize=(10, 5))  
    fig1=plt.subplot(1, 2, 1)          
    plotCM(cmlogistic,"confusion matrix: FFT based logistic classifier","off",fig1.axes)   
    fig2=plt.subplot(1, 2, 2)     
    plotCM(cmknn,"confusion matrix: FFT based KNN classifier","off",fig2.axes) 
    plt.tight_layout(pad=0.4, w_pad=0, h_pad=1.0)     
 
    plt.savefig("E:/wulingfei/logistic_classification/logistic_classify5.png", bbox_inches="tight")

Logistic classify5.png

上图中对角线代表正确诊断的比例。对角线方格的颜色越深,其他位置方格的颜色越浅,分类器效果越好。可以看出,logistic在pop和rock两类上的表现明显比起KNN分类器明显较弱。这意味着如果我们要继续使用logitic,应该在这两类音乐的feature提取上再下更多功夫。

个人工具
名字空间
操作
导航
工具箱