# install.packages(c("lars", "mcmcse"))  # if needed

# ----- load data -----
data(diabetes, package = "lars")

# diabetes disease progression after one year
y <- scale(diabetes$y, center = TRUE, scale = FALSE) |> as.numeric()

# possible explanatory attributes
x <- scale(diabetes$x, center = TRUE, scale = TRUE)
n <- nrow(x)
p <- ncol(x)
colnames(x) <- colnames(diabetes$x)

######################################################################
# Bayesian lasso: MH-within-Gibbs skeleton
#   y | beta, sig2 ~ N(x beta, sig2 I)
#   beta_j iid ~ Laplace(0, lambda) with f(b) = (lambda/2) exp{-lambda |b|}
#   sig2 ~ Inv-Gamma(a, b) with density ∝ (sig2)^-(a+1) exp{-b/sig2}
#
# Gibbs:  sig2 | beta, y is conjugate
# MH:     beta_j | beta_-j, sig2, y
######################################################################

set.seed(1)

# ----- hyperparameters -----
a <- 2
b <- 2
lambda <- 1

# ----- MCMC settings -----
n_iter <- 20000
burn   <- 2000
thin   <- 1
s_beta <- 0.1   # RW step size for beta_j proposals (tune)

# ----- initialize -----
beta_curr <- rep(0, p)
sig2_curr <- 1

# ----- full conditional for sig2 -----
rcond_sig2 <- function(beta, a, b) {
  mu_hat <- as.numeric(x %*% beta)
  rss <- sum((y - mu_hat)^2)
  shape <- a + n / 2
  rate  <- b + 0.5 * rss
  1 / rgamma(1, shape = shape, rate = rate)
}

# ----- FILL IN: proposal + log acceptance ratio for ONE beta_j update -----
rprop_beta_j <- function(beta_j_curr, s_beta) {
  # return beta_j_new
}

log_r_beta_j <- function(beta_j_new, beta_j_curr, j,
                         beta_curr, sig2_curr, lambda) {
  # return log acceptance ratio for updating beta_j
}

# ----- storage -----
keep_idx <- which(seq_len(n_iter) > burn & ((seq_len(n_iter) - burn) %% thin == 0))
n_keep <- length(keep_idx)
beta_draws <- matrix(NA_real_, n_keep, p)
sig2_draws <- numeric(n_keep)

# ----- MH-within-Gibbs loop -----
k <- 0
acc_beta <- integer(p)  # coordinatewise accept counts
for (t in 1:n_iter) {
  
  # --- update beta one coordinate at a time via MH ---
  for (j in 1:p) {
    beta_j_curr <- beta_curr[j]
    beta_j_new  <- rprop_beta_j(beta_j_curr, s_beta)
    
    log_r <- log_r_beta_j(beta_j_new, beta_j_curr, j,
                          beta_curr, sig2_curr, lambda)
    
    if (log(runif(1)) < min(0, log_r)) {
      beta_curr[j] <- beta_j_new
      acc_beta[j] <- acc_beta[j] + 1L
    }
  }
  
  # --- update sig2 via Gibbs ---
  sig2_curr <- rcond_sig2(beta_curr, a, b)
  
  # --- save ---
  if (t %in% keep_idx) {
    k <- k + 1
    beta_draws[k, ] <- beta_curr
    sig2_draws[k]   <- sig2_curr
  }
}

# ----- FILL IN: coordinatewise acceptance rates -----

# ----- FILL IN: coordinatewise ess -----

# ----- FILL IN: estimates -----