R/bayes_approx.R

In CompRiskRel: Reliability and Competing Risks Analysis under Hybrid Censoring

#' Bayesian Inference via Gibbs Sampling for Censored Data Models
#'
#' @param log_posterior Function taking parameter vector and returning log posterior density.
#' @param init_par Initial parameter vector.
#' @param n_sim Total number of MCMC iterations.
#' @param burn_in Number of initial iterations to discard as burn-in.
#'
#' @return S3 object of class \code{bayes_fit} with parameter chains, posterior means, and credible intervals.
#' @export
#'
#' @examples
#' log_post <- function(th) {
#'   if (th[1] <= 0) return(-Inf)
#'   dexp(th[1], rate = 1, log = TRUE) +
#'     sum(dexp(c(0.5, 1.2, 0.8), rate = th[1], log = TRUE))
#' }
#' bayes_gibbs_censored(
#'   log_post, init_par = c(1.0),
#'   n_sim = 1000, burn_in = 200
#' )
bayes_gibbs_censored <- function(log_posterior, init_par, n_sim = 5000, burn_in = 1000) {
  p_dim <- length(init_par)
  chain <- matrix(0, nrow = n_sim, ncol = p_dim)
  chain[1, ] <- init_par
  
  curr_par <- init_par
  for (i in 2:n_sim) {
    for (j in 1:p_dim) {
      cand <- curr_par
      cand[j] <- stats::rnorm(1, mean = curr_par[j], sd = 0.1)
      if (cand[j] > 0) {
        log_alpha <- log_posterior(cand) - log_posterior(curr_par)
        if (!is.na(log_alpha) && log(stats::runif(1)) < log_alpha) {
          curr_par[j] <- cand[j]
        }
      }
    }
    chain[i, ] <- curr_par
  }
  
  post_chain <- chain[(burn_in + 1):n_sim, , drop = FALSE]
  post_means <- colMeans(post_chain)
  ci_lower <- apply(post_chain, 2, function(x) stats::quantile(x, 0.025))
  ci_upper <- apply(post_chain, 2, function(x) stats::quantile(x, 0.975))
  
  res <- list(
    chain = post_chain,
    post_means = post_means,
    ci_lower = ci_lower,
    ci_upper = ci_upper,
    method = "Gibbs Sampling"
  )
  class(res) <- "bayes_fit"
  return(res)
}

#' Bayesian Inference via Metropolis-Hastings MCMC Algorithm
#'
#' @param log_posterior Function evaluating log-posterior density.
#' @param init_par Initial parameter values.
#' @param n_sim Total MCMC sample size.
#' @param burn_in Discarded burn-in samples.
#' @param proposal_sd Proposal standard deviation vector.
#'
#' @return S3 object of class \code{bayes_fit}.
#' @export
#'
#' @examples
#' log_post <- function(th) {
#'   if (th[1] <= 0) return(-Inf)
#'   dexp(th[1], rate = 1, log = TRUE) +
#'     sum(dexp(c(0.5, 1.2, 0.8), rate = th[1], log = TRUE))
#' }
#' bayes_mh_censored(
#'   log_post, init_par = c(1.0),
#'   n_sim = 1000, burn_in = 200
#' )
bayes_mh_censored <- function(log_posterior, init_par, n_sim = 5000, burn_in = 1000, proposal_sd = 0.1) {
  p_dim <- length(init_par)
  if (length(proposal_sd) == 1) proposal_sd <- rep(proposal_sd, p_dim)
  
  chain <- matrix(0, nrow = n_sim, ncol = p_dim)
  chain[1, ] <- init_par
  curr_par <- init_par
  curr_lp <- log_posterior(curr_par)
  
  n_accept <- 0
  for (i in 2:n_sim) {
    cand <- curr_par + stats::rnorm(p_dim, mean = 0, sd = proposal_sd)
    cand_lp <- log_posterior(cand)
    if (!is.na(cand_lp) && !is.nan(cand_lp) && log(stats::runif(1)) < (cand_lp - curr_lp)) {
      curr_par <- cand
      curr_lp <- cand_lp
      n_accept <- n_accept + 1
    }
    chain[i, ] <- curr_par
  }
  
  post_chain <- chain[(burn_in + 1):n_sim, , drop = FALSE]
  post_means <- colMeans(post_chain)
  ci_lower <- apply(post_chain, 2, function(x) stats::quantile(x, 0.025))
  ci_upper <- apply(post_chain, 2, function(x) stats::quantile(x, 0.975))
  
  res <- list(
    chain = post_chain,
    post_means = post_means,
    ci_lower = ci_lower,
    ci_upper = ci_upper,
    acceptance_rate = n_accept / n_sim,
    method = "Metropolis-Hastings M-H Algorithm"
  )
  class(res) <- "bayes_fit"
  return(res)
}

