斯坦福机器学习实现与分析之五(高斯判别分析)

高斯判别分析(GDA)简介


  首先,高斯判别分析的作用也是用于分类。对于两类样本,其服从伯努利分布,而对每个类中的样本,假定都服从高斯分布,则有:

yBernouli(ϕ)

x|y=0N(μ0,Σ)

x|y=1N(μ1,Σ)

  这样,根据训练样本,估计出先验概率以及高斯分布的均值和协方差矩阵(注意这里两类内部高斯分布的协方差矩阵相同),即可通过如下贝叶斯公式求出一个新样本分别属于两类的概率,进而可实现对该样本的分类。

 p(y|x)=p(x|y)p(y)p(x)

y=argmaxyp(y|x)=argmaxyp(x|y)p(y)p(x)=argmaxyp(x|y)p(y) 

 

GDA详细推导


  那么高斯判别分析的核心工作就是估计上述未知量ϕ,μ0,μ1,Σ。如何来估计这些参数?又该最大似然估计上场了。其对数似然函数为:

l(ϕ,μ0,μ1,Σ)=logi=1mp(x(i),y(i))=logi=1mp(x(i)|y(i))p(y(i))=i=1mlogp(x(i)|y(i))+i=1mlogp(y(i))=i=1mlog(p(x(i)|y(i)=0)1y(i)p(x(i)|y(i)=1)y(i))+i=1mlogp(y(i))=i=1m(1y(i))logp(x(i)|y(i)=0)+i=1my(i)logp(x(i)|y(i)=1)+i=1mlogp(y(i))

  注意此函数第一部分只和μ0,Σ有关,第二部分只和μ1,Σ有关,最后一部分只和ϕ有关。最大化该函数,首先求ϕ,先对其求偏导数:

l(ϕ,μ0,μ1,Σ)ϕ=i=1mlogp(y(i))ϕ=i=1mlogϕy(i)(1ϕ)1y(i))ϕ=i=1my(i)logϕ+(1y(i))log(1ϕ)ϕ=i=1m(y(i)1ϕ(1y(i))11ϕ)=i=1m(I(y(i)=1)1ϕI(y(i)=0)11ϕ)

  此处I为指示函数。令其为0,可求解出:

ϕ=I(y(i)=1)I(y(i)=0)+I(y(i)=1)=I(y(i)=1)m 

  同样地,对μ0求偏导数:

l(ϕ,μ0,μ1,Σ)μ0=i=1m(1y(i))logp(x(i)|y(i)=0)μ0=i=1m(1y(i))(log1(2π)n|Σ|12(x(i)μ0)TΣ1(x(i)μ0))μ0=i=1m(1y(i))(Σ1(x(i)μ0))=i=1mI(y(i)=0)Σ1(x(i)μ0) 

  令其为0,可求解得:

μ0=i=1mI(y(i)=0)x(i)i=1mI(y(i)=0) 

  根据对称性可直接得出:

μ1=i=1mI(y(i)=1)x(i)i=1mI(y(i)=1) 

  下面对Σ求偏导数,由于似然函数只有前面两部分与Σ有关,则将前两部分改写如下:

i=1m(1y(i))logp(x(i)|y(i)=0)+i=1my(i)logp(x(i)|y(i)=1)=i=1m(1y(i))(log1(2π)n|Σ|12(x(i)μ0)TΣ1(x(i)μ0))+i=1my(i)(log1(2π)n|Σ|12(x(i)μ1)TΣ1(x(i)μ1))=i=1mlog1(2π)n|Σ|12i=1m(x(i)μy(i))TΣ1(x(i)μy(i))=i=1m(n2log(2π)12log(|Σ|))12i=1m(x(i)μy(i))TΣ1(x(i)μy(i))

  进而有:

l(ϕ,μ0,μ1,Σ))Σ=12i=1m(1|Σ||Σ|Σ1)12i=1m(x(i)μy(i))(x(i)μy(i))TΣ1Σ=m2Σ112i=1m(x(i)μy(i))(x(i)μy(i))T(Σ2))

  这里推导用到了:

|Σ|Σ=|Σ|Σ1

Σ1Σ=Σ2 

  令其为0,从而求得:

Σ=1mi=1m(x(i)μy(i))(x(i)μy(i))T 

  上面的推导似乎很复杂,但其结果却是非常简洁。通过上述公式,所有的参数都已经估计出来,需要判断一个新样本x时,可分别使用贝叶斯求出p(y=0|x)和p(y=1|x),取概率更大的那个类。

  实际计算时,我们只需要比大小,那么贝叶斯公式中分母项可以不计算,由于2个高斯函数协方差矩阵相同,则高斯分布前面那相同部分也可以忽略。实际上,GDA算法也是一个线性分类器,根据上面推导可以知道,GDA的分界线(面)的方程为:

        (1ϕ)exp((xμ0)TΣ1(xμ0)=ϕexp((xμ1)TΣ1(xμ1)

  取对数展开后化解,可得:

 2xTΣ1(μ1μ0)=μ1TΣ1μ1μ0TΣ1μ0+logϕlog(1ϕ)

  若A=2Σ1(μ1μ0)=(a1,a2,...,an)b=μ1TΣ1μ1μ0TΣ1μ0+logϕlog(1ϕ),则

 a1x1+a2x2+...+anxn=b

  这就是GDA算法的线性分界面。

 

GDA实现


  这里也采用前面讲逻辑回归生成的数据来进行实验,直接load进来进行处理,详见逻辑回归。GDA训练代码如下: 

复制代码
 1 %mu=[mu0 mu1]
 2 function [mu sigma phi]=GDA_train(Sample)
 3     [m, n] = size(Sample); %m个样本,每个n维
 4     Y = Sample(:, end);
 5     X = [Sample(:,1:end-1)];
 6     
 7     idx = find(Y==0);
 8     mu(:,1)=mean(X(idx,:));
 9     
10     idx2 = find(Y==1);    
11     mu(:,2)=mean(X(idx2,:));
12     
13     phi = size(idx2,1)/m;
14     
15     sigma = zeros(n-1);
16     for i = 1:m
17         x = X(i, :)';
18         muc = mu(:, Y(i) + 1);
19         sigma = sigma + (x - muc) * (x - muc)';
20     end
21     sigma = sigma / m;    
22 end
View Code
复制代码

   测试代码: 

复制代码
 1 function GDA_test
 2     clear
 3     close all
 4     clc
 5     load('log_data.mat');
 6     [mu sigma phi]=GDA_train(Sample);
 7     
 8     %显示结果,以下代码不通用,样本维数增加时显示不可用
 9     figure,
10     idx = find(Sample(:,3)==1);
11     plot(Sample(idx,1), Sample(idx,2), 'g*');hold on
12     idx = find(Sample(:,3)==0);
13     plot(Sample(idx,1), Sample(idx,2), 'ro');hold on
14     
15     [t1 t2]=meshgrid(min(Sample(:,1)):.1:max(Sample(:,1)), min(Sample(:,2)):.1:max(Sample(:,2)));
16     G1=Gaussian(mu(:,1),sigma,t1,t2,1-phi);
17     G2=Gaussian(mu(:,2),sigma,t1,t2,phi);
18     contour(t1,t2,G1);hold on
19     contour(t1,t2,G2);hold on
20     
21     A=2*inv(sigma)*(mu(:,2)-mu(:,1));
22     b=mu(:,2)'*inv(sigma)*mu(:,2) - mu(:,1)'*inv(sigma)*mu(:,1) + log(phi) - log(1-phi);
23     x1=min(Sample(:,1)):.1:max(Sample(:,1));
24     x2=(b-A(1)*x1)/A(2);
25     plot(x1,x2,'m')
26 end
27 
28 function G = Gaussian(mu,sigma,t1,t2,phi)
29     for i=1:size(t1,1)
30         for j=1:size(t1,2)
31             x=[t1(i,j);t2(i,j)];
32             z=(x-mu)'*inv(sigma)*(x-mu);
33             G(i,j)=phi*exp(-z)/0.001;
34         end
35     end    
36 end
View Code
复制代码

  训练结果如下,训练样本中,正负样本均为100个,故ϕ=0.5

 

  改变正负样本数量,即相当于改变先验概率,则实验结果如下(相应的ϕ的值显示在图像标题):


       

 

算法分析


  1.与逻辑回归的关系

    根据上面的结果以及贝叶斯公式,可有

 p(y=1|x)=p(x|y=1)p(y=1)p(x)=N(μ1,Σ)ϕN(μ0,Σ)(1ϕ)+N(μ1,Σ)ϕ=1/(1+N(μ0,Σ))N(μ1,Σ)1ϕϕ)

    而

N(μ0,Σ)N(μ1,Σ)=exp{(xμ0)TΣ1(xμ0)(xμ1)TΣ1(xμ1)}=exp{2(μ1μ0)TΣ1x+(μ0TΣμ0μ1TΣμ1)}

    那么,令

2Σ1(μ1μ0)=(θ1,θ2,...,θn)Tθ0=μ0TΣμ0μ1TΣμ1+log1ϕϕ

    则

p(y=1|x)=11+exp(θ0+θ1x1+θ2x2+...+θnxn)

    这不就是逻辑回归的形式么?

    在推导逻辑回归的时候,我们并没有假设类内样本是服从高斯分布的,因而GDA只是逻辑回归的一个特例,其建立在更强的假设条。故两者效果比较:

    a.逻辑回归是基于弱假设推导的,则其效果更稳定,适用范围更广

    b.数据服从高斯分布时,GDA效果更好

    c.当训练样本数很大时,根据中心极限定理,数据将无限逼近于高斯分布,则此时GDA的表现效果会非常好

 

  2.为何要假设两类内部高斯分布的协方差矩阵相同?

    从直观上讲,假设两个类的高斯分布协方差矩阵不同,会更加合理(在混合高斯模型中就是如此假设的),而且可推导出类似上面简洁的结果。

    假定两个类有相同协方差矩阵,分析具有以下几点影响:

    A.当样本不充分时,使用不同协方差矩阵会导致算法稳定性不够;过少的样本甚至导致协方差矩阵不可逆,那么GDA算法就没法进行

    B.使用不同协方差矩阵,最终GDA的分界面不是线性的,同样也推导不出GDA的逻辑回归形式

 

  3.使用GDA时对训练样本有何要求?

    首先,正负样本数的比例需要符合其先验概率。若是预先明确知道两类的先验概率,那么可使用此概率来代替GDA计算的先验概率;若是完全不知道,则可以公平地认为先验概率为  50%。

    其次,样本数必须不小于样本特征维数,否则会导致协方差矩阵不可逆,按照前面分析应该是多多益善。

         

 

posted @   过客冲冲  阅读(11803)  评论(16编辑  收藏  举报
编辑推荐:
· Linux系列:如何用 C#调用 C方法造成内存泄露
· AI与.NET技术实操系列(二):开始使用ML.NET
· 记一次.NET内存居高不下排查解决与启示
· 探究高空视频全景AR技术的实现原理
· 理解Rust引用及其生命周期标识(上)
阅读排行:
· 阿里最新开源QwQ-32B,效果媲美deepseek-r1满血版,部署成本又又又降低了!
· 单线程的Redis速度为什么快?
· SQL Server 2025 AI相关能力初探
· 展开说说关于C#中ORM框架的用法!
· AI编程工具终极对决:字节Trae VS Cursor,谁才是开发者新宠?
点击右上角即可分享
微信分享提示