## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ## ----setup-------------------------------------------------------------------- library(BayesianLasso) ## ----------------------------------------------------------------------------- # Simulate data set.seed(123) n <- 100 p <- 10 X <- matrix(rnorm(n * p), nrow = n) beta <- c(rep(2, 3), rep(0, p - 3)) y <- X %*% beta + rnorm(n) # Run the modified Hans sampler output_Hans <- Modified_Hans_Gibbs( X = X, y = y, a1=2, b1=1, u1=2, v1=1, nsamples=1000, beta_init= rep(1,10), lambda_init=1, sigma2_init=1, verbose=100) colMeans(output_Hans$mBeta) # Run the modified PC sampler output_PC <- Modified_PC_Gibbs( X = X, y = y, a1=2, b1=1, u1=2, v1=1, nsamples=1000, lambda_init=1, sigma2_init=1, verbose=100) colMeans(output_PC$mBeta) ## ----eval=FALSE--------------------------------------------------------------- # effective_sample_size <- function(samples) # { # # Using the spectral method # N <- length(samples) # # # Compute spectral density at zero using spectrum0.ar (uses AR smoothing) # spectral_density_zero <- spectrum0.ar(samples)$spec # # # Compute effective sample size using variance ratio # sample_variance <- var(samples) # tau <- spectral_density_zero / sample_variance # Integrated autocorrelation time # ESS <- N / tau # if (ESS > N) { # ESS = N # } # # return(ESS) # } ## ----generate-table-1, eval=FALSE--------------------------------------------- # # Load libraries # library(BayesianLasso) # # This code requires bayesreg # # if (requireNamespace("monomvn", quietly = TRUE)) { # monomvn::monomvn(...) # use the function safely # } else { # message("monomvn package not installed; skipping example.") # } # if (requireNamespace("bayeslm", quietly = TRUE)) { # bayeslm::bayeslm(...) # use the function safely # } else { # message("bayeslm package not installed; skipping example.") # } # if (requireNamespace("rstan", quietly = TRUE)) { # rstan::rstan(...) # use the function safely # } else { # message("rstan package not installed; skipping example.") # } # if (requireNamespace("bayesreg", quietly = TRUE)) { # bayesreg::bayesreg(...) # use the function safely # } else { # message("bayesreg package not installed; skipping example.") # } # # # # Example: dataset_name <- "diabetes2" # # if (dataset_name=="diabetes2") # { # if (!requireNamespace("lars", quietly = TRUE)) { # install.packages("lars") # } # data(diabetes, package = "lars") # # y = diabetes$y # x = diabetes$x # inds = 1:ncol(x) # # # Normalizing and scaling the dataset by function normalize() # norm = normalize(y,x, scale = TRUE) # x = norm$mX # x <- model.matrix(~.^2, data=data.frame(x=x))[,-1] # y <- norm$vy # # } # # if (dataset_name=="Kakadu2") # { # if (!requireNamespace("Ecdat", quietly = TRUE)) { # install.packages("Ecdat") # } # dat <- data("Kakadu", package = "Ecdat") # # # # Get y vector and X matrix # y <- as.vector(dat$income) # x <- dat[,c(2:21,23)] # # x <- model.matrix(~.^2,data=x)[,-1] # } # # # Make sure Crime.csv (or comData.Rdata) is downloaded manually before running this # if (dataset_name=="Crime") # { # rdata_path <- system.file("extdata", "comData.Rdata", package = "BayesianLasso") # load(rdata_path) # # Show only if file exists # crime_path <- "Crime.csv" # if (file.exists(crime_path)) { # crime_data <- read.csv(crime_path) # } else { # message("Please download 'communities.data' from the UCI Machine Learning Repository and convert to Crime.csv.") # } # # # mX <- t(t(X) %>% na.omit()) # mX <- mX[,-which(colnames(X)%in%c("ownHousQrange","rentUpperQ"))] # varnames <- colnames(mX) # vy <- Y[,"murders"] # x <- mX # y <- vy # inds = 1:ncol(x) # } # # # # Set prior hyperparameter constants # a1 = 1.0E-2 # Prior shape for sigma2 # b1 = 1.0E-2 # Prior scale for sigma2 # u1 = 1.0E-2 # Prior shape for lambda2 # v1 = 1.0E-2 # Prior scale for lambda2 # # # Initial values for lambda2 and sigma2 # lambda2_init = 10 # lambda_init = sqrt(lambda2_init) # sigma2_init = 1 # # # Number of samples to run the MCMC # nburn = 1000 # nsamples = 5000 # inds_use = (nburn + 1):nsamples # N = length(inds_use) # # # To store elapsed time and results of the PC sampler # vtime_val_PC = c() # results_PC <- NULL # to be initialized after the first run # # # Running the modified PC sampler 5 times and taking the average of the results across runs. # for (i in 1:5) { # time_val <- system.time({ # res_PC <- Modified_PC_Gibbs( # X = x, y = y, a1, b1, u1, v1, # nsamples, # lambda_init = lambda_init, # sigma2_init = sigma2_init, # verbose = 1000 # ) # })[3] # # vtime_val_PC[i] <- time_val # # # Initialize accumulators after first run # if (is.null(results_PC)) { # results_PC <- list( # mBeta = res_PC$mBeta, # vsigma2 = res_PC$vsigma2, # vlambda2 = res_PC$vlambda2 # ) # } else { # results_PC$mBeta <- results_PC$mBeta + res_PC$mBeta # results_PC$vsigma2 <- results_PC$vsigma2 + res_PC$vsigma2 # results_PC$vlambda2 <- results_PC$vlambda2 + res_PC$vlambda2 # } # } # # # Take averages # mBeta = (res_PC$mBeta)/5 # vsigma2 = (res_PC$vsigma2)/5 # vlambda2 = (res_PC$vlambda2)/5 # time_val_PC = mean(vtime_val_PC) # # # Compute effective sample sizes # vESS <- c() # for(j in 1:p) { # vESS[j] <- effective_sample_size(mBeta[inds_use,j]) # } # Ef_PC = median(vESS)/time_val_PC # # ESS_sigma2_PC = effective_sample_size(as.vector(vsigma2[inds_use])) # Ef_sigma2_PC = ESS_sigma2_PC/time_val_PC # # ESS_lambda2_PC = effective_sample_size(as.vector(vlambda2[inds_use])) # Ef_lambda2_PC = ESS_lambda2_PC/time_val_PC # # stat_vec_PC = c( # 100*median(vESS)/N, # Ef_PC, # 100*ESS_sigma2_PC/N, # Ef_sigma2_PC, # 100*ESS_lambda2_PC/N, # Ef_lambda2_PC, # time_val_PC) # # name_vals = c("mix_beta", "eff_beta", "mix_sigma2", "eff_sigma2", "mix_lambda2", "eff_lambda2", "time") # names(stat_vec_PC) = name_vals # # # # To store elapsed time and results of the Hans sampler # vtime_val_Hans = c() # results_Hans <- NULL # to be initialized after the first run # # # Running the modified Hans sampler 5 times and taking the average of the results across runs. # for (i in 1:5) { # time_val <- system.time({ # res_Hans <- Modified_Hans_Gibbs( # X = x, y = y, a1, b1, u1, v1, # nsamples, # beta_init = as.vector(colMeans(mBeta)), # lambda_init = lambda_init, # sigma2_init = sigma2_init, # verbose = 1000 # ) # })[3] # # vtime_val_Hans[i] <- time_val # # # Initialize accumulators after first run # if (is.null(results_Hans)) { # results_Hans <- list( # mBeta = res_Hans$mBeta, # vsigma2 = res_Hans$vsigma2, # vlambda2 = res_Hans$vlambda2 # ) # } else { # results_Hans$mBeta <- results_Hans$mBeta + res_Hans$mBeta # results_Hans$vsigma2 <- results_Hans$vsigma2 + res_Hans$vsigma2 # results_Hans$vlambda2 <- results_Hans$vlambda2 + res_Hans$vlambda2 # } # } # # # Take averages # mBeta = (results_Hans$mBeta)/5 # vsigma2 = (results_Hans$vsigma2)/5 # vlambda2 = (results_Hans$vlambda2)/5 # time_val_Hans = mean(vtime_val_Hans) # # # Compute effective sample sizes # vESS <- c() # for(j in 1:p) { # vESS[j] <- effective_sample_size(mBeta[inds_use,j]) # } # Ef_Hans = median(vESS)/time_val_Hans # # ESS_sigma2_Hans = effective_sample_size(as.vector(vsigma2[inds_use])) # Ef_sigma2_Hans = ESS_sigma2_Hans/time_val_Hans # # ESS_lambda2_Hans = effective_sample_size(as.vector(vlambda2[inds_use])) # Ef_lambda2_Hans = ESS_lambda2_Hans/time_val_Hans # # stat_vec_Hans = c( # 100*median(vESS)/N, # Ef_Hans, # 100*ESS_sigma2_Hans/N, # Ef_sigma2_Hans, # 100*ESS_lambda2_Hans/N, # Ef_lambda2_Hans, # time_val_Hans) # # name_vals = c("mix_beta", "eff_beta", "mix_sigma2", "eff_sigma2", "mix_lambda2", "eff_lambda2", "time") # names(stat_vec_Hans) = name_vals