eb.ran.beta.mix_function(ymat, xmat, ind.r, beta, bv, ev, mispat)
{
#input##############################################################	
# ymat: outcome matrix & ymat[i,] for ith subject                  #
# xmat: common time covariates s.t. xmat*beta=y                    #
# beta: the list of coefficients estimates for each class          #
# ind.r: the list of column indices corresponding to random effects#
# bv  : the list of covariance matrix for latent variables         #
# ev  : the list of within subject covariance matrix               #
#mispat: missing pattern ( 1 for missing)                          # 
####################################################################

## output ##########################################################
#eb.beta: EB estimates of random growth factors w.r.t. classes     #
####################################################################

#*************************** Main *********************************#
   # class number specified in GGMM
	n.class <- length(bv)
	
	# sample size 
	n.obs <- nrow(ymat)
	
	# number of repeated measures
	repn <- ncol(ymat)
	
	# create a list of indices for fixed effect 
	len.beta <- length(beta[[1]])
	betaind <- c(1:length(beta[[1]]))
	ind.f <- ind.r
	for(j in 1:n.class) {
		ind.f[[j]] <- setdiff(betaind, ind.r[[j]])
	}
	
	# create lists of objects for placing class-specific Empirical Bayes (eb) 
	# estimates for growth factors & the associated covariance matrices 
	eb.beta <- as.list(c(1:n.class))
	for(j in 1:n.class) {
		eb.beta[[j]] <- matrix(0, n.obs, len.beta)/0
		bv[[j]] <- as.matrix(bv[[j]])
	}
	
	# compute class-specific eb estimates at level 1 
	for(i in 1:n.obs) {
		# indices for nonmissing components in ymat[i,] 
		s <- sort.list(mispat[i,  ])
		s <- s[1:(repn - sum(mispat[i,  ]))]
		l.s <- length(s)
		
		# note that the following calculation is limited to subjects 
		# with sufficient number of nonmissing 
		if(l.s > (len.beta - 1)) {
			for(j in 1:n.class) {
				u1 <- ginverse(t(xmat[s, ind.r[[j]]]) %*% xmat[s, ind.r[[
					j]]]) %*% t(xmat[s, ind.r[[j]]])
				newv <- ginverse(u1 %*% ev[[j]][s, s] %*% t(u1)) + 
					ginverse(bv[[j]][ind.r[[j]], ind.r[[j]]])
				inv.newv <- ginverse(newv)
				
				# six cases are considered to ensure compatability 
				# among the following matrix operations 
				
					if(length(ind.r[[j]]) > 1 && length(ind.f[[j]]) > 1) {
					eb.beta[[j]][i, ind.r[[j]]] <- t(u1 %*% (ymat[
						i, s] - xmat[s, ind.f[[j]]] %*% beta[[
						j]][ind.f[[j]]])) %*% ginverse(u1 %*% ev[[
						j]][s, s] %*% t(u1)) %*% inv.newv + t(
						beta[[j]][ind.r[[j]]]) %*% ginverse(bv[[
						j]][ind.r[[j]], ind.r[[j]]]) %*% inv.newv
				}
				
	       		if(length(ind.r[[j]]) > 1 && length(ind.f[[j]])==1) {
					eb.beta[[j]][i, ind.r[[j]]] <- t(u1 %*% c(ymat[
						i, s] - xmat[s, ind.f[[j]]] * beta[[j]][
						ind.f[[j]]])) %*% ginverse(u1 %*% ev[[
						j]][s, s] %*% t(u1)) %*% inv.newv + t(
						beta[[j]][ind.r[[j]]]) %*% ginverse(bv[[
						j]][ind.r[[j]], ind.r[[j]]]) %*% inv.newv
				}
				
	    	       if(length(ind.r[[j]]) > 1 && length(ind.f[[j]]) < 1) {
					eb.beta[[j]][i, ind.r[[j]]] <- t(u1 %*% (ymat[
						i, s])) %*% ginverse(u1 %*% ev[[j]][s,
						s] %*% t(u1)) %*% inv.newv + t(beta[[j]][
						ind.r[[j]]]) %*% ginverse(bv[[j]][ind.r[[
						j]], ind.r[[j]]]) %*% inv.newv
				}
					
					if(length(ind.r[[j]]) ==1 && length(ind.f[[j]]) > 1) {
					eb.beta[[j]][i, ind.r[[j]]] <- t(u1 %*% c(ymat[
						i, s] - xmat[s, ind.f[[j]]] %*% beta[[
						j]][ind.f[[j]]])) %*% ginverse(u1 %*% ev[[
						j]][s, s] %*% t(u1)) %*% inv.newv + (
						beta[[j]][ind.r[[j]]]) * ginverse(bv[[
						j]][ind.r[[j]], ind.r[[j]]]) * inv.newv
				}
				   if(length(ind.r[[j]]) ==1 && length(ind.f[[j]]) ==1) {
					eb.beta[[j]][i, ind.r[[j]]] <- t(u1 %*% (ymat[
						i, s] - xmat[s, ind.f[[j]]] * beta[[j]][
						ind.f[[j]]])) %*% ginverse(u1 %*% ev[[
						j]][s, s] %*% t(u1)) %*% inv.newv + (
						beta[[j]][ind.r[[j]]]) * ginverse(bv[[
						j]][ind.r[[j]], ind.r[[j]]]) * inv.newv
				}
				
     				if(length(ind.r[[j]]) ==1 && length(ind.f[[j]]) <1) {
	    			eb.beta[[j]][i, ind.r[[j]]] <- t(u1 %*% (ymat[
						i, s] )) %*% ginverse(u1 %*% ev[[
						j]][s, s] %*% t(u1)) %*% inv.newv + (
						beta[[j]][ind.r[[j]]]) * ginverse(bv[[
						j]][ind.r[[j]], ind.r[[j]]]) * inv.newv
				}

			}
		}
	}
	eb.beta
#*********************************** End ****************************************#	

}