EM算法R代码及疑问.docx

《EM算法R代码及疑问.docx》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、学习了EM的理论后想写一下练练手。在网上找到了混合高斯的R代码模拟数据是样本集5000个，前2000个是以3为均值，1为方差的高斯分布，后3000个是以-2为均值，2为方差的高斯分布。1.# 模拟数据  2.miu1 <- 3  3.miu2 <- -2  4.sigma1 <- 1  5.sigma2 <- 2  6.alpha1 <- 0.4  7.alpha2 <- 0.6  8.# 生成两种高斯分布的样本  9.n <- 5000  10.x <- rep(0,n)  11.n1 <- floor(n*alpha1)  12.n2 <- n - n1  13.x[1:n1] <- r

2、norm(n1)*sigma1 + miu1  14.x[(n1+1):n] <- rnorm(n2)*sigma2 + miu2  15.hist(x,freq=F)  16.lines(density(x),col='red')  17.# 设置初始值  18.m <- 2  19.miu <- runif(m)  20.sigma <- runif(m)  21.alpha <- c(0.2,0.8)  22.prob <- matrix(rep(0,n*m),ncol=m)  23.# 迭代  24.for (step in 1:100){  25.# E步，已知参数求概率  26.

3、for (j in 1:m){  27.        prob[,j]<- sapply(x,dnorm,miu[j],sigma[j])  28.    }  29.    sumprob <- rowSums(prob)  30.    prob<- prob/sumprob     ####求出样本属于各类的概率  31.  32.  33.    oldmiu <- miu  34.    oldsigma <- sigma  35.    oldalpha <- alpha  36.  37.  1.## M步 更新参数  2. for (j in 1:m){  3.      

4、  p1 <- sum(prob[ ,j])  4.        p2 <- sum(prob[ ,j]*x)  5.        miu[j] <- p2/p1  6.        alpha[j] <- p1/n  7.        p3 <- sum(prob[ ,j]*(x-miu[j])^2)  8.        sigma[j] <- sqrt(p3/p1)  9.    }  10.  11.  12.## 结束迭代的条件  13.epsilo <- 1e-4  14.    if (sum(abs(miu-oldmiu))

5、sum(abs(sigma-oldsigma))

6、z=j)/∑P(xi

7、zj)计算。但是在看理论的时候E步计算的应该是P(z=j

8、xi)，也就是z=j的后验概率。用贝叶斯公式，后验概率可以从先验概率计算得到。（这玩意怎么编辑公式？

9、只能截图了）其中P(x

10、z)是高斯分布的概率密度，P(z)就是代码中的alpha。于是我修改了一下上面代码中E步的prob参数计算的方法：1.for (j in 1:m){  2.   prob[,j]<- sapply(x,dnorm,miu[j],sigma[j])  3.}  4.  5.for(i in (1:dim(prob)[1])){  6.   px<-alpha%*%prob[i,]  7.   prob[i,]<-(prob[i,]*alpha)/px  8.}  这样计算的prob才是z=j的后验概率。分别贴出来修改前后的计算结果：修改前：step35miu2.8439

11、96-2.195416sigma1.0973231.849515alpha0.44222490.5577751 到第35次迭代时收敛，得到的均值是2.843996，-2.195416（真实值是3，-2），方差1.097323，1.849515（真实值是1，2），alpha 0.4422249，0.5577751（真实值0.4,0.6）。修改后的：step57miu-1.9246973.005561sigma2