【the EM algorithm】自己动手，丰衣足食。

1. 首先，用MATLAB生成符合双峰正态分布的随机数。

r=0.5;
mu1=50;
sigma1=20;
mu2=200;
sigma2=20;
x=zeros(10000,1);
for i=1:10000
r1=rand;
x(i,1)=(mu2+sigma2*randn)*heaviside(r1-r)+(mu1+sigma1*randn)*heaviside(r-r1);
end
hist(x)

2. 然后就是废话少说，上代码。

import java.io.BufferedReader;
import java.io.FileReader;
import java.io.IOException;public class MatlabData {double[] d = new double[10005];double[] r = new double[10005];int N=10000;private double normal(double mu, double sigma, double x){double exp = -(x-mu)*(x-mu)/(2*sigma*sigma);return 1.0/(sigma*Math.sqrt(2*Math.PI))*Math.pow(Math.E, exp);}//http://blog.163.com/huai_jing@126/blog/static/17186198320119231094873/private void EM(double pai, double mu1, double mu2, double s1, double s2){int i,cnt=0;for(i=1;i<=N;i++)r[i] = 0;double eps = 1e-6;double opai=Integer.MIN_VALUE, omu1=Integer.MIN_VALUE, omu2=Integer.MIN_VALUE, os1=Integer.MIN_VALUE, os2=Integer.MIN_VALUE;do{prt("====================== " + (++cnt) + "\n");for(i=1;i<=N;i++){r[i] = pai*normal(mu2, s2, d[i])/((1-pai)*normal(mu1, s1, d[i])+pai*normal(mu2, s2, d[i]));}double up=0, down=0;for(i=1;i<=N;i++){up+=((1-r[i])*d[i]);down+=(1-r[i]);}omu1 = mu1;mu1 = up/down;up = 0;for(i=1;i<=N;i++){up+=((1-r[i])*(d[i]-omu1)*(d[i]-omu1));}os1 = s1;s1 = Math.sqrt(up/down);up=0; down=0;for(i=1;i<=N;i++){up+=((r[i])*d[i]);down+=(r[i]);}omu2 = mu2;mu2 = up/down;up = 0;for(i=1;i<=N;i++){up+=((r[i])*(d[i]-omu2)*(d[i]-omu2));}os2 = s2;s2 = Math.sqrt(up/down);opai = pai;pai = down / N;prt("mu1 : "+mu1+"\n");prt("mu2 : "+mu2+"\n");prt("s1  : "+s1+"\n");prt("s2  : "+s2+"\n");prt("pai : "+pai+"\n");}while(Math.abs(mu1-omu1)>=eps||Math.abs(mu2-omu2)>=eps||Math.abs(s1-os1)>=eps||Math.abs(s2-os2)>=eps||Math.abs(pai-opai)>=eps);}public static void main(String[] args) throws IOException{MatlabData m = new MatlabData();String iptPath = "data/data.txt";BufferedReader br = new BufferedReader(new FileReader(iptPath));String temp;int i = 0;while((temp = br.readLine())!=null){m.d[++i] = Double.parseDouble(temp);}br.close();prt("end of read : "+i+"!\n");m.EM(0.1, 125.2, 125.4, 6023.4, 6023.4);prt("end of the EM !\n");}private static void prt(String string){System.out.print(string);}
}

结果如下：

end of read : 10000!
====================== 1
...
====================== 9
mu1 : 50.04014214593736
mu2 : 199.9913405589886
s1  : 20.07997207253577
s2  : 19.995272104999646
pai : 0.5016524895149193
end of the EM !

可以跟第一步的参数对比一下，还是蛮准的嘛！

参考文献：

matlab中如何生成符合双峰正态分布的随机数

EM算法（expectation-maximization algorithm）

【the EM algorithm】自己动手，丰衣足食。相关推荐

【转载】（EM算法）The EM Algorithm
(EM算法)The EM Algorithm EM是我一直想深入学习的算法之一,第一次听说是在NLP课中的HMM那一节,为了解决HMM的参数估计问题,使用了EM算法.在之后的MT中的词对齐中也用到了. ...
从0开始编写dapper核心功能、压榨性能、自己动手丰衣足食
我偶然听说sqlsugar的性能比dapper强.对此我表示怀疑(由于我一直使用的dapper存在偏见吧),于是自己测试了sqlsugar.freesql.dapper发现他们的给我的结果是 sqls ...
第八讲：期望最大化算法(EM algorithm)
在前面的若干讲义中,我们已经讲过了期望最大化算法(EM algorithm),使用场景是对一个高斯混合模型进行拟合(fitting a mixture of Gaussians).在本章里面,我们要给 ...
一键免费下载全网在线视频素材,自己动手丰衣足食
最近经常分享一个下载视频和音乐的下载方法,今天看到这个震惊了. 一键免费下载全网在线视频素材,自己动手丰衣足食
xlnt踩坑记2_自己动手丰衣足食
可以参考鄙人上一篇博客 xlnt踩坑记1 然后就这样我爆肝一下午之后终于搞到了xlnt库,他的dll和lib 我开始认识到了--当我开始搞一些比较偏的Project的时候,真的就要靠自己动手丰衣足食 ...
dealunay triangulation 之自己动手丰衣足食
总算搞到了个简化的watson 算法.自己动手丰衣足食.这个是标准的2维watson算法,顺便介绍一下这个算法的特点.和delaunay三角化的定义.对于平面上任意给定的点集,存在一种唯一的三角化,满 ...
（EM算法）The EM Algorithm
EM是我一直想深入学习的算法之一,第一次听说是在NLP课中的HMM那一节,为了解决HMM的参数估计问题,使用了EM算法.在之后的MT中的词对齐中也用到了.在Mitchell的书中也提到EM可以用于贝叶 ...
EM Algorithm
(From yesterday to now, I have learning EM for 4 times, and each time I will gain some new ideas and ...
小白日志——扫灰、加内存条、装系统自己动手丰衣足食
2015-4-4 周围女生的笔记本常常用着用着就很卡,卡了也忍着,实在忍不下去了,如我这个小白,就决定给电脑扫扫灰.小白没有经验,只有一颗"勇猛"的心,遂昨天在网上搜罗一筐资料后, ...

【the EM algorithm】自己动手，丰衣足食。

【the EM algorithm】自己动手，丰衣足食。相关推荐

最新文章

热门文章