Logistic回归多分类之鸢尾花

来源：互联网发布：新倩女幽魂网游mac版编辑：程序博客网时间：2024/04/25 19:43

作者：金良（golden1314521@gmail.com） csdn博客： http://blog.csdn.net/u012176591

- 多项逻辑回归模型原理
- 鸢尾花数据可视化
- 算法实现代码
- 混淆矩阵
- 进一步封装

1.多项逻辑回归模型原理

逻辑回归模型是二分类模型，用于二分类问题。可以将其推广为多项逻辑回归模型（multi-nominal logistic regression model），用于多分类。

假设类别Y的取值集合为{1,2,⋯,K}，那么多项逻辑回归模型是

P (y = k | x) = exp ( w k \cdot x ) 1 + \sum K - 1 k = 1 exp ( w k \cdot x ) k = 1, 2, \dots, K - 1 P (y = K | x) = 1 1 + \sum K - 1 k = 1 exp ( w k \cdot x )

似然函数为

\prod i = 1 N \prod k = 1 K P (y i = k | x i) y k i

其中

P(yi=k|xi)为模型在输入样本

xi时将其判为类别

k的概率，

yki起到指示函数的作用，当

k等于样本

xi的标签类别时为1，其余均为0.

对似然函数取对数，然后取负，得到L(w1,w2,⋯,wK−1)(简记为L(w))，最终要训练出的模型参数w1,w2,⋯,wK−1要使得L(w)的值取最小。

$L (w) = - 1 N \sum i = 1 N \sum k = 1 K y k i log P (y i = k | x i) = - 1 N \sum i = 1 N (\sum k = 1 K - 1 y k i log exp ( w k \cdot x i ) 1 + \sum K - 1 j = 1 exp ( w j \cdot x i ) + y K i log 1 1 + \sum K - 1 j = 1 exp ( w j \cdot x i )) = - 1 N \sum i = 1 N (\sum k = 1 K - 1 y k i (w k \cdot x i - log (1 + \sum j = 1 K - 1 exp (w j \cdot x i))) - y K i log (1 + \sum j = 1 K - 1 exp (w j \cdot x i))) = 1 N \sum i = 1 N (\sum k = 1 K y k i log (1 + \sum j = 1 K - 1 exp (w j \cdot x i)) - \sum k = 1 K - 1 y k i w k \cdot x i)$

考虑到过拟合的发生，对L(w)加上一个正则化项

∑K−1k=1|wk|22σ2

则L(w)写成

$L (w) = 1 N \sum i = 1 N (\sum k = 1 K y k i log (1 + \sum j = 1 K - 1 exp (w j \cdot x i)) - \sum k = 1 K - 1 y k i w k \cdot x i) + \sum k = 1 K - 1 | w k | 2 2 σ 2$

对L(w)关于wk求梯度，得到

$\partial L ( w ) w k = 1 N \sum i = 1 N exp ( w k \cdot x i ) 1 + \sum K - 1 k = 1 exp ( w k \cdot x i ) \cdot x i - 1 N \sum i = 1 N y k i x i + w k σ 2 = 1 N \sum i = 1 N P (y i = k | x i) \cdot x i - 1 N \sum i = 1 N y k i x i + w k σ 2$

第一项1N∑Ni=1P(yi=k|xi)⋅xi 可视为类别k的的后验期望值，第二项1N∑Ni=1ykixi 可视为类别k的先验期望值，第三项是正则化项，缓解过拟合。

接下来用梯度下降法对参数 w1,w2,⋯,wK−1 进行修正即可。

2.鸢尾花数据可视化

这里写图片描述

作图代码如下：

