ICode9

精准搜索请尝试: 精确搜索
首页 > 编程语言> 文章详细

EM算法估计混合高斯模型(GMM)

2019-06-28 09:52:09  阅读:466  来源: 互联网

标签:EM xi 高斯 GMM sum obsX theta 身高 zj


文章目录

1、EM原理

EM本质是上是极大似然估计(MLE)概率模型有时即含有观测变量,又含有隐变量,如果概率模型的变量都是观测变量,只要show出测量数据,可以直接用极大似然估计法,或者用贝叶斯估计法估计模型参数。但是当模型含有隐变量时,就不能简单的用这些估计方法,EM算法就是含有隐变量的概率模型参数的极大似然估计法。
举个例子啥是隐变量,假设学校有100人,我们可以测出这100人的身高数据,我们都知道人的身高有时依赖于性别,但是不知道男女比例,这时男女比可以作为隐变量。

2、啥是混合高斯模型

先上表达式:
p(xj)=k=1MN(xjμk,Σk)P(Gk),j=1,2,,N p(\boldsymbol{x_j})=\sum_{k=1}^{M}N(\boldsymbol{x_j}|\boldsymbol{\mu_k},\Sigma_k)P(G_k),j=1,2,\dots,N p(xj​)=k=1∑M​N(xj​∣μk​,Σk​)P(Gk​),j=1,2,…,N
上式是混合高斯模型的pdf,M表示高斯模型的数量(很多时候可视作分类簇数量),P(Gk)P(G_k)P(Gk​)表示事件k发生的概率,μk,Σk\boldsymbol{\mu_k},\Sigma_kμk​,Σk​分别表示第k个高斯模型的均值和协方差矩阵。

3、EM算法求解GMM

