##
if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=10000} else {R=10}
set.seed(66)
nvar=5 ## number of coefficients
nlgt=1000 ## number of cross-sectional units
nobs=10 ## number of observations per unit
nz=2 ## number of regressors in mixing distribution
## set hyper-parameters
## B=ZDelta + U
Z=matrix(c(rep(1,nlgt),runif(nlgt,min=-1,max=1)),nrow=nlgt,ncol=nz)
Delta=matrix(c(-2,-1,0,1,2,-1,1,-.5,.5,0),nrow=nz,ncol=nvar)
iota=matrix(1,nrow=nvar,ncol=1)
Vbeta=diag(nvar)+.5*iota%*%t(iota)
## simulate data
lgtdata=NULL
for (i in 1:nlgt)
{ beta=t(Delta)%*%Z[i,]+as.vector(t(chol(Vbeta))%*%rnorm(nvar))
X=matrix(runif(nobs*nvar),nrow=nobs,ncol=nvar)
prob=exp(X%*%beta)/(1+exp(X%*%beta))
unif=runif(nobs,0,1)
y=ifelse(unif<prob,1,0)
lgtdata[[i]]=list(y=y,X=X,beta=beta)
}
Data=list(Dat=lgtdata,Demo=Z)
out=rhierBinLogit(Data=list(lgtdata=lgtdata,Z=Z),Mcmc=list(R=R))
cat("Deltadraws ",fill=TRUE)
mat=apply(out$Deltadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
mat=rbind(as.vector(Delta),mat); rownames(mat)[1]="delta"; print(mat)
cat("Vbetadraws ",fill=TRUE)
mat=apply(out$Vbetadraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
mat=rbind(as.vector(Vbeta),mat); rownames(mat)[1]="Vbeta"; print(mat)
if(0){
td=as.vector(Delta)
par(mfrow=c(2,2))
matplot(out$Deltadraw[,(1:nvar)],type="l")
abline(h=td[1:nvar],col=(1:nvar))
matplot(out$Deltadraw[,((nvar+1):(2*nvar))],type="l")
abline(h=td[(nvar+1):(2*nvar)],col=(1:nvar))
matplot(out$Vbetadraw[,c(1,7,13,19,25)],type="l")
abline(h=1.5)
matplot(out$Vbetadraw[,-c(1,7,13,19,25)],type="l")
abline(h=.5)
}
Run the code above in your browser using DataLab