I have implemented a Gibbs Sampler for the Bayesian Elastic Net (BEN) according to this paper on Penalized Regression by Kyung et al.
In this paper, they execute a simulation study that has been used in other papers on Penalized Regression (LASSO, Bridge, Ridge) to compare the performance of the proposed models.

Here are details of the simulation taken from the above mentioned paper:

We simulate data from the true model $$ y=X\beta+\sigma\epsilon \quad\epsilon_i\,{\raise.17ex\hbox{$\scriptstyle\sim$}}\,\text{iid}\,N(0,1) $$ We simulate data sets with $n=20$ to fit models and $n=200$ to compare prediction errors of proposed models with eight predictors. We let $\beta=(3,1.5,0,0,2,0,0,0)$ and $\sigma=3$. The pairwise correlation between $x_i$ and $x_j$ was set to be $corr(i,j)=0.5^{|i-j|}$.
Later on they say, that for the prediction error, they calculate the average mean squared error based on 50 replications. By average they mean the median in this case.

To simulate this data and calculate the MSE I've used following code in R:

# Number of observations
n.train <- 20
n.test  <- 200
# Error variance
sigma <- 3
# Pairwise correlation of X
cor <- 0.5
# Number of predictors
p <- 8
# Create training and test data set (package QRM and mvtnorm required)
Z <- equicorr(p, rho=cor)
X.train <- rmvnorm(n.train,sigma=Z)
X.test  <- rmvnorm(n.test,sigma=Z)
# Create error 
error.train <- rnorm(n.train,mean=0,sd=1)
error.test  <- rnorm(n.test,mean=0,sd=1)
# Create beta
beta.true <- c(3,1.5,0,0,2,0,0,0)
# Create both responses
Y.train <- X.train %*% beta.true + sigma*error.train
Y.test  <- X.test %*% beta.true + sigma*error.test

# Fit the training data set with the BEN Gibbs Sampler
beta.ben <- BEN(X.train,Y.train, iter=11000, burn = 1000)
# Calculate the predicted response
Y.pred   <- X.test %*% beta.ben
# Calculate the mean squared error (MSE)
MSE      <- sum((Y.train - Y.pred)^2)/n.train

My problem is that my results are not even close to comparable to the ones in the paper which makes me doubt my simulation study "setup".
As one of the authors of the paper has uploaded the Gibbs Sampler code and I could check if I did something wrong, I know that the problem doesn't lie there.

So my questions are:

  1. Does anybody have experience with this kind of simulation study and can check if I did something wrong?
  2. Is the MSE I calculate the same as the one used in the paper? In researching on this topic I found many different ways to calculate the MSE and it was also sometimes used but actually the mean squared prediction error was meant. For example the Wikipedia article on MSE alone lists three variations.

I don't need help with coding, rather more information on how this simulation is typically excecuted so I can figure out what I'm doing wrong.

  • $\begingroup$ Welcome to the site, @user27436. I can't tell if this is primarily a question asking for help w/ coding or w/ the ideas involved. If the former, this Q would be off-topic on CV (see our help page), but on-topic on Stack Overflow; if the latter, it belongs here. Can you edit to clarify? If you do want help primarily w/ coding, flag your Q & we'll migrate it for you (please don't cross-post, though). $\endgroup$ Jun 28, 2013 at 17:56
  • 1
    $\begingroup$ Hello gung, thanks for the comment. The question is primarily about the ideas involved and not the coding as far as I can see. I edited the question to make it clearer. $\endgroup$ Jul 2, 2013 at 10:10

1 Answer 1


From what I understood, the correlation matrix is a Toeplitz symmetric matrix with $0.5$ on the 2nd diagonal, $0.5^2$ on the 3rd diagonal etc... But you built your correlation matrix with function equicorr, which creates a matrix with $1$ on the diagonal and $0.5$ anywhere else.

Wouldn't it be better to generate your correlation matrix this way?:

Z <- toeplitz(cor^(0:(p-1)))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.