QAEnsemble (version 1.0.0)

psrfdiagnostic: 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: \(\hat{R} = \frac{S}{\sum_{i=1}^{K} \frac{s_i}{K}}\), where \(S\) is the distance between the upper and lower values of the \(100 (1 - \alpha)\%\) interval for the pooled samples, \(s_i\) is the distance between the upper and lower values of the \(100 (1 - \alpha)\%\) interval for the \(i^{\textrm{th}}\) chain, and \(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.

Usage

psrfdiagnostic(my_samples_burn_in, alpha = 0.05)

Value

A list object is returned that contains two vectors and one matrix: p_s_r_f_vec, L_vec, and d_matrix.

p_s_r_f_vec: this is the vector of potential scale reduction factors in order of the parameters

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

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)

Arguments

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.

alpha

the alpha value here corresponds to the 100(1 - alpha)% credible intervals to be estimated, with the default value as alpha = 0.05

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