R/psrfdiagnostic.R

In QAEnsemble: Ensemble Quadratic and Affine Invariant Markov Chain Monte Carlo

#' Potential Scale Reduction Factor computation
#'
#'@description This function computes the potential scale reduction factor for
#'each parameter to formally test the convergence of the MCMC sampling to the
#'estimated posterior distribution which was developed by Gelman
#'and Brooks (1998). This potential scale reduction factor is based on empirical
#'interval lengths with the following formula:
#'\eqn{\hat{R} = \frac{S}{\sum_{i=1}^{K} \frac{s_i}{K}}}, where \eqn{S} is the
#'distance between the upper and lower values of the \eqn{100 (1 - \alpha)\%}
#'interval for the pooled samples, \eqn{s_i} is the distance between the upper
#'and lower values of the \eqn{100 (1 - \alpha)\%} interval for the
#'\eqn{i^{\textrm{th}}} chain, and \eqn{K} is the total number of chains used.
#'When the potential scale reduction factor is close to 1 for all the estimated
#'parameters, this indicates that the MCMC sampling converged to the estimated
#'posterior distribution for each parameter.
#'
#' @param my_samples_burn_in This parameter is a matrix of parameter samples
#' returned from the Ensemble MCMC algorithm 'ensemblealg', the matrix
#' dimensions are given by
#' (Number of parameters) x (Number of chains) x (Number of iterations - Number
#' of burn in iterations).
#' It is recommended to burn-in the parameter samples from the starting
#' iterations before running the 'psrfdiagnostic' to assess the convergence.
#' @param alpha the alpha value here corresponds to the 100(1 - alpha)% credible
#' intervals to be estimated, with the default value as alpha = 0.05
#'
#' @return A list object is returned that contains two vectors and one matrix:
#' p_s_r_f_vec, L_vec, and d_matrix.\cr
#'
#' p_s_r_f_vec: this is the vector of potential scale reduction factors in order
#' of the parameters \cr
#'
#' L_vec: this is the vector of distances between the upper and lower values of
#' the 95% interval for the pooled samples and these distances are in order of
#' the parameters \cr
#'
#' d_matrix: this is the matrix of distances between the upper and lower values
#' of the 95% interval for the samples in each of the chains, the matrix
#' dimensions are given by
#' (Number of parameters) x (Number of chains)
#' @export
#'
#' @references Brooks SP and Gelman A (1998) General methods for monitoring
#' convergence of iterative simulations. J Comp Graph Stat 7(4):434-455.
#'
#' @examples #Take 100 random samples from a multivariate normal distribution
#' #with mean c(1, 2) and covariance matrix matrix(c(1, 0.75, 0.75, 1), nrow = 2, ncol = 2)
#' #for each of four chains.
#'
#' my_samples_example = array(0, dim=c(2, 4, 100))
#'
#' for(j in 1:4)
#' {
#'   for(i in 1:100)
#'   {
#'     my_samples_example[,j,i] = solve(matrix(c(1, 0.75, 0.75, 1), nrow = 2, ncol = 2))%*%
#'     rnorm(2, mean = 0, sd = 1) +  matrix(c(1, 2), nrow = 2, ncol = 1, byrow = TRUE)
#'   }
#' }
#'
#' #The potential scale reduction factors for each parameter are close to 1
#' psrfdiagnostic(my_samples_example)$p_s_r_f_vec
psrfdiagnostic <- function(my_samples_burn_in, alpha = 0.05)
{
  my_size = dim(my_samples_burn_in)

  num_par = my_size[1]

  num_chains = my_size[2]

  num_iter_left = my_size[3]

  my_samples_burn_in_collapse = matrix(NA, nrow = num_par, ncol = num_chains*num_iter_left)

  for(i in 1:num_par)
  {
    my_samples_burn_in_collapse[i,] = as.vector(my_samples_burn_in[i,,])
  }

  p_s_r_f_vec = matrix(NA, nrow = 1, ncol = num_par)
  L_vec = matrix(NA, nrow = 1, ncol = num_par)
  d_matrix = matrix(NA, nrow = num_par, ncol = num_chains)

  for(i in 1:num_par)
  {
    L_vec[1,i] = stats::quantile(my_samples_burn_in_collapse[i,], probs = (1 - (alpha/2))) - stats::quantile(my_samples_burn_in_collapse[i,], probs = alpha/2)

    for(j in 1:num_chains)
    {
      d_matrix[i,j] = stats::quantile(my_samples_burn_in[i,j,], probs = (1 - (alpha/2))) - stats::quantile(my_samples_burn_in[i,j,], probs = alpha/2)
    }

    p_s_r_f_vec[1,i] = L_vec[1,i]/mean(d_matrix[i,])
  }

  return(list("p_s_r_f_vec" = p_s_r_f_vec, "L_vec" = L_vec, "d_matrix" = d_matrix))

}