EM算法具體過程看前一篇博客

一、matlab實現

1.matlab代碼

close all; clear; clc;%% M=3; % 高斯數量 N=600; % 數據樣本總數 th=0.000001; % 聚合閥值 K=2; % 輸出信號保留% 待生成數據的參數 a_real =[2/3;1/6;1/6]; mu_real=[3 4 6;5 3 7]; cov_real(:,:,1)=[5 0;0 0.2]; cov_real(:,:,2)=[0.1 0;0 0.1]; cov_real(:,:,3)=[0.1 0;0 0.1]; % 這里生成的數據全部符合標準 x=[ mvnrnd( mu_real(:,1) , cov_real(:,:,1) , round(N*a_real(1)) )' ,...mvnrnd( mu_real(:,2) , cov_real(:,:,2) , round(N*a_real(2)) )' ,...mvnrnd( mu_real(:,3) , cov_real(:,:,3) , round(N*a_real(3)) )' ];figure(1),plot(x(1,:),x(2,:),'.')%% EM % 參數初始化 a=[1/3,1/3,1/3]; %各類的比例（權重） mu=[1 2 3; %均值初始化2 1 4]; cov(:,:,1)=[1 0; %協方差初始化0 1]; cov(:,:,2)=[1 0;0 1]; cov(:,:,3)=[1 0;0 1];t=inf; count=0; figure(2),hold on while t>=tha_old = a;mu_old = mu;cov_old= cov; rznk_p=zeros(M,N);%生成M行N列零矩陣for cm=1:Mmu_cm=mu(:,cm);cov_cm=cov(:,:,cm);for cn=1:N %計算Pi(x)p_cm=exp(-0.5*(x(:,cn)-mu_cm)'/cov_cm*(x(:,cn)-mu_cm));rznk_p(cm,cn)=p_cm;endrznk_p(cm,:)=rznk_p(cm,:)/sqrt(det(cov_cm));endrznk_p=rznk_p*(2*pi)^(-K/2); %E step%開始求rznk 相當于Pr(i|Xt)rznk=zeros(M,N);%r(Zpikn=zeros(1,M);%r(Zpikn_sum=0;for cn=1:Nfor cm=1:M%計算p(x|*)概率分布pikn(1,cm)=a(cm)*rznk_p(cm,cn); % pikn_sum=pikn_sum+pikn(1,cm);endfor cm=1:Mrznk(cm,cn)=pikn(1,cm)/sum(pikn);endend%求rank結束 % M stepnk=zeros(1,M);for cm=1:Mfor cn=1:Nnk(1,cm)=nk(1,cm)+rznk(cm,cn);endenda=nk/N;rznk_sum_mu=zeros(M,1);% 求均值MU%nk(cm) 相當于ni%rznk_sum_mu ni*Xtfor cm=1:Mrznk_sum_mu=0;for cn=1:Nrznk_sum_mu=rznk_sum_mu+rznk(cm,cn)*x(:,cn);endmu(:,cm)=rznk_sum_mu/nk(cm);end% 求協方差COV for cm=1:Mrznk_sum_cov=zeros(K,K);for cn=1:N%求協方差（Hi-u)^2rznk_sum_cov=rznk_sum_cov+rznk(cm,cn)*(x(:,cn)-mu(:,cm))*(x(:,cn)-mu(:,cm))';endcov(:,:,cm)=rznk_sum_cov/nk(cm);endt=max([norm(a_old(:)-a(:))/norm(a_old(:));norm(mu_old(:)-mu(:))/norm(mu_old(:));norm(cov_old(:)-cov(:))/norm(cov_old(:))]);temp_f=sum(log(sum(pikn)));plot(count,temp_f,'r+')count=count+1; end %while hold off f=sum(log(sum(pikn)));% 輸出結果 a mu covfigure(3), hold on plot(x(1,:),x(2,:),'k.'); plot(mu_real(1,:),mu_real(2,:),'*c'); plot(mu(1,:),mu(2,:),'+r'); hold offfigure(4), hold on for i=1:N[max_temp,ind_temp]=max(rznk(:,i));if ind_temp==1plot(x(1,i),x(2,i),'k.');endif ind_temp==2plot(x(1,i),x(2,i),'b.');endif ind_temp==3plot(x(1,i),x(2,i),'r.');end end%fcm聚類 [center, U, OBJ_FCN]=fcm(x',3); figure(5), hold on for i=1:N[max_temp,ind_temp]=max(U(:,i));if ind_temp==1plot(x(1,i),x(2,i),'k.');endif ind_temp==2plot(x(1,i),x(2,i),'b.');endif ind_temp==3plot(x(1,i),x(2,i),'r.');end endplot(center(:,1),center(:,2),'c*')hold off