#' Importance Sampling for Posterior Estimation
#'
#' @param log_target Function evaluating log target density (unnormalized log posterior).
#' @param proposal_pdf Density function of proposal distribution.
#' @param rproposal Random generation function for proposal distribution.
#' @param n_sim Sample size of proposal draws.
#'
#' @return List with posterior mean estimates, normalized importance weights, and effective sample size.
#' @export
#'
#' @examples
#' importance_sampling_censored(
#'   log_target = function(th) dexp(th, rate = 2, log = TRUE),
#'   proposal_pdf = function(th) dexp(th, rate = 1),
#'   rproposal = function(n) rexp(n, rate = 1),
#'   n_sim = 1000
#' )
importance_sampling_censored <- function(log_target, proposal_pdf, rproposal, n_sim = 10000) {
  draws <- rproposal(n_sim)
  if (is.vector(draws)) draws <- matrix(draws, ncol = 1)
  
  log_w <- apply(draws, 1, log_target) - log(apply(draws, 1, proposal_pdf))
  max_lw <- max(log_w, na.rm = TRUE)
  w <- exp(log_w - max_lw)
  w[is.na(w) | is.nan(w)] <- 0
  w_norm <- w / sum(w)
  
  est_mean <- colSums(draws * w_norm)
  ess <- 1 / sum(w_norm^2)
  
  res <- list(
    draws = draws,
    weights = w_norm,
    post_means = est_mean,
    effective_sample_size = ess,
    method = "Importance Sampling"
  )
  class(res) <- "bayes_fit"
  return(res)
}

#' Lindley Asymptotic Approximation for Posterior Expectations
#'
#' @param log_lik Log-likelihood function of parameter vector.
#' @param log_prior Log-prior density function of parameter vector.
#' @param init_par Maximum likelihood estimate or posterior mode.
#'
#' @return List with Lindley posterior estimates and standard errors.
#' @export
#'
#' @examples
#' lindley_approx_censored(
#'   log_lik = function(th) -sum((c(1.2, 0.8, 1.5) - th[1])^2),
#'   log_prior = function(th) dexp(th[1], rate = 1, log = TRUE),
#'   init_par = c(1.1)
#' )
lindley_approx_censored <- function(log_lik, log_prior, init_par) {
  p_dim <- length(init_par)
  eps <- 1e-4
  
  # Numerical 2nd derivatives of log-likelihood (Hessian matrix L_22)
  L22 <- matrix(0, p_dim, p_dim)
  f0 <- log_lik(init_par)
  for (i in 1:p_dim) {
    for (j in 1:p_dim) {
      p_ij <- init_par; p_ij[i] <- p_ij[i] + eps; p_ij[j] <- p_ij[j] + eps
      p_i <- init_par; p_i[i] <- p_i[i] + eps
      p_j <- init_par; p_j[j] <- p_j[j] + eps
      L22[i, j] <- (log_lik(p_ij) - log_lik(p_i) - log_lik(p_j) + f0) / (eps^2)
    }
  }
  
  # Invert negative Hessian to get covariance matrix sigma
  sigma <- tryCatch(solve(-L22), error = function(e) diag(0.01, p_dim))
  
  # First derivative of log prior
  rho <- numeric(p_dim)
  for (i in 1:p_dim) {
    p_plus <- init_par; p_plus[i] <- p_plus[i] + eps
    p_minus <- init_par; p_minus[i] <- p_minus[i] - eps
    rho[i] <- (log_prior(p_plus) - log_prior(p_minus)) / (2 * eps)
  }
  
  # Lindley expansion for posterior mean
  lindley_est <- init_par + as.vector(sigma %*% rho)
  se <- sqrt(pmax(diag(sigma), 1e-6))
  
  res <- list(
    post_means = lindley_est,
    se = se,
    sigma = sigma,
    method = "Lindley Approximation"
  )
  class(res) <- "bayes_fit"
  return(res)
}