bayes_lasso.R

library('runjags')


lasso.bayes <- function(y, X, n.iter = 10000) {
    ms <- '
    model {
        sigma ~ dunif(0, 100)
        tau <- pow(sigma, -2)

        lambda ~ dunif(0, 10)
        b[1] ~ dnorm(0, .001)

        for (i in 2:p) {
            b[i] ~ ddexp(0, lambda)
        }
        mu <- X %*% b

        for (i in 1:n) {
            y[i] ~ dnorm(mu[i], tau)
        }
    }' # for the paper, see Park & Casella (2004)
    params <- c('b', 'sigma', 'lambda')

    jags.dat <- list('n' = length(y), 'y' = y, 'X' = X, 'p' = ncol(X))
    run.jags(ms, monitor = params, data = jags.dat, n.chains = 4, method = "rjparallel")
}

y <- rnorm(100)
X <- cbind(1, rnorm(100), runif(100, 0, 100), rbeta(100, 1, 9))

ss <- lasso.bayes(y, X, n.iter = 100)