Hi, i would to know, if someone have ever write the code to estimate the
parameter (mixing proportion, mean, a var/cov matrix) of a mixture of two
multivariate normal distribution. I wrote it and it works (it could find
mean and mixing proportion, if I fix the var/cov matrix), while if I fix
anything, it doesn't work. My suspect is that when the algorithm iterates
the var/cov matrix, something wrong happens, but i don't know how solve this
problem.
I enclose my script, if someone would give it a look, I would appreciate so
much.
thank you
daniele
#########################
#Start Script
#########################
library(mvtnorm)
libray(scatterplot3d)
library(MASS)
n=100
m1=c(5,-5)
m2=c(-3,3)
s1=matrix(c(2,1,1,3), 2,2)
s2=matrix(c(4,1,1,6), 2,2)
alpha=0.3
c1 <- mvrnorm(round(n*alpha),m1,s1)
c2 <- mvrnorm(round(n*(1-alpha)),m2,s2)
allval <- rbind(c1,c2)
x <- allval[sample(n,n),]
mixm<-
function(x,m1,m2,s1,s2, alpha)
{
f1<-dmvnorm(x, m1, s1, log=FALSE)
f2<-dmvnorm(x, m2, s2, log=FALSE)
f=alpha*f1+(1-alpha)*f2
f
}
plot(x)
den<-mixm(x,m1,m2,s1,s2,alpha)
scatterplot3d(x[,1],x[,2],den, highlight.3d=TRUE, col.axis="blue",
col.grid="lightblue",angle=120, pch=20)
em2mn<- function(y)
{
n<-length(y[,1])
p<-matrix(0,n,1)
f1<-matrix(0,n,1)
f2<-matrix(0,n,1)
tau<-matrix(0,n,2)
eps<-0.0001
mu01<-c(0,0)
mu02<-c(0,0)
sd01<-matrix(0,2,2)
sd02<-matrix(0,2,2)
cov1<-matrix(0,2,2)
cov2<-matrix(0,2,2)
# 1 inizializzare i valori
alpha0= runif(1,0,1)
for (j in 1:2) {
mu01[j] <- runif(1,min=quantile(y[,j], probs =0.25),
max=quantile(y[,j], probs =0.75))
mu02[j] <- runif(1,min=quantile(y[,j], probs =0.25),
max=quantile(y[,j], probs =0.75))
}
sd01<- var(y[1:round(n/2),])
sd02<- var(y[round(n/2):n,])
#sd01<-s1
#sd02<-s2
#prima iterazione
iter<-0
f1<-dmvnorm(y, mu01, sd01, log=FALSE)
f2<-dmvnorm(y, mu02, sd02, log=FALSE)
p=alpha0*f1+(1-alpha0)*f2
scatterplot3d(y[,1], y[,2], p , highlight.3d=TRUE, col.axis="blue",
col.grid="lightblue",main="val iniz mistura normali
multivariata",
angle=120, pch=20)
#verosimiglianza iniziale
l0<-sum(log(p))
l1<-l0
alpha<-alpha0
mu1<-mu01
mu2<-mu02
sd1<-sd01
sd2<-sd02
for (iter in 1:itermax)
{
#passo E
for (i in 1:n) {
tau[i,1]<-(alpha*f1[i])/p[i]
tau[i,2]<-((1-alpha)*f2[i])/p[i]
}
#passo M
alpha= mean(tau[,1])
mu1=colSums(tau[,1]*y)/sum(tau[,1])
mu2=colSums(tau[,2]*y)/sum(tau[,2])
ycen1<-(y-mu1)
ycen2<-(y-mu2)
cov1<-matrix(0,2,2)
cov2<-matrix(0,2,2)
for (i in 1:n){
cov1<-cov1+ (tau[i,1]*(ycen1[i,])%*%t(ycen1[i,]))
cov2<-cov2+ (tau[i,2]*(ycen2[i,])%*%t(ycen2[i,]))
}
# w1<-sqrt(tau[,1])
# w2<-sqrt(tau[,2])
# ywei1<-w1*ycen1
# ywei2<-w2*ycen2
sd1<-cov1/sum(tau[,1])
sd2<-cov2/sum(tau[,2])
f1<-dmvnorm(y,mu1,sd1,log=FALSE)
f2<-dmvnorm(y,mu2,sd2,log=FALSE)
p<-alpha*f1+(1-alpha)*f2
L<-sum(log(p))
cat(iter,L,"\t",alpha,"\t",mu1,"\t",mu2,"\n")
if (abs(L1-L)<eps) break
L1<-L
}
denfin<-mixm(x,mu1,mu2,sd1,sd2,alpha)
scatterplot3d(x[,1],x[,2],denfin, highlight.3d=TRUE, col.axis="blue",
col.grid="lightblue",main="densità valori
finali",angle=120, pch=20)
return(list(alpha0=alpha0,alpha=alpha,mu1=mu1,mu2=mu2,sd1=sd1,sd2=sd2))
}
em2mn(x)
##############################
end script
##############################
--
Dr. Daniele Riggi, PhD student
University of Milano-Bicocca
Department of Statistics
Building U7, Via Bicocca degli Arcimboldi, 8
20126 Milano, Italy
cell. +39 328 3380690
mailto: daniele.riggi@gmail.com
[[alternative HTML version deleted]]