先上似然函数:
L(θ)=i=1Nzp(xi,z;θ) L(\theta)=\prod_{i=1}^N\sum_zp(x_i,z;\theta) L(θ)=i=1∏N​z∑​p(xi​,z;θ)
取对数后:
l(θ)=logi=1Nzp(xi,z;θ)=i=1Nlogzp(xi,z;θ) l(\theta)=log\prod_{i=1}^N\sum_zp(x_i,z;\theta)=\sum_{i=1}^Nlog\sum_zp(x_i,z;\theta) l(θ)=logi=1∏N​z∑​p(xi​,z;θ)=i=1∑N​logz∑​p(xi​,z;θ)
由于log函数满足凹函数性质,由一系列操作可以得到:
l(θ)i=1NzQi(zj)logp(xi,z;θ)Qi(zj) l(\theta)\ge\sum_{i=1}^N\sum_zQ_i(z_j) log\frac{p(x_i,z;\theta)}{Q_i(z_j) } l(θ)≥i=1∑N​z∑​Qi​(zj​)logQi​(zj​)p(xi​,z;θ)​
p(xi,z;θ)Qi(zj)\frac{p(x_i,z;\theta)}{Q_i(z_j) }Qi​(zj​)p(xi​,z;θ)​为常数时,等号成立。
所以E-step:
Qi(zj)=p(zjxi;θ) Q_i(z_j)=p(z_j|x_i;\theta) Qi​(zj​)=p(zj​∣xi​;θ)
M-step就是在E-step上使上述函数值的期望取得最大时参数θ\thetaθ的取值:
θ:=argmaxl(θ \theta:=argmax l(\theta) θ:=argmaxl(θ)
接下来就是将GMM的pdf代入到EM算法步骤中:
第一步:对第i个样本对第j个高斯分布的贡献率:
Qi(z=j)=P(z=jxi;μ,Σ) Q_i(z=j)=P(z=j|x_i;\mu,\Sigma) Qi​(z=j)=P(z=j∣xi​;μ,Σ)
第二步:根据E-step中的Q估计μ,Σ\mu,\Sigmaμ,Σ:
μl=iNQixiiNQi \mu_l=\frac{\sum_i^N Q_i x_i}{\sum_i^NQ_i } μl​=∑iN​Qi​∑iN​Qi​xi​​
Σj=iNQj(i)(x(i)μj)(x(i)μj)TiNQj(i) \Sigma_j=\frac{\sum_i^N Q_j^{(i)} (x^{(i)}-\mu_j)(x^{(i)}-\mu_j)^T}{\sum_i^N Q_j^{(i)}} Σj​=∑iN​Qj(i)​∑iN​Qj(i)​(x(i)−μj​)(x(i)−μj​)T​
对于j事件的发生概P(Gj)P(G_j)P(Gj​):
P(Gj)=iNQj(i)N P(G_j)=\frac{\sum_i^N Q_j^{(i)} }{N} P(Gj​)=N∑iN​Qj(i)​​

4、实例

假设从学校学生中随机选取100位学生,测量这100位学生身高,既有男学生也有女生,现在已知测量数据。并且知道男性女性的身高均服从高斯分布N(μ1,σ1),N(μ2,σ2)N(\mu_1,\sigma1),N(\mu2,\sigma2)N(μ1​,σ1),N(μ2,σ2),
试估计参数μ1,σ1,μ2,σ2\mu_1,\sigma1,\mu2,\sigma2μ1​,σ1,μ2,σ2以及男女比例pmalepfemalep_{male}和p_{female}pmale​和pfemale​
第一步产生这一百个学生的身高数据:

N=100;                %样本数目 
mu1=170;              %男生身高均值
sigma1=0.5;           %男生身高均方根
mu2=160;              %女生身高均值
sigma2=0.1;           %女生身高均方根
pm=0.6;               %男生所占总数比例
pf=1-pm;              %女生所占总数比例
maleDatas=normrnd(mu1,sigma1,[1,N*pm]);%产生男生身高数据服从高斯随机分布
femaleDates=normrnd(mu2,sigma2,[1,N*pf]);%产生女生身高数据服从高斯随机分布
obsX=[maleDatas,femaleDates];
randIdx=randperm(N);
obsX=obsX(randIdx);    %混合男女身高数据
maxiter=100;           %最大迭代次数

第二步,用EM算法进行参数估计:

%% EM算法解GMM模型
% 第一步,初始化参数
estpm=0.5*ones(1,maxiter);                %待估男生比例
estpf=1-estpm;                            %待估女生比例
estmu1=mean(obsX)*ones(1,maxiter);        %待估男生身高均值
estsigma1=sqrt(std(obsX))*ones(1,maxiter);%待估男生身高均方根
estmu2=estmu1;                            %待估女生身高均值
estsigma2=estsigma1+rand*ones(1,maxiter); %待估男女生身高均方根
Q=zeros(2,N);                             %初始化贡献率
for i=2:maxiter
   % 第二步,E-step,计算贡献率
   for j=1:N
        k(1,j)=estpm(i-1)/(sqrt(2*pi)*estsigma1(i-1))*exp(-(obsX(j)-estmu1(i-1))^2/(2*estsigma1(i-1)^2));
        k(2,j)=estpf(i-1)/(sqrt(2*pi)*estsigma2(i-1))*exp(-(obsX(j)-estmu2(i-1))^2/(2*estsigma2(i-1)^2));
        p(j)=k(1,j)+k(2,j);
        Q(1,j)=k(1,j)/p(j);               %计算每个样本点对男生高斯分布贡献率
        Q(2,j)=k(2,j)/p(j);               %计算每个样本点对女生高斯分布贡献率
   end
   % 第三步,M-step,更新参数
   nk=sum(Q,2);
   estmu1(i)=Q(1,:)*obsX'/nk(1);
   estmu2(i)=Q(2,:)*obsX'/nk(2);
   estsigma1(i)=sqrt(sum(Q(1,:).*(obsX-estmu1(i)).^2)/nk(1));
   estsigma2(i)=sqrt(sum(Q(2,:).*(obsX-estmu2(i)).^2)/nk(2));
   estpm(i)=nk(1)/N;
   estpf(i)=nk(2)/N;
end

结果:
1
在这里插入图片描述
可以看出经过100次迭代后,GMM参数估计得到了很好的估计。

标签:EM,xi,高斯,GMM,sum,obsX,theta,身高,zj
来源: https://blog.csdn.net/baidu_36161424/article/details/93906383

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有