import numpy as npfrom mpl_toolkits.mplot3d import Axes3Dimport matplotlib.pyplot as pltimport itertools as itimport matplotlib as mpl#mpl.rcParams['xtick.labelsize'] = 6 #mpl.rcParams['ytick.labelsize'] = 6mpl.rcParams['axes.labelsize'] = 6 #设置所有坐标的刻度的字体大小attributes = ['SepalLength','SepalWidth','PetalLength','PetalWidth'] #鸢尾花的四个属性名comb = list(it.combinations([0,1,2,3],3)) # 每幅图挑出三个属性datas = []labels = []with open('Iris.txt','r') as f:     for line in f:         linedata =  line.split(' ')        datas.append(linedata[:-1]) #前4列是4个属性的值        labels.append(linedata[-1].replace('\n',''))  #最后一列是类别datas = np.array(datas)datas = datas.astype(float) #将二维的字符串数组转换成浮点数数组labels = np.array(labels)kinds = list(set(labels)) #3个类别的名字列表fig = plt.figure()plt.title(u'鸢尾花数据的可视化显示',{'fontname':'STFangsong','fontsize':10})for k in range(4):    ax = fig.add_subplot(int('22'+str(k+1)), projection='3d')    for label,color,marker in zip(kinds,['r','b','k','y','c','g'],['s','*','o','^']):        data = datas[labels == label ]        for i in range(data.shape[0]):            if i == 0: #对第一个点标记一个label属性，用于图例的显示，注意label不能设置多次，否则会重复显示                ax.plot([data[i,comb[k][0]]], [data[i,comb[k][1]]], [data[i,comb[k][2]]], color=color,marker = marker,markersize=3,label = label,linestyle = 'None')             else:                ax.plot([data[i,comb[k][0]]], [data[i,comb[k][1]]], [data[i,comb[k][2]]], color=color,marker = marker,markersize=3,linestyle = 'None')    ax.legend(loc="upper left",prop={'size':8},numpoints=1) #图例的放置位置和字体大小    ax.set_xlabel(attributes[comb[k][0]],fontsize=8) #设置坐标的标签为属性名    ax.set_ylabel(attributes[comb[k][1]],fontsize=8)    ax.set_zlabel(attributes[comb[k][2]],fontsize=8)    ax.set_xticks(np.linspace(np.min(datas[:,comb[k][0]]),np.max(datas[:,comb[k][0]]),3)) #设置坐标的刻度    ax.set_yticks(np.linspace(np.min(datas[:,comb[k][1]]),np.max(datas[:,comb[k][1]]),3))       ax.set_zticks(np.linspace(np.min(datas[:,comb[k][2]]),np.max(datas[:,comb[k][2]]),3)) plt.subplots_adjust(wspace = 0.1,hspace=0.1) #子图之间的空间大小savefig('Iris.png',dpi=700,bbox_inches='tight')plt.show()

3.算法实现代码

attributes = ['SepalLength','SepalWidth','PetalLength','PetalWidth'] #鸢尾花的四个属性名datas = []labels = []with open('Iris.txt','r') as f:     for line in f:         linedata =  line.split(' ')        datas.append(linedata[:-1]) #前4列是4个属性的值        labels.append(linedata[-1].replace('\n',''))  #最后一列是类别datas = np.array(datas)datas = datas.astype(float) #将二维的字符串数组转换成浮点数数组labels = np.array(labels)kinds = list(set(labels)) #3个类别的名字列表means = datas.mean(axis=0) #各个属性的均值stds = datas.std(axis=0) #各个属性的标准差N,M =  datas.shape[0],datas.shape[1]+1 #N是样本数，M是参数向量的维K = 3 #K是类别数data = np.ones((N,M))data[:,1:] = (datas - means)/stds   #对原始数据进行标准差归一化W = np.zeros((K-1,M)) # 存储参数矩阵priorEs = np.array([1.0/N*np.sum(data[labels == kinds[i]],axis=0) for i in range(K-1)]) #各个属性的先验先验期望值liklist = []for it in range(1000):    lik = 0 # 当前的对数似然函数值    for k in range(K-1):#似然函数值的第一部分        lik -= np.sum(np.dot(W[k],data[labels == kinds[k]].transpose()))    lik += 1.0/N*np.sum(np.log(np.sum(np.exp(np.dot(W,data.transpose())),axis = 0) +1)) #似然函数的第二部分    liklist.append(lik) #    wx = np.exp(np.dot(W,data.transpose()))    probs = np.divide(wx,1+np.sum(wx,axis=0).transpose()) #K-1*N的矩阵    posteriorEs = 1.0/N*np.dot(probs,data) #各个属性的后验期望值    gradients = posteriorEs - priorEs + 1.0/100 *W #梯度，最后一项是高斯项，防止过拟合    W  -= gradients #对参数进行修正#probM每行三个元素，分别表示data中对应样本被判给三个类别的概率probM = np.ones((N,K))probM[:,:-1] = np.exp(np.dot(data,W.transpose()))probM /= np.array([np.sum(probM,axis = 1)]).transpose() #得到概率predict = np.argmax(probM,axis = 1).astype(int)   #取最大概率对应的类别#rights列表储存代表原始标签数据的序号，根据labels数据生成rights = np.zeros(N) rights[labels == kinds[1]] =1rights[labels == kinds[2]] =2rights = rights.astype(int)print '误判的样本的个数：%d\n'%np.sum(predict != rights) #误判的个数