2.結果顯示

二、Java實現

編譯工具為：eclipse
1.Java代碼

package cn.sxt.oo2;/*** 一維情況下的EM算法實現*1、求期望（e-step)*2、期望最大化（估值）（M-step)*3、循環以上兩部直到收斂*/ public class MyEM {private static final double[] points={1.0,1.3,2.2,2.6,2.8,5.0,7.3,7.4,7.5,7.7,7.9};private static double[][] w;//權值private static double[] means = {7.7,2.3};//均值private static double[] variances= {1,1};//方差private static double[] probs = {0.5,0.5};//每個類的概率；這里默認選擇k=2了；/*** 高斯分布計算公式，也就是先驗概率*///p(x|c)的計算private static double gaussianPro(double point,double mean,double variance){double prob = 0.0;prob = (1/(Math.sqrt(2*Math.PI)*Math.sqrt(variance)))*Math.exp(-(point-mean)*(point-mean)/(2*variance));return prob;}/*** E-step的主要邏輯*/private static double[][] countPostprob(double[] means,double[] variances,double[] points,double[] probs){int clusterNum = means.length;int pointNum = points.length;double[][] postProbs = new double[clusterNum][pointNum]; double[] denominator = new double[pointNum];for(int m = 0;m <pointNum;m++){denominator[m] = 0.0;for(int n = 0;n<clusterNum;n++){denominator[m]+=(gaussianPro(points[m], means[n], variances[n])*probs[n]);}}for(int i = 0;i<clusterNum;i++){for(int j = 0;j<pointNum;j++){postProbs[i][j]=(gaussianPro(points[j], means[i], variances[i])*probs[i])/(denominator[j]);}}return postProbs;}private static void eStep(){w = countPostprob(means, variances, points, probs);}/*** M-step的主要邏輯之一：由E-step得到的期望，重新計算均值*/private static double[] guessMean(double[][] w,double[] points){int wLength = w.length;double[] means = new double[w.length];double[] wi = new double[wLength];for (int m = 0; m < wLength; m++) {wi[m] = 0.0;for(int n = 0; n<points.length;n++){wi[m] += w[m][n];}}for(int i = 0;i<w.length;i++){means[i] = 0.0;for(int j = 0;j<points.length;j++){means[i]+=(w[i][j]*points[j]);}means[i] /= wi[i];}return means;}/*** M-step的主要邏輯之一：由E-step得到的期望，重新計算方差*/private static double[] guessVariance(double[][] w,double[] points){int wLength = w.length;double[] means = new double[w.length];double[] variances = new double[wLength];double[] wi = new double[wLength];for (int m = 0; m < wLength; m++) {wi[m] = 0.0;for(int n = 0; n<points.length;n++){wi[m] += w[m][n];}}means = guessMean(w, points);for(int i = 0;i<wLength;i++){variances[i] = 0.0;for(int j = 0;j<points.length;j++){variances[i] +=(w[i][j]*(points[j]-means[i])*(points[j]-means[i])); }variances[i] /= wi[i];}return variances;}/*** M-step的主要邏輯之一：由E-step得到的期望，重新計算概率**/private static double[] guessProb(double[][] w){int wLength = w.length;double[] probs = new double[wLength];for(int i = 0;i<wLength;i++){probs[i] = 0.0;for(int j = 0;j<w[i].length;j++){probs[i]+=w[i][j];}probs[i] /=w[i].length;}return probs;}private static void mStep(){means = guessMean(w, points);variances = guessVariance(w, points);probs = guessProb(w);}/*** 計算前后兩次迭代的參數的差值* */private static double threshold(double[] bef_values,double[] values){double diff = 0.0;for(int i = 0 ; i < values.length;i++){diff += (values[i]-bef_values[i]);}return Math.abs(diff);}public static void main(String[] args)throws Exception{int k = 2;w = new double[k][points.length];double[] bef_means;double[] bef_var;do{bef_means = means;bef_var = variances;eStep();mStep();}while(threshold(bef_means, means)<0.01&&threshold(bef_var, variances)<0.01);for(double prob:probs)System.out.println(prob); } }

2.結果顯示

總結

以上是生活随笔為你收集整理的EM算法matlab和Java实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。