Dear all, I need help with EM algorithm. I am modeling this alogirthm based on " Incomplete Data in Generalized Linear Models" by Joseph G Ibrahim. i have half way through developing the R code based on this this paper, I need little help in tweaking my code furthure. ---------------------------------------------------------------------------- a=read.table("C:/sudhi/ex3.csv",header=T,sep=",") #Set design matrix n=1558 X1=as.matrix(rep(1,n)) X2=as.matrix((a[,1])) X3=as.matrix((a[,2])) X4=as.matrix((a[,3])) X5=as.matrix((a[,4])) X6=as.matrix((a[,5])) X7=as.matrix((a[,6])) X8=as.matrix((a[,7])) X=cbind(X1,X2,X3,X4,X5,X6,X7,X8) #Start Newton-Raphson Method #Set Initial parameter B0=0.001 B1=0.001 B2=0.001 B3=0.001 B4=0.001 B5=0.001 B6=0.001 B7=0.001 #Parameters B=matrix(c(B0,B1,B2,B3,B4,B4,B6,B7),nrow=8,ncol=1,byrow=F) #Start Iteration Diff=1 AA = 0 while(Diff>.000001){ AA= AA +1 plot_colors <- c("blue","red","forestgreen") plot(B, type="l", lwd=2, col=plot_colors[1], add=TRUE) split.screen(c(1,1)) #Probability for each comb P=(1+exp(-(X%*%B)))^-1 #Set V matrix ni=1 V=diag(n) for(i in 1:n) { V[i,i]=diag(ni%*%P[i]%*%(1-P[i])) } #Set S matrix #Response R=a[,8] S=matrix(data=1,nrow=n,ncol=1) for (i in 1:n) {S[i]=R[i]-ni*P[i]} new.B=B+solve(t(X)%*%V%*%X)%*%(t(X)%*%S) Diff=max(abs(new.B-B)) B=new.B} #Check the second derivative -diag(t(X)%*%V%*%X) #Check the first derivative t(X)%*%S #Maximum Likelihood Estimates --------------------------------------------------- Please see the pdf file attached for Joseph G Ibrahim's logic Sincerely Sudhi http://r.789695.n4.nabble.com/file/n4167861/EM_on_Logit_Model_-Ibrahim.PDF EM_on_Logit_Model_-Ibrahim.PDF -- View this message in context: http://r.789695.n4.nabble.com/EM-Algorithm-for-missing-data-tp4167861p4167861.html Sent from the R help mailing list archive at Nabble.com.