R/efficient.R

In pmwg: Particle Metropolis Within Gibbs

#' Check for efficient proposals if necessary
#'
#' Takes a mix proportion vector (3 x float) and the efficient proposal mu and
#' sigma. If efficient proposals are to be used (mix_proportion[3] > 0) then
#' test the efficient proposal values to see whether they are not null and
#' appropriate.
#'
#' @param efficient_mu The mu value for the efficient proposals
#' @param efficient_sig2 The sigma value for the efficient proposals
#' @param mix_proportion A vector of floats between 0 and 1 and summing to 1
#'   which give the proportion of particles to generate from the population
#'   level parameters, the individual random effects and the conditional
#'   parameters respectively
#'
#' @return nothing, stops operation on incorrect combination of parameters.
#' @keywords internal
check_efficient <- function(mix_proportion, efficient_mu, efficient_sig2) {
  if (mix_proportion[3] != 0) {
    if (is.null(efficient_mu) || is.null(efficient_sig2)) {
      stop(
        paste0(
          "Mu and sigma from efficient conditional ",
          "proposals must be provided for mix_proportion[3] > 0"
        )
      )
    }
  }
}


#' Create distribution parameters for efficient proposals
#'
#' From the existing samples, create a proposal distribution for drawing
#' efficient samples from.
#'
#' @param x The current pmwgs object
#'
#' @return A list containing the mu and sigma for the proposal distribution.
#' @keywords internal
create_efficient <- function(x) {
  proposal_means <- array(dim = c(x$n_pars, x$n_subjects))
  proposal_sigmas <- array(dim = c(x$n_pars, x$n_pars, x$n_subjects))
  for (s in 1:x$n_subjects) {
    cparms <- conditional_parms(
      s,
      extract_samples(x)
    )
    proposal_means[, s] <- cparms$cmeans
    proposal_sigmas[, , s] <- cparms$cvars
  }
  list(
    efficient_mu = proposal_means,
    efficient_sig2 = proposal_sigmas
  )
}


#' Obtain the efficent mu and sigma from the adaptation phase draws
#'
#' @param s current subject number
#' @param samples A list containing previous samples
#'
#' @return A list containing the conditional mean and variances for this subject
#' @keywords internal
conditional_parms <- function(s, samples) {
  tmudim <- dim(samples$theta_mu)
  n_par <- tmudim[1]
  n_iter <- tmudim[2]
  pts2_unwound <- apply(
    samples$theta_sig,
    3,
    unwind
  )
  all_samples <- rbind(
    samples$alpha[, s, ],
    samples$theta_mu[, ],
    pts2_unwound
  )
  mu_tilde <- apply(all_samples, 1, mean)
  sigma_tilde <- stats::var(t(all_samples))
  condmvn <- condMVNorm::condMVN(
    mean = mu_tilde,
    sigma = sigma_tilde,
    dependent.ind = 1:n_par,
    given.ind = (n_par + 1):length(mu_tilde),
    X.given = c(
      samples$theta_mu[, n_iter],
      unwind(samples$theta_sig[, , n_iter])
    )
  )
  list(cmeans = condmvn$condMean, cvars = condmvn$condVar)
}


#' Test that the sampler has successfully adapted
#'
#' @param pmwgs The full pmwgs object with all samples
#' @param n_unique The number of unique samples to look for in random effects
#'   for each subject.
#' @param i The number for the current iteration of the sampler
#'
#' @return A string representing successful/unsuccessful adaptation. Can be one
#'   of c("success", "continue", "increase")
#'
#' @keywords internal
test_sampler_adapted <- function(pmwgs, n_unique, i) {
  if (i < n_unique) {
    return("continue")
  }
  test_samples <- extract_samples(pmwgs, stage = "adapt")
  if (check_adapted(test_samples$alpha, unq_vals = n_unique)) {
    message("Enough unique values detected: ", n_unique)
    message("Testing proposal distribution creation")
    attempt <- try({
      lapply(
        X = 1:pmwgs$n_subjects,
        FUN = conditional_parms,
        samples = test_samples
      )
    })
    if (class(attempt) == "try-error") {
      warning("An problem was encountered creating proposal distribution")
      warning("Increasing required unique values and continuing adaptation")
      return("increase")
    }
    else {
      message("Successfully adapted after ", i, "iterations - stopping early")
      return("success")
    }
  }
  return("continue")
}


#' Check whether the adaptation phase has successfully completed
#'
#' @param samples The subject mean samples with which we are working
#' @param unq_vals The number of unique values for each subject
#'
#' @return A boolean TRUE or FALSE depending on the result of the test
#' @keywords internal
check_adapted <- function(samples, unq_vals = 20) {
  # Only need to check uniqueness for one parameter
  first_par <- samples[1, , ]
  # Split the matrix into a list of vectors by subject
  # Needed for the case where every sample is unique for all subjects
  first_par_list <- split(first_par, seq(NROW(first_par)))
  # Get unique pars (new accepted particles) and check length for
  # all subjects is greater than unq_vals
  all(lapply(lapply(first_par_list, unique), length) > unq_vals)
}