# components to which each variance covariance matrix is estimated
print(X)
# variance covariance matrices for the components
print(Ds)
# model matrix for each component
Zs=lapply(data.frame(X),function(x){model.matrix(~x-1)})
print(Zs)
# heteroscedastic errors (sigma's)
sigmas=c(1,2)
Sigma=diag(c(Zs[[1]]%*%sigmas))
print(Sigma)
# expected variance covariance matirx
V=array(0,c(8,8))
for(i in 1:3){
	V=V+Zs[[i]]%*%Ds[[i]]%*%t(Zs[[i]])
}
V=V+Sigma
print(V)
# Scaling parameter for gamma distribution
Nu=10
# observed data matrix
Y0=matrix(rnorm(8000),1000)%*%chol(V)/sqrt(rgamma(1000,Nu/2,Nu/2))

# model fitting
# Covariance matrix type:
#  1: full matrix
#  2: diagonal different variances
#  3: diagonal same variance
res=lmm(Y0,Zs[[1]],Zs,covType=c(1,2,3))

# estimated sigma's
print(res[[1]])
# estimated variance covariance matrices for the components
print(res[[2]])
# estimated Nu
print(res[[7]])