Hi, R users

  I am trying to use EM algorithmto estimate a latent class model of discrete choice. The basic model is a logit model which has two variables X=(X1,X2) and the utility is defined as v = Beta*X. There are 3 classes of individuals each has its own Beta values, so Beta is a 3*2 matrix. For each individual, (there are 1000), he make a choice between two randomly generated choice alternatives and has to choose one. The proportion of classes is set to (0.2,0.5,0.3). based on this seting, the proportion parameters alpha always end up with one of them goes to 0 in EM estimation. I checked the code and can not find what's the problem. Maybe this is a bit off topic in R help, but I code it in R and hence in hope that some one here can figure it out. 

Plus, I have a normal mixture model with success EM result. 

Thanks in advance.

Following is the code:

#=========================

X = array(,dim=c(1000,2))#first alternative
b= c(2.3,0.7,0.3,0.7,0.1,0.8)
alpha = c(0.3,0.5,0.2)#proportion

X[,1] = runif(1000,min=5,max=20)
X[,2] = runif(1000,min=5,max=20)
##bi = 0.3
Y = array(,dim=c(1000,2))#second alternative

Y[,1] = runif(1000,min=5,max=20)
Y[,2] = runif(1000,min=5,max=20)

V11 = X[1:300,]%*%b[1:2] 
V12 = Y[1:300,]%*%b[1:2] 
V21 = X[301:800,]%*%b[3:4]
V22 = Y[301:800,]%*%b[3:4]
V31 = X[801:1000,]%*%b[5:6]
V32 = Y[801:1000,]%*%b[5:6]

V1 = rbind(V11,V21,V31)+rnorm(1000)
V2 = rbind(V12,V22,V32)+rnorm(1000)

oo = exp(V1)+exp(V2)

P1 = exp(V1)/oo
P2 = exp(V2)/oo
D = ifelse(V1>V2,0,1)

#second part of Q function
Q2 = function(Beta,H){
 probs = logProbInd(Beta)
 li = sum(H*probs)
 return(li)
}


logProbInd=function(Beta){#X, Y, D values take as is from environment
 dim(Beta) = c(2,3)
 Beta = t(Beta)
 probs = matrix(,nrow = 1000, ncol = 3)
 for(i in 1:3){
  v1 = X%*%Beta[(i-1)*2+1:2]
  v2 = Y%*%Beta[(i-1)*2+1:2]
  p1 = exp(v1)/(exp(v1)+exp(v2))
  p2 = exp(v2)/(exp(v1)+exp(v2))
  probs[,i] = ifelse (D==0,log(p1),log(p2))
  
 }
 
 return (probs)
}

#H [individuals][class]
E_step = function(alpha,Beta){#calc posterior of H
 
 tmpH = matrix(,nrow = 1000,ncol =3)
 lprobs = logProbInd(Beta)
 for(i in 1:3){#classes
  tmpH[,i] = alpha[i]*exp(lprobs[,i])
 }
 H = tmpH /apply(tmpH,1,sum)
 return( H)
}

M_step = function(H,Beta){
 #first part use direct estimation
 aita = apply(H,2,sum)/1000
 
 opt.c = optim(Beta,Q2,H=H,method="BFGS",control = list(fnscale = -1))
 
 lik = opt.c$value
 return(c(aita,opt.c$par,lik))
}


#EM loops
alf = c(0.33,0.33,0.33)
Bt = seq(0.1,0.6,by=0.1)

sc = rep(-8000,5)
i=1

while(T){
 H= E_step(alpha=alf,Beta=Bt)
 theta = M_step(H=H,Beta=Bt)
 print(theta)
 alf = theta[1:3]
 Bt = theta[4:9]
 
 #check convergence
 sc[(i%%5) +1] = theta[10]
 i=i+1
 #if((theta[10] - mean(sc) ) < 0.0005) break
}


	[[alternative HTML version deleted]]