该算法的收敛性如下图所示：
这里写图片描述

4.混淆矩阵

#输出混淆矩阵K = np.shape(list(set(rights).union(set(predict))))[0]matrix = np.zeros((K,K))for righ,predic in zip(rights,predict):    matrix[righ,predic] += 1print "%-12s" % ("类别"),for i in range(K):    print "%-12s" % (kinds[i]),  print '\n',for i in range(K):    print "%-12s" % (kinds[i]),     for j in range(K):        print "%-12d" % (matrix[i][j]),     print '\n',

这里写图片描述

上述混淆矩阵（confusion matrix）的行表示真实类别，列代表模型输出类别。可知，类setosa 的样本没有误判为其它类别，而且其它类别的样本没有误判为类 setosa 的，也就是说，我们训练出的模型没有将类setosa 的样本与其它类的样本混淆；4个versicolor 样本被误判为’virginica’ ， 4个’virginica’样本被误判为versicolor 。总体看来，分类效果还是不错的。

5.进一步封装

def LogisticRegression(datas,labels):    kinds = list(set(labels)) #3个类别的名字列表    means = datas.mean(axis=0) #各个属性的均值    stds = datas.std(axis=0) #各个属性的标准差    N,M =  datas.shape[0],datas.shape[1]+1 #N是样本数，M是参数向量的维    K = np.shape(list(set(labels)))[0] #K是类别数    data = np.ones((N,M))    data[:,1:] = (datas - means)/stds   #对原始数据进行标准差归一化       W = np.zeros((K-1,M)) # 存储参数矩阵    priorEs = np.array([1.0/N*np.sum(data[labels == kinds[i]],axis=0) for i in range(K-1)]) #各个属性的先验先验期望值    for it in range(1000):        wx = np.exp(np.dot(W,data.transpose()))        probs = np.divide(wx,1+np.sum(wx,axis=0).transpose()) #K-1*N的矩阵        posteriorEs = 1.0/N*np.dot(probs,data) #各个属性的后验期望值        gradients = posteriorEs - priorEs + 1.0/100 *W #梯度，最后一项是高斯项，防止过拟合        W  -= gradients #对参数进行修正    #probM每行三个元素，分别表示data中对应样本被判给三个类别的概率    probM = np.ones((N,K))    probM[:,:-1] = np.exp(np.dot(data,W.transpose()))    probM /= np.array([np.sum(probM,axis = 1)]).transpose() #得到概率    return probMif __name__ == "__main__":    datas = []    labels = []    with open('Iris.txt','r') as f:         for line in f:             linedata =  line.split(' ')            datas.append(linedata[:-1]) #前4列是4个属性的值            labels.append(linedata[-1].replace('\n',''))  #最后一列是类别    datas = np.array(datas)    datas = datas.astype(float) #将二维的字符串数组转换成浮点数数组    labels = np.array(labels)    kinds = list(set(labels)) #3个类别的名字列表    #LogisticRegression可以接受不同个数的类别数据，这里用元数据的子集（两个类别的样本）来测试    m,n = 1,2    sel = np.array(list(where(labels == kinds[m])[0])+(list(where(labels == kinds[n])[0])))    slabels = labels[sel]    sdatas = datas[sel]    probM = LogisticRegression(sdatas,slabels)    print probM

1 0