# -----------------------------------------------------------------------------------
# SSVS Setup
# -----------------------------------------------------------------------------------
ss0 <- 0.01
ss1 <- 100

# coef_store resulting from estimation without shrinkage
s0 <- as.numeric(ss0*apply(coef_store,2,sd))
s1 <- as.numeric(ss1*apply(coef_store,2,sd))

# initial values
delta_draw <- matrix(1,K,1)
V_pr <- diag(as.numeric(delta_draw*s1 + (1-delta_draw)*s0))
V_pr.inv <- solve(V_pr)

sigma2_draw <- 10

# store parameters
coef_store <- matrix(NA, nrow = nsave, ncol = K)
sigma2_store <- matrix(NA, nrow = nsave, ncol = 1)
delta_store <- matrix(NA,nsave,K)

# start Gibbs sampling
t.start <- Sys.time()

# -----------------------------------------------------------------------------------
# SSVS MCMC
# -----------------------------------------------------------------------------------

for(irep in 1:ntot){
  # --------------------------------------------------------------------------------
  # draw p(theta | .)
  V_po <- try(solve(XpX / sigma2_draw + V_pr.inv), silent=TRUE)
  if (is(V_po,"try-error")) V_po <- ginv(XpX / sigma2_draw + V_pr.inv)
  theta_po <- matprod(V_po, (crossprod(X, y) / sigma2_draw))

  theta_draw <- try(theta_po + t(chol(V_po)) %*% rnorm(K,0,1), silent=T)
  if (is(theta_draw, "try-error")) theta_draw <- mvrnorm(1, theta_po, V_po)
  Xtheta <- matprod(X,theta_draw)

  # --------------------------------------------------------------------------------
  # draw p(sigma2 | .)
  a_po <- a_pr + N/2
  b_po <- b_pr + crossprod(y - Xtheta)/2
  sigma2_draw <- 1/rgamma(1, a_po, b_po)

  # --------------------------------------------------------------------------------
  # X: SSVS prior
  for(j in 1:K){
    #Xa: check if low or high s fits better to theta_draw
    p0 <- dnorm(theta_draw[j], 0, sqrt(s0[j])) * 0.5
    p1 <- dnorm(theta_draw[j], 0, sqrt(s1[j])) * 0.5
    prob.delta <- p1/(p0+p1)
    if(is.nan(prob.delta)){prob.delta <- 0}

    #Xb: assign value to delta
    if(prob.delta > runif(1)) delta_draw[j] <- 1 else delta_draw[j] <- 0
    # delta_draw[j] <- sample(0:1,1,prob=c(1-prob.delta,prob.delta))
  }

  #Xc: re-estimate VCV with new deltas
  V_pr <- diag(as.numeric(delta_draw*s1 + (1-delta_draw)*s0))
  V_pr.inv <- solve(V_pr)

  # --------------------------------------------------------------------------------
  # store draws
  if(irep > nburn){
    coef_store[irep-nburn,] <- theta_draw
    sigma2_store[irep-nburn] <- sigma2_draw
    delta_store[irep-nburn,] <- delta_draw
  }
  pb$tick()
}

storage_ssvs <- list("coef_store"=coef_store,"sigma2_store"=sigma2_store,"time"=as.numeric(difftime(Sys.time(),t.start,units = "secs")))