Can anybody help me to correct code the EM algorithm for two-component Laplace mixture models?
library(LaplacesDemon)
library(base)
EM_TwoMixtureLaplace = function(alpha.1, mu1, mu2, lambda1, lambda2, x, maxiter=1000, tol=1e-5)
{
diff=1
iter=0
while (diff>tol & iter<maxiter) {
## E-step: compute omega:
d1=dlaplace(x, mu1, lambda1)
d2=dlaplace(x, mu2, lambda2)
omega1=d1alpha.1/(d1p+d2*(1-alpha.1))
omega2=1-omega1
# compute density in two groups
## M-step: update p, mu and sd
resid1=x-mu1
resid2=x-mu2
p.new=mean(omega1)
mu1.new=sum(omega1median(x))
mu2.new=sum(omega2median(x))
lambda1.new=sum(omega1abs(resid1)) / sum(omega1)
lambda2.new=sum(omega2abs(resid2)) / sum(omega2)
## calculate diff to check convergence
diff=sqrt(sum((mu1.new-mu1)^2+(mu2.new-mu2)^2+(lambda1.new-lambda1)^2+(lambda2.new-lambda2)^2))
alpha.1=p.new;
mu1=mu1.new;
mu2=mu2.new;
lambda1=lambda1.new;
lambda2=lambda2.new;
iter=iter+1;
cat("Iter", iter, ": mu1=", mu1.new, ", mu2=",mu2.new, ", lambda1=",lambda1.new,
", lambda2=",lambda2.new, ", alpha.1=", p.new, ", diff=", diff, "\n")
}
}
simulation
n=1000
alpha1=0.3
x1=rlaplace(nalpha1,-2,1)
x2=rlaplace(n(1-alpha1),10,1)
componentLaplace2.M=x= c(x1,x2)
median(x)
hist(x, 50)
ggplot(data.frame(componentLaplace2.M), aes(x = componentLaplace2.M)) +
geom_histogram(fill = "purple", binwidth = 0.1)
initial values for EM
alpha.1=0.5
mu1=quantile(x, 0.1);
mu2=quantile(x, 0.9)
lambda1=1
lambda2=3
c(alpha.1, mu1, mu2, lambda1, lambda2)
EM_TwoMixtureLaplace(alpha.1, mu1, mu2, lambda1, lambda2, x)