快捷導(dǎo)航

python實(shí)現(xiàn)高斯判別分析算法的例子

更新時(shí)間：2019年12月09日 08:43:51 作者：xiaolewennofollow

今天小編就為大家分享一篇python實(shí)現(xiàn)高斯判別分析算法的例子，具有很好的參考價(jià)值，希望對(duì)大家有所幫助。一起跟隨小編過來看看吧

高斯判別分析算法（Gaussian discriminat analysis）

高斯判別算法是一個(gè)典型的生成學(xué)習(xí)算法（關(guān)于生成學(xué)習(xí)算法可以參考我的另外一篇博客）。在這個(gè)算法中，我們假設(shè)p(x|y)p(x|y)服從多元正態(tài)分布。

注：在判別學(xué)習(xí)算法中，我們假設(shè)p(y|x)p(y|x)服從一維正態(tài)分布，這個(gè)很好類比，因?yàn)樵谀Ｐ椭休斎霐?shù)據(jù)XX通常是擁有很多維度的，所以對(duì)于XX的條件概率建模時(shí)要取多維正態(tài)分布。

多元正態(tài)分布

多元正態(tài)分布也叫多元高斯分布，這個(gè)分布的兩個(gè)參數(shù)分別是平均向量μ∈Rnμ∈Rn和一個(gè)協(xié)方差矩陣∑∈Rn×n∑∈Rn×n

關(guān)于協(xié)方差矩陣的定義；假設(shè)XX是由nn個(gè)標(biāo)量隨機(jī)變量組成的列向量，并且μkμk是第kk個(gè)元素的期望值，即μk=E(Xk)μk=E(Xk)，那么協(xié)方差矩陣被定義為

下面是一些二維高斯分布的概率密度圖像：

最右邊的圖像展現(xiàn)的二維高斯分布的均值是零向量（2x1的零向量），協(xié)方差矩陣Σ=IΣ=I(2x2的單位矩陣)，像這樣以零向量為均值以單位陣為協(xié)方差的多維高斯分布稱為標(biāo)準(zhǔn)正態(tài)分布，中間的圖像以零向量為均值，Σ=0.6IΣ=0.6I;最右邊的圖像中Σ=2IΣ=2I，觀察發(fā)現(xiàn)當(dāng)ΣΣ越大時(shí)，高斯分布越“鋪開”，當(dāng)ΣΣ越小時(shí)，高斯分布越“收縮”。

讓我們看一些其他例子對(duì)比發(fā)現(xiàn)規(guī)律

上圖中展示的三個(gè)高斯分布對(duì)應(yīng)的均值均為零向量，協(xié)方差矩陣分別對(duì)應(yīng)與下面三個(gè)

最左邊的圖像是我們熟悉的標(biāo)準(zhǔn)二維正態(tài)分布，然后我們觀察到當(dāng)我們?cè)黾应拨驳姆侵鲗?duì)角元素時(shí)，概率密度圖像沿著45°線(x1=x2x1=x2)“收縮”，從對(duì)應(yīng)的等高線輪廓圖可以跟清楚的看到這一點(diǎn)：

通過對(duì)比右邊和中間的兩幅圖發(fā)現(xiàn)，通過減少主對(duì)角元素可以讓概率密度圖像變得“收縮”，不過是在相反的方向上。

高斯判別分析模型

當(dāng)我們處理輸入特征是連續(xù)隨機(jī)變量xx時(shí)的分類問題時(shí)，我們可以使用高斯判別分析模型（GDA），用多元正態(tài)分布模型來描述p(x|y)p(x|y),模型的具體數(shù)學(xué)表達(dá)式是這樣的：

通過最大化似然函數(shù)ll可以得到上面四個(gè)參數(shù)的估計(jì)值：

我們用圖像直觀的描述一下算法處理的結(jié)果：

python的實(shí)現(xiàn)demo 如下：

第57的高斯概率密度函數(shù)用矩陣運(yùn)算寫有bug沒跑通，又因?yàn)閷?shí)驗(yàn)數(shù)據(jù)只有二維，于是在紙上對(duì)上文中矩陣運(yùn)算公式進(jìn)行了化簡(jiǎn)至最后結(jié)果寫在了函數(shù)里。如有疑問可以拿出筆來演算一下。

#GDA
#author:Xiaolewen
import matplotlib.pyplot as plt
from numpy import *

#Randomly generate two cluster data of Gaussian distributions
mean0=[2,3]
cov=mat([[1,0],[0,2]])
x0=random.multivariate_normal(mean0,cov,500).T #The first class point which labael equal 0
y0=zeros(shape(x0)[1])
#print x0,y0
mean1=[7,8]
cov=mat([[1,0],[0,2]])
x1=random.multivariate_normal(mean1,cov,300).T
y1=ones(shape(x1)[1]) #The second class point which label equals 1
#print x1,y1

x=array([concatenate((x0[0],x1[0])),concatenate((x0[1],x1[1]))])
y=array([concatenate((y0,y1))])
m=shape(x)[1]
#print x,y,m
#Caculate the parameters:\phi,\u0,\u1,\Sigma
phi=(1.0/m)*len(y1)
#print phi
u0=mean(x0,axis=1) 
#print u0
u1=mean(x1,axis=1)
#print u1

xplot0=x0;xplot1=x1 #save the original data to plot 
x0=x0.T;x1=x1.T;x=x.T
#print x0,x1,x
x0_sub_u0=x0-u0
x1_sub_u1=x1-u1
#print x0_sub_u0
#print x1_sub_u1
x_sub_u=concatenate([x0_sub_u0,x1_sub_u1])
#print x_sub_u

x_sub_u=mat(x_sub_u)
#print x_sub_u

sigma=(1.0/m)*(x_sub_u.T*x_sub_u)
#print sigma

#plot the discriminate boundary ,use the u0_u1's midnormal
midPoint=[(u0[0]+u1[0])/2.0,(u0[1]+u1[1])/2.0]
#print midPoint
k=(u1[1]-u0[1])/(u1[0]-u0[0])
#print k
x=range(-2,11)
y=[(-1.0/k)*(i-midPoint[0])+midPoint[1] for i in x]



#plot contour for two gaussian distributions
def gaussian_2d(x, y, x0, y0, sigmaMatrix):
 return exp(-0.5*((x-x0)**2+0.5*(y-y0)**2))
delta = 0.025
xgrid0=arange(-2, 6, delta)
ygrid0=arange(-2, 6, delta)
xgrid1=arange(3,11,delta)
ygrid1=arange(3,11,delta)
X0,Y0=meshgrid(xgrid0, ygrid0) #generate the grid
X1,Y1=meshgrid(xgrid1,ygrid1)
Z0=gaussian_2d(X0,Y0,2,3,cov)
Z1=gaussian_2d(X1,Y1,7,8,cov)

#plot the figure and add comments
plt.figure(1)
plt.clf()
plt.plot(xplot0[0],xplot0[1],'ko')
plt.plot(xplot1[0],xplot1[1],'gs')
plt.plot(u0[0],u0[1],'rx',markersize=20)
plt.plot(u1[0],u1[1],'y*',markersize=20)
plt.plot(x,y)
CS0=plt.contour(X0, Y0, Z0)
plt.clabel(CS0, inline=1, fontsize=10)
CS1=plt.contour(X1,Y1,Z1)
plt.clabel(CS1, inline=1, fontsize=10)
plt.title("Gaussian discriminat analysis")
plt.xlabel('Feature Dimension (0)')
plt.ylabel('Feature Dimension (1)')
plt.show(1)

這是最終的擬合結(jié)果：