# \donttest{
# Generate hierarchical linear data
R = 1000
N = 2000
nobs = 5 #number of observations
nvar = 3 #columns
nz = 2
Z = matrix(runif(N*nz),ncol=nz)
Z = t(t(Z)-apply(Z,2,mean))
Delta = matrix(c(1,-1,2,0,1,0), ncol = nz)
tau0 = 0.1
iota = c(rep(1,nobs))
tcomps=NULL
a = diag(1, nrow=3)
tcomps[[1]] = list(mu=c(-5,0,0),rooti=a)
tcomps[[2]] = list(mu=c(5, -5, 2),rooti=a)
tcomps[[3]] = list(mu=c(5,5,-2),rooti=a)
tpvec = c(.33,.33,.34)
ncomp=length(tcomps)
regdata=NULL
betas=matrix(double(N*nvar),ncol=nvar)
tind=double(N)
for (reg in 1:N) {
tempout=bayesm::rmixture(1,tpvec,tcomps)
if (is.null(Z)){
betas[reg,]= as.vector(tempout$x)
}else{
betas[reg,]=Delta%*%Z[reg,]+as.vector(tempout$x)}
tind[reg]=tempout$z
X=cbind(iota,matrix(runif(nobs*(nvar-1)),ncol=(nvar-1)))
tau=tau0*runif(1,min=0.5,max=1)
y=X%*%betas[reg,]+sqrt(tau)*rnorm(nobs)
regdata[[reg]]=list(y=y,X=X,beta=betas[reg,],tau=tau)
}
Prior1=list(ncomp=ncomp)
keep=1
Mcmc1=list(R=R,keep=keep)
Data1=list(list(regdata=regdata,Z=Z))
returns = s_max(R = R, N = N, Data = Data1, s = 1, iterations = 2)
returns
# }
Run the code above in your browser using DataLab