【ML学习笔记】17：多元正态分布下极大似然估计最小错误率贝叶斯决策

简述多元正态分布下的最小错误率贝叶斯

如果特征的值向量服从d元正态分布，即其概率密度函数为：

即其分布可以由均值向量和对称的协方差矩阵

唯一确定。

如果认为样本的特征向量在类内服从多元正态分布：

即对于每个类i，具有各自的类内的均值向量和协方差矩阵。

如之前所学，最小错误率贝叶斯的判别函数的原始形式是：

类条件概率密度服从多元正态分布，带入，得：

因为是比较大小用的，去掉与类号i无关的项：

而用于分类的决策面是：

在实际问题下的措施

①反向判别

这次要做的还是用[身高，体重，鞋码]->[性别]的数据，认为类别男女出现的次数一样，即先验概率都是0.5，则可以在判别函数中去除这一项：

判别函数仅仅是比较大小用的，为了减少计算机计算量，改用反向判别的函数：

反向判别函数则要取函数值小的那一方作为预测的类别。

②分类面采样以对ROC上采样点做估计

对于多元正态分布，用理论公式的方式(多重积分)去计算错误率是很困难的，所以绘制ROC曲线时，我采用改变分类面，获得多个采样分类面，在采样分类面上用频率估计错误率。

改变采样分类面的方式可以改变分类面的常数阈值，即把0改成一定区间上的正负值b：

代码实现

#-*-coding:utf-8-*-
from numpy import *
import operator
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D #绘制3D坐标的函数#求二元正态分布概率密度函数值
#传入的向量x=[特征i值,特征j值]
def gauss2(x,miu,sgm):r=mat([x[0]-miu[0],x[1]-miu[1]])multi=r*sgm.I*r.Tmulti=float(multi) #1乘1矩阵取内容k=exp(-multi/2) #.I求逆,.T转置k/=2*math.pi*linalg.det(sgm)**(1/2) #linalg.det求行列式的值return k#从文件中读取数据
def getData():fr1=open(r'boynew.txt')fr0=open(r'girlnew.txt')arrayOLines1=fr1.readlines() #读取文件arrayOLines0=fr0.readlines()#特征矩阵dataMat1=[]dataMat0=[]#男同学for line in arrayOLines1:line=line.strip() #strip()去掉首尾空格listFromLine=line.split() #按空白字符分割成列表#string变floatfor i in range(len(listFromLine)):listFromLine[i]=float(listFromLine[i])#加入特征矩阵dataMat1.append(listFromLine)#女同学for line in arrayOLines0:line=line.strip() #strip()去掉首尾空格listFromLine=line.split() #按空白字符分割成列表#string变floatfor i in range(len(listFromLine)):listFromLine[i]=float(listFromLine[i])#加入特征矩阵dataMat0.append(listFromLine)return dataMat1,dataMat0#训练(不带归一化)
def TraBys(i,j):dataMat1,dataMat0=getData() #获取数据dataMat1=mat(dataMat1)dataMat0=mat(dataMat0)#求男女同学的均值和协方差矩阵miu1,miu0,sgm1,sgm0=getMiuSgm(i,j,dataMat1,dataMat0)return miu1,miu0,sgm1,sgm0#训练(带归一化)
def TraBysStd(i,j):dataMat1,dataMat0=getData() #获取数据'''#生成标签向量labelVec=[1 for k in range(len(dataMat1))]labelVec0=[0 for k in range(len(dataMat0))]labelVec.extend(labelVec0)'''#男女特征集放到一起dataMat=dataMat1dataMat.extend(dataMat0)#并转换为矩阵对象dataMat,保留了dataMat1和dataMat0!dataMat=mat(dataMat)dataMat1=mat(dataMat1)dataMat0=mat(dataMat0)#用合矩阵dataMat0返回归一化之后的矩阵#主要是要获得原范围ranges,原最小值minVals#这两个参数使用单一的dataMat1和dataMat0都没法表征总体!newDataMat,ranges,minVals=autoNorm(dataMat)#还要把男女同学自己的矩阵归一化#才能传进去求均值和协方差矩阵newDataMat1,_,_=autoNorm(dataMat1)newDataMat0,_,_=autoNorm(dataMat0)#维度修正,实在是无奈之举newDataMat1=reshape(newDataMat1,(len(dataMat1),-1))newDataMat0=reshape(newDataMat0,(len(dataMat0),-1))#求男女同学的均值和协方差矩阵miu1,miu0,sgm1,sgm0=getMiuSgm(i,j,newDataMat1,newDataMat0)#除了返回这四个参数,还要返回原范围和最小值#这是为了对外来特征向量归一化return miu1,miu0,sgm1,sgm0,ranges,minVals#对外来特征向量x归一化
def rebuidx(x,ranges,minVals):return (x-minVals)/ranges#获取男女同学的均值和协方差矩阵
def getMiuSgm(i,j,dataMat1,dataMat0):#计算男生i特征均值Imiu1=0.0for r in range(len(dataMat1)):Imiu1+=dataMat1[r,i] #注意矩阵不能像数组那样取值!Imiu1=Imiu1/len(dataMat1)#计算女生i特征均值Imiu0=0.0for r in range(len(dataMat0)):Imiu0+=dataMat0[r,i]Imiu0=Imiu0/len(dataMat0)#计算男生j特征均值Jmiu1=0.0for r in range(len(dataMat1)):Jmiu1+=dataMat1[r,j]Jmiu1=Jmiu1/len(dataMat1)#计算女生j特征均值Jmiu0=0.0for r in range(len(dataMat0)):Jmiu0+=dataMat0[r,j]Jmiu0=Jmiu0/len(dataMat0)#男生的均值向量miu1=[Imiu1,Jmiu1]#女生的均值向量miu0=[Imiu0,Jmiu0]#求男生的协方差矩阵sgm1=mat(zeros((2,2)))for k in range(len(dataMat1)):x=mat([dataMat1[k,i]-miu1[0],dataMat1[k,j]-miu1[1]])sgm1+=x.T*xsgm1/=len(dataMat1)#求女生的协方差矩阵sgm0=mat(zeros((2,2)))for k in range(len(dataMat0)):x=mat([dataMat0[k,i]-miu0[0],dataMat0[k,j]-miu0[1]])sgm0+=x.T*xsgm0/=len(dataMat0)return miu1,miu0,sgm1,sgm0 #返回两均值,两方差#绘制两个类条件概率密度图
def show(i,j,miu1,miu0,sgm1,sgm0):'''#类似于标准差的东西用来确定采样点起始和终止范围sqrtsgm1=linalg.det(sgm1)**(1/2) #行列式的(1/2)次幂,类似于标准差sqrtsgm0=linalg.det(sgm0)**(1/2) #行列式的(1/2)次幂,类似于标准差'''#采样间距if i==0: #如果特征i是身高,那么采样区间是140~200Ilft=140Irgt=200xlab='Height'elif i==1: #如果特征i是体重,那么采样区间是30~80Ilft=30Irgt=80xlab='Weight'else: #如果特征i是鞋码,那么采样区间是32~46Ilft=32Irgt=46xlab='Shoe Size'if j==0:Jlft=140Jrgt=200ylab='Height'elif j==1:Jlft=30Jrgt=80ylab='Weight'else:Jlft=32Jrgt=46ylab='Shoe Size'#linspace创建等差数列,在I/J方向上取采样点I1=linspace(Ilft,Irgt,50)J1=linspace(Jlft,Jrgt,50)I0=linspace(Ilft,Irgt,50)J0=linspace(Jlft,Jrgt,50)#映射以扩充为二维采样面上的点X1,Y1=meshgrid(I1,J1)X0,Y0=meshgrid(I0,J0)#用正态分布概率密度函数得到采样点的Z值序列Z1=[]for k in range(len(J1)):for p in range(len(I1)):Z1.append(gauss2([X1[k][p],Y1[k][p]],miu1,sgm1))Z1=reshape(Z1,(len(J1),len(I1)))Z0=[]for k in range(len(J0)):for p in range(len(I0)):Z0.append(gauss2([X0[k][p],Y0[k][p]],miu0,sgm0))Z0=reshape(Z0,(len(J0),len(X0)))#用这些坐标点绘制图像fig1=plt.figure()#创建一个绘图对象ax=Axes3D(fig1)#用这个绘图对象创建一个Axes对象(有3D坐标)#绘制男生i,j特征的类条件概率密度函数ax.plot_surface(array(X1),array(Y1),Z1,rstride=1,cstride=1,cmap=plt.cm.coolwarm)#绘制女生i,j的类条件概率密度函数ax.plot_surface(array(X0),array(Y0),Z0,rstride=1,cstride=1,cmap=plt.cm.jet)plt.title("The probability density function of the class condition")#给三个坐标轴注明ax.set_xlabel(xlab,color='r')  ax.set_ylabel(ylab,color='g')  ax.set_zlabel('Probability Density',color='b')plt.show()#不同的特征,范围可能不同,作差得到的值可能差别很大
#如果它们的重要程度一样,显然不应如此
#可以将不同取值范围的特征值数值归一化到0~1之间
def autoNorm(dataMat):#下面参数0表示从列中取最大最小minVals=dataMat.min(0) #每列的最小值存到minVals里maxVals=dataMat.max(0) #每列的最大值存到maxVals里ranges=maxVals-minVals #每列最大最小值之差,存到ranges里newDataMat=zeros(shape(dataMat)) #用来存归一化后的矩阵m=dataMat.shape[0] #取第0维即行数newDataMat=dataMat-tile(minVals,(m,1)) #把最小值重复m成m行,用原值减去newDataMat=newDataMat/tile(ranges,(m,1)) #把减完的每个偏差值除以自己列最大和最小的差值return newDataMat,ranges,minVals #返回归一化之后的矩阵,原范围,原最小值#计算dxd矩阵(二次项携带矩阵)
def Wi(sgm):return (-1/2)*sgm.I#计算d维列向量(一次项系数矩阵)
def wi(miu,sgm):return sgm.I*miu.T#计算常数项,还要传入先验概率
def omgi(miu,sgm,Ppre):return float((-1/2)*miu*sgm.I*miu.T)-\(1/2)*log(linalg.det(sgm))+\log(Ppre)#新的计算判别函数的方式
#因为先验概率都是0.5,不妨去掉这一项,然后乘以系数-2作反向判别
def newgx(miu1,miu0,sgm1,sgm0,x):#暂存第一项s1=(x-miu1)*sgm1.I*(x-miu1).Ts0=(x-miu0)*sgm0.I*(x-miu0).T#取矩阵内容s1=float(s1)s0=float(s0)#求出两个反向判别gxrvs_g1x=s1+log(linalg.det(sgm1))rvs_g0x=s0+log(linalg.det(sgm0))#因为反向判别,所以反向返回return rvs_g0x,rvs_g1x#这个函数存在问题,但是没有发现问题所在之处
#用了上面的newgx去代替它,效果不错
#传入训练好的参数(两均值,两标准差),以及投入测试的特征向量x(mat化后)
#计算两类的判别函数值g1(x)和g0(x)
def gx(miu1,miu0,sgm1,sgm0,x):  #二次项携带矩阵W1=Wi(sgm1)W0=Wi(sgm0)#一次项系数矩阵w1=wi(miu1,sgm1)w0=wi(miu0,sgm0)#常数项omg1=omgi(miu1,sgm1,0.5) #认为男女先验概率都是0.5omg0=omgi(miu0,sgm0,0.5)#二次/一次项值1x1方阵s1_2=x*W1*x.Ts1_1=w1.T*x.Ts0_2=x*W0*x.Ts0_1=w0.T*x.T#提取方阵值s1_2=float(s1_2)s1_1=float(s1_1)s0_2=float(s0_2)s0_1=float(s0_1)#计算判别函数值g1x=s1_2+s1_1+omg1g0x=s0_2+s0_1+omg0#返回求得的判别函数值return g1x,g0x#对归一化后的x向量作判别,决策面方程g1(x)-g0(x)=b
#默认b=0,之所以设定b可改变是为了绘制ROC曲线
def TstBys(miu1,miu0,sgm1,sgm0,x,b=0):#列表变矩阵,否则不能转置x=mat(x)miu1=mat(miu1)miu0=mat(miu0)#计算判别函数值#g1x,g0x=gx(miu1,miu0,sgm1,sgm0,x)g1x,g0x=newgx(miu1,miu0,sgm1,sgm0,x)#返回判别类的标号if g1x-g0x > b:return 1else:return 0#用测试集做错误率测试,传入特征号i和j
def errTst(i,j):fr1=open(r'boytst.txt') #测试集fr0=open(r'girltst.txt')arrayOLines1=fr1.readlines() #读取文件arrayOLines0=fr0.readlines()#特征矩阵dataMat=[]#标签向量labelVec=[]#男同学for line in arrayOLines1:line=line.strip() #strip()去掉首尾空格listFromLine=line.split() #按空白字符分割成列表#string变floatfor r in range(len(listFromLine)):listFromLine[r]=float(listFromLine[r])#加入特征矩阵dataMat.append(listFromLine)#加入标签向量labelVec.append(1)#女同学for line in arrayOLines0:line=line.strip() #strip()去掉首尾空格listFromLine=line.split() #按空白字符分割成列表#string变floatfor r in range(len(listFromLine)):listFromLine[r]=float(listFromLine[r])#加入特征矩阵dataMat.append(listFromLine)#加入标签向量labelVec.append(0)errCount=0.0 #错误分类次数计数#极大似然法训练出四个参数(两均值,两协方差矩阵)#miu1,miu0,sgm1,sgm0,minVals,ranges=TraBysStd(i,j)miu1,miu0,sgm1,sgm0=TraBys(i,j)#对于测试集的每个样本(下标)for k in range(len(labelVec)):#特征矩阵第k行的第[i,j]号特征向量做测试''' x=dataMat[k] #先把第k行全取过来x=rebuidx(x,minVals,ranges) #对x作归一化,这时不是列表了!x=[x[0,i],x[0,j]] #然后再变成只取i号和j号特征值的'''x=[dataMat[k][i],dataMat[k][j]]#如果预测的类别号和真实的标签向量中存的不一样if TstBys(miu1,miu0,sgm1,sgm0,x,0)!=labelVec[k]:errCount+=1 #记录分类错误print len(labelVec),"次测试错误率是:",errCount/len(labelVec)return dataMat,labelVec #返回测试集,用来求后面的两类错误率#对于传入的b值,改变超平面的位置
#从而获取假阳性率(误识率)和召回率(1-拒识率)
#假设1男性为阳性,0女性为阴性
#用测试集做错误率测试,传入特征号i和j
def getROCxy(dataMat,labelVec,i,j,b):FPR=0.0 #假阳性数->率RECALL=0.0 #召回数->率FPplusTN=0.0 #存阴性样本数TPplusFN=0.0 #存阳性样本数#极大似然法训练出四个参数(两均值,两协方差矩阵)miu1,miu0,sgm1,sgm0=TraBys(i,j)#对于测试集的每个样本(下标)for k in range(len(labelVec)):#特征矩阵第k行的第[i,j]号特征向量做测试x=[float(dataMat[k][i]),float(dataMat[k][j])]#本次预测出来的类号yc=TstBys(miu1,miu0,sgm1,sgm0,x,b)#如果实际是阳性样本if labelVec[k]==1:TPplusFN+=1 #阳性样本数+1else: #如果是阴性样本FPplusTN+=1 #阴性样本数+1#阳性样本 && 预测为阳性if labelVec[k]==1 and yc==1:RECALL+=1 #真阳性数+1#阴性样本 && 预测为阳性elif labelVec[k]==0 and yc==1:  FPR+=1 #假阳性数+1#print FPR,FPplusTN,RECALL,TPplusFN#假阳性率FPR/=FPplusTN#召回(真阳性)率RECALL/=TPplusFNreturn FPR,RECALL#演示用的函数,传入两个特征编号i和j
def Go(i,j):miu1,miu0,sgm1,sgm0=TraBys(i,j) #训练获得4个参数show(i,j,miu1,miu0,sgm1,sgm0) #绘制类条件概率密度图dataMat,labelVec=errTst(i,j) #求测试错误率print '正在采样绘制ROC...'#绘制ROC曲线用的采样点X=[]Y=[]#分类面采样阈值从-50到60,每次移动1for b in range(-50,60,1):#获取假阳率和真阳率FPR,RECALL=getROCxy(dataMat,labelVec,i,j,b)#分别加入到坐标采样列表中X.append(FPR)Y.append(RECALL)#绘制ROC曲线plt.plot(X,Y,"g-",linewidth=2)#横纵坐标名称,标题名称plt.xlabel('FPR')plt.ylabel('RECALL')if i==0:xlab='Height'elif i==1:xlab='Weight'else:xlab='Shoe Size'if j==0:ylab='Height'elif j==1:ylab='Weight'else:ylab='Shoe Size'plt.title('ROC curve (if use '+xlab+' and '+ylab+' as feature)')plt.grid(True) #显示网格plt.show() #显示

测试

直接用封装好的演示用的函数演示功能：

import bayes as bs
bs.Go(1,2)

类条件概率密度函数图：

关闭窗口后开始采样绘制ROC：

绘制好的ROC曲线：

与一元情况的比较

在上篇中，仅使用身高或者仅使用体重做判别，得到的分类器效果都很不好。从常理上也能解释，身高或者体重都不是鉴别男女的有效手段。建立好这个分类器后，可以看看这两个特征联合的分类效果如何。

看似不那么合适的两个特征，联合后的分类效果好了很多。