# Simulated example proposed by Nott and Kohn (2005).

n=50
p=15
X = matrix(,n,p)                               # data matrix 
varcov = diag(rep(1,10))                       # variance-covariance matrix
X[,1:10] = rmnorm(n, mean = rep(0,10), varcov) # X1 ... X10 covariates
E = rmnorm(n, mean = rep(0,5), varcov[1:5,1:5])
X[,11:15] = X[,1:5]%*%c(0.3,0.5,0.7,0.9,1.1)%*%t(rep(1,5)) + E # X11.. X15 
e = rnorm(n,0,2.5)                                             # error variance
Y = 2*X[,1]- X[,5] + 1.5*X[,7] + X[,11] + 0.5*X[,13] + e       # response
X = scale(X,TRUE,TRUE)  # standardize the data to have zero mean and unit variance
sign=(1:p)[cor(Y,X)<0]  # make the Pearson correlation to be positive
X[,sign] = -X[,sign]


#=============================================================================
#============== Simulation Study =============================================
# save the results
rmse.y=matrix(,100,7)
rmse.b=matrix(,100,7)
gamma1 = matrix(,100,p) 
gamma2 = matrix(,100,p) 
lambda.post = matrix(,100,4)
gamma3 = matrix(,100,p)
gamma4 = matrix(,100,p)
gamma5 = matrix(,100,p)

# loop for 100 data sets
for(jj in 1:100){
  X = data[[jj]]$X   # data is a list with 100 data sets
  Y = data[[jj]]$Y
  
### BLVS
  dat= blvs(X, Y, lambda=0.446, alpha2=0.0001, gamma2=10000, nburn=1000, ndraw=10000) 
  beta.post = apply((dat[,1:p]*dat[,(p+1):(2*p)]),2,median) # posterior est.
  Yhat = X%*%beta.post
  rmse.y[jj,1] = mean((Y-Yhat)^2)                      #RMSE of fitted values
  rmse.b[jj,1] = mean((beta.post-btrue)^2)             #RMSE of beta est.  
  gamma1[jj,] =  rbind(apply(dat[,(p+1):(2*p)],2,mean))#posterior inclus. probs
  
### BLVS-G  
  dat =  blvs(X, Y, r=6.144, delta=1/13.768, alpha2=0.0001, gamma2=10000, nburn=1000, ndraw=10000) 
  beta.post = apply((dat[,1:p]*dat[,(p+1):(2*p)]),2,median)# posterior est.
  Yhat = X%*%beta.post
  rmse.y[jj,2] = mean((Y-Yhat)^2)                      #RMSE of fitted values
  rmse.b[jj,2] = mean((beta.post-btrue)^2)             #RMSE of beta est.  
  gamma2[jj,] =  apply(dat[,(p+1):(2*p)],2,mean)       #posterior inclus. probs
  lambda.post[jj,1] = mean(dat[,(2*p+2)])              #posterior mean of lambda

### NG
  # monomvn R package is required
  y=Y; x=X
  bols=lm(y~x+0)$coef
  M = sum(bols^2)/p
  reg.blas <- blasso(x, y, T=11000, RJ = FALSE, case="ng", rd=c(2,M/2), icept=FALSE, verb=0, normalize=FALSE)
  beta.post = apply(reg.blas$beta[1001:11000,],2,mean) # posterior est.
  Yhat = X%*%beta.post
  rmse.y[jj,3] = mean((Y-Yhat)^2)                       #RMSE of fitted values
  rmse.b[jj,3] = mean((beta.post-btrue)^2)              #RMSE of beta est.  
  
### NG  with RJ
  # monomvn R package is required
  reg.blas <- blasso(x, y, T=11000, RJ = TRUE, case="ng", rd=c(2,M/2), icept=FALSE, verb=0, normalize=FALSE)
  beta.post = apply(reg.blas$beta[1001:11000,],2,median)# posterior est.
  Yhat = X%*%beta.post
  rmse.y[jj,4] = mean((Y-Yhat)^2)                        #RMSE of fitted values
  rmse.b[jj,4] = mean((beta.post-btrue)^2)               #RMSE of beta est. 
  s = summary(reg.blas, burnin=1000)
  gamma3[jj,] =  as.vector(s$bn0)                       #posterior inclus. probs
  lambda.post[jj,2] = mean(sqrt(reg.blas$lambda2[1001:11000])) #poster mean of lambda


### HS
  # monomvn R package is required
  reg.blas <- blasso(x, y, T=11000, RJ = FALSE, case="hs", icept=FALSE, verb=0, normalize=FALSE)
  beta.post = apply(reg.blas$beta[1001:11000,],2,mean)  # posterior est.
  Yhat = X%*%beta.post
  rmse.y[jj,5] = mean((Y-Yhat)^2)                       #RMSE of fitted values
  rmse.b[jj,5] = mean((beta.post-btrue)^2)              #RMSE of beta est.
  
### HS with RJ
  # monomvn R package is required
  reg.blas <- blasso(x, y, T=11000, RJ = TRUE, case="hs", icept=FALSE, verb=0, normalize=FALSE)
  beta.post = apply(reg.blas$beta[1001:11000,],2,median)# posterior est.
  Yhat = X%*%beta.post
  rmse.y[jj,6] = mean((Y-Yhat)^2)                       #RMSE of fitted values
  rmse.b[jj,6] = mean((beta.post-btrue)^2)              #RMSE of beta est.
  s = summary(reg.blas, burnin=1000)
  gamma4[jj,] =  as.vector(s$bn0)                       #posterior inclus. probs
  lambda.post[jj,3] = mean(sqrt(reg.blas$lambda2[1001:11000])) #poster mean of lambda

 #### rescaled spikeslabGAM
  # spikeSlabGAM R package is required
  d = data.frame(X)
  d$Y = with(d,Y)
  mcmc <- list(nChains = 1, burnin = 1000, chainLength = 10000, thin = 1, reduceRet = TRUE)
  m1 <- spikeSlabGAM(Y ~ lin(X1)+lin(X2)+lin(X3)+lin(X4)+lin(X5)+lin(X6)+lin(X7)+lin(X8)+lin(X9)+lin(X10)+lin(X11)+lin(X12)+lin(X13)+lin(X14)+lin(X15), data = d, mcmc = mcmc)
  beta.post = as.numeric(m1$postMeans$beta[-1,1]) 
  a=2*sqrt(diag(t(X)%*%X)/n)
  beta.post = beta.post/a                            # posterior est.
  Yhat = X%*%beta.post
  rmse.y[jj,5] = mean((Y-Yhat)^2)                    # RMSE of fitted values
  rmse.b[jj,5] = mean((beta.post-btrue)^2)           # RMSE of beta est.
  s=summary(m1)
  gamma5[jj,] = as.numeric(s$postMeans$pV1[,1])      #posterior inclus. probs
}
#=============================================================================