Package {LugsailGR}

Type: Package
Title: Generalized Gelman-Rubin Diagnostic and Effective Sample Size for MCMC
Version: 0.1.0
Date: 2026-07-26
Description: Provides generalized univariate and multivariate 'Gelman-Rubin' convergence diagnostics, effective sample size ('ESS') estimates, and principled termination thresholds for Markov chain Monte Carlo ('MCMC') simulations, based on Vats and Knudson (2021) <doi:10.1214/20-STS812>. The package incorporates replicated lugsail batch means variance estimators to construct stable convergence statistics for single and multiple chains. Additionally, it offers comprehensive tools for evaluating 'MCMC' output generated from user-supplied probability density functions ('PDF') or log-likelihoods, including implementations for censored data models under right, left, interval, 'Type-I', 'Type-II', progressive, and hybrid censoring schemes.
Depends: R (≥ 3.5.0)
License: GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
Encoding: UTF-8
LazyData: true
Imports: stats, graphics
Suggests: testthat (≥ 3.0.0)
RoxygenNote: 7.3.3
NeedsCompilation: no
Packaged: 2026-07-26 17:38:52 UTC; shikhar tyagi
Author: Shikhar Tyagi ORCID iD [aut, cre], Arvind Pandey [aut], Bhupendra Singh [aut], Vrijesh Tripathi [aut]
Maintainer: Shikhar Tyagi <shikhar1093tyagi@gmail.com>
Repository: CRAN
Date/Publication: 2026-08-05 07:40:09 UTC

Calculate Minimum Required Effective Sample Size (M_alpha_eps_p)

Description

Computes the theoretical minimum effective sample size required to obtain a confidence region with relative volume epsilon and coverage 1 - alpha, based on Vats, Flegal and Jones (2019) and Vats and Knudson (2021).

Usage

calc_ess_bound(p = 1, alpha = 0.05, epsilon = 0.1)

Arguments

p

Integer, number of parameters.

alpha

Numeric, significance level (default 0.05 for 95% confidence).

epsilon

Numeric, relative volume tolerance (default 0.10).

Value

Numeric scalar, required minimum ESS.

Examples

calc_ess_bound(p = 1, alpha = 0.05, epsilon = 0.10)
calc_ess_bound(p = 5, alpha = 0.05, epsilon = 0.05)

Calculate Replicated Lugsail Batch Means Estimator

Description

Computes the replicated lugsail batch means variance estimator (univariate) or time-average covariance matrix (multivariate) for MCMC chains as described by Vats and Knudson (2021).

Usage

calc_lugsail_bm(arr, b = NULL)

Arguments

arr

3D numeric array of dimension (n, p, m) representing n steps, p parameters, and m chains.

b

Batch size. If NULL, defaults to floor(sqrt(n)).

Value

A list containing:

tau_L

Scalar variance estimate if p = 1, or p x p covariance matrix if p > 1.

b

Batch size used.

a

Number of batches.

Calculate Gelman-Rubin Termination Threshold (delta_epsilon)

Description

Computes the principled Gelman-Rubin diagnostic termination threshold delta_epsilon based on the target ESS bound M_alpha_eps_p and number of chains m.

Usage

calc_psrf_cutoff(m = 1, M_val = 1537)

Arguments

m

Integer, number of chains.

M_val

Numeric, minimum required ESS obtained from calc_ess_bound.

Value

Numeric scalar, target threshold delta_epsilon.

Examples

M_bound <- calc_ess_bound(p = 1, alpha = 0.05, epsilon = 0.10)
calc_psrf_cutoff(m = 3, M_val = M_bound)

Gelman-Rubin Convergence Diagnostics for Censored Data Models

Description

Evaluates lugsail Gelman-Rubin convergence diagnostics, effective sample size, and parameter estimates for MCMC chains fitted to censored data models under various censoring schemes (right, left, interval, Type-I, Type-II, progressive, and hybrid).

Usage

censor_gr_test(
  pdf_fn,
  cdf_fn,
  surv_fn = NULL,
  data,
  start_par,
  censor_type = c("right", "left", "interval", "type1", "type2", "progressive", "hybrid"),
  status = NULL,
  censor_info = NULL,
  n_iter = 2000,
  n_chains = 3,
  burn_in = 0.5,
  scale = 0.5,
  alpha = 0.05,
  epsilon = 0.1
)

Arguments

pdf_fn

Function 'function(x, par)' returning probability density values for variable 'x' given parameters 'par'.

cdf_fn

Function 'function(x, par)' returning cumulative distribution function values.

surv_fn

Optional function 'function(x, par)' returning survival function values. If 'NULL', computed as '1 - cdf_fn(x, par)'.

data

Numeric vector or matrix of observed sample data.

start_par

Initial parameter vector.

censor_type

Character string specifying censoring scheme. Options include: "right" (Right censoring), "left" (Left censoring), "interval" (Interval censoring), "type1" (Type-I censoring), "type2" (Type-II censoring), "progressive" (Progressive Type-II censoring), or "hybrid" (Hybrid censoring).

status

Numeric vector or matrix indicating event/censoring status (e.g. 1 = observed, 0 = censored).

censor_info

Optional list specifying parameters for censoring schemes (e.g. cutoff time T_censor, number of failures r_censor, or removal pattern R_pattern).

n_iter

Integer, total number of MCMC iterations per chain (default 2000).

n_chains

Integer, number of parallel chains (default 3).

burn_in

Numeric, burn-in fraction (default 0.50).

scale

Numeric or vector, proposal scale for random walk Metropolis-Hastings (default 0.5).

alpha

Numeric, significance level (default 0.05).

epsilon

Numeric, relative volume tolerance (default 0.10).

Value

An object of class "lugsail_gr" containing Gelman-Rubin diagnostics, effective sample sizes, cutoff thresholds, parameter estimates, and MCMC samples.

References

Vats, D. and Knudson, C. (2021). Revisiting the Gelman–Rubin Diagnostic. Statistical Science, 36(4), 518–529. doi:10.1214/20-STS812.

Examples

# Example: Right-censored exponential distribution
set.seed(123)
n_obs <- 30
true_rate <- 0.5
times <- rexp(n_obs, rate = true_rate)
censor_time <- 2.0
obs_times <- pmin(times, censor_time)
status <- as.numeric(times <= censor_time)

f_exp <- function(x, rate) dexp(x, rate = rate)
F_exp <- function(x, rate) pexp(x, rate = rate)
S_exp <- function(x, rate) 1 - pexp(x, rate = rate)

fit <- censor_gr_test(pdf_fn = f_exp, cdf_fn = F_exp, surv_fn = S_exp,
                     data = obs_times, start_par = c(rate = 0.8),
                     censor_type = "right", status = status,
                     n_iter = 1000, n_chains = 3)
print(fit)

Format Input Chains into a 3D Array

Description

Format Input Chains into a 3D Array

Usage

format_chains(x)

Arguments

x

Input object: vector, matrix, 3D array, or list of matrices/data.frames.

Value

3D numeric array of dimension (n, p, m).

Upgraded Univariate and Multivariate Gelman-Rubin Diagnostic

Description

Computes the upgraded Gelman-Rubin convergence statistics (\hat{R}_L and \hat{R}^p_L), Effective Sample Size (\widehat{\mathrm{ESS}}), and principled termination thresholds (\delta_{\epsilon}) proposed by Vats and Knudson (2021) using replicated lugsail batch means variance estimators.

Usage

lugsail_gr(x, alpha = 0.05, epsilon = 0.1, b = NULL, multivariate = TRUE)

Arguments

x

Input MCMC chain output. Can be a numeric vector (single univariate chain), matrix of dimension (n, p) (single multivariate chain), 3D array of dimension (n, p, m) (multiple multivariate chains), or a list of length m containing matrices or data.frames.

alpha

Numeric, significance level for the confidence region (default 0.05 for 95% confidence).

epsilon

Numeric, relative volume tolerance for target precision (default 0.10).

b

Optional batch size for lugsail batch means. If NULL, defaults to floor(sqrt(n)).

multivariate

Logical, if TRUE (default), computes multivariate Gelman-Rubin diagnostic \hat{R}^p_L and multivariate ESS when p > 1.

Value

An object of class "lugsail_gr" containing:

psrf

Vector of univariate lugsail potential scale reduction factors \hat{R}_L for each parameter.

mpsrf

Multivariate lugsail potential scale reduction factor \hat{R}^p_L.

ess

Vector of univariate effective sample sizes for each parameter.

mess

Multivariate effective sample size \widehat{\mathrm{ESS}}_p.

delta_eps

Principled target threshold \delta_{\epsilon}.

M_alpha_eps_p

Minimum required effective sample size bound M_{\alpha, \epsilon, p}.

converged

Logical, indicating if diagnostic has converged (mpsrf <= delta_eps).

means

Posterior sample mean vector across chains.

sd

Posterior standard deviation vector.

s2

Sample variance (scalar or matrix S).

tau_L

Replicated lugsail variance/covariance matrix estimate \hat{T}_L.

n

Chain length (number of iterations per chain).

m

Number of parallel chains.

p

Number of parameter dimensions.

alpha

Significance level.

epsilon

Relative volume tolerance.

References

Vats, D. and Knudson, C. (2021). Revisiting the Gelman–Rubin Diagnostic. Statistical Science, 36(4), 518–529. doi:10.1214/20-STS812.

Gelman, A. and Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7(4), 457–472. doi:10.1214/ss/1177011136.

Brooks, S. P. and Gelman, A. (1998). General methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics, 7(4), 434–455. doi:10.1080/10618600.1998.10474787.

Examples

# Example 1: Univariate MCMC output (3 chains of 1000 iterations)
set.seed(123)
chain1 <- rnorm(1000, mean = 0, sd = 1)
chain2 <- rnorm(1000, mean = 0.05, sd = 1)
chain3 <- rnorm(1000, mean = -0.05, sd = 1)
res1 <- lugsail_gr(list(chain1, chain2, chain3))
print(res1)

# Example 2: Multivariate MCMC output (2 chains of 500 iterations, 2 parameters)
mat1 <- matrix(rnorm(1000), ncol = 2)
mat2 <- matrix(rnorm(1000), ncol = 2)
res2 <- lugsail_gr(list(mat1, mat2))
summary(res2)

Evaluate Gelman-Rubin Diagnostics for User-Supplied Target Density Function

Description

Runs Markov chain Monte Carlo (MCMC) sampling for a user-supplied probability density function (PDF) or log-likelihood function of variable and parameters with data, and returns upgraded lugsail Gelman-Rubin convergence test statistics, effective sample sizes, parameter estimates, and termination decisions.

Usage

mcmc_gr_test(
  target_pdf,
  data = NULL,
  start_par,
  n_iter = 2000,
  n_chains = 3,
  burn_in = 0.5,
  scale = 0.5,
  alpha = 0.05,
  epsilon = 0.1,
  ...
)

Arguments

target_pdf

Function taking parameter vector par as its first argument and optional data. It should return the probability density value or log-density value for par.

data

Optional data vector, matrix, data.frame, or list passed to target_pdf.

start_par

Initial parameter vector. If named, parameter names will be preserved.

n_iter

Integer, total number of MCMC iterations per chain (default 2000).

n_chains

Integer, number of parallel MCMC chains to generate (default 3).

burn_in

Numeric, fraction of initial burn-in samples to discard (default 0.50).

scale

Numeric or vector, proposal standard deviation for random walk Metropolis-Hastings sampler (default 0.5).

alpha

Numeric, significance level for convergence bound (default 0.05).

epsilon

Numeric, relative volume tolerance (default 0.10).

...

Additional parameters passed to target_pdf.

Value

An object of class "lugsail_gr" containing Gelman-Rubin test statistics, effective sample sizes, cutoff thresholds, parameter estimates, standard errors, and MCMC chain samples.

References

Vats, D. and Knudson, C. (2021). Revisiting the Gelman–Rubin Diagnostic. Statistical Science, 36(4), 518–529. doi:10.1214/20-STS812.

Examples

# Example 1: Normal distribution target PDF with data
set.seed(42)
true_mean <- 3
my_data <- rnorm(50, mean = true_mean, sd = 1)

# Log-posterior density function for parameter mu
log_post <- function(mu, data) {
  sum(dnorm(data, mean = mu, sd = 1, log = TRUE))
}

# Run diagnostic test
fit <- mcmc_gr_test(target_pdf = log_post, data = my_data, start_par = c(mu = 0), 
                    n_iter = 1000, n_chains = 3)
print(fit)

Plot Method for lugsail_gr Objects

Description

Produces diagnostic plots for MCMC chains, including trace plots, running Gelman-Rubin statistic \hat{R}_L plots, and density estimates.

Usage

## S3 method for class 'lugsail_gr'
plot(x, type = c("all", "trace", "running_gr", "density"), ...)

Arguments

x

An object of class "lugsail_gr".

type

Character string specifying plot type: "all" (default), "trace", "running_gr", or "density".

...

Additional graphical parameters.

Value

No return value, called for side effects (generating plots).

Examples

set.seed(123)
chain1 <- rnorm(500)
chain2 <- rnorm(500)
fit <- lugsail_gr(list(chain1, chain2))
plot(fit, type = "running_gr")

Print Method for lugsail_gr Objects

Description

Print Method for lugsail_gr Objects

Usage

## S3 method for class 'lugsail_gr'
print(x, ...)

Arguments

x

An object of class "lugsail_gr".

...

Further arguments passed to print.

Value

Invisibly returns x.

Summary Method for lugsail_gr Objects

Description

Summary Method for lugsail_gr Objects

Usage

## S3 method for class 'lugsail_gr'
summary(object, ...)

Arguments

object

An object of class "lugsail_gr".

...

Further arguments passed to summary.

Value

A summary table object.

Sample Titanic Survival Dataset

Description

Subset of the Titanic passenger dataset used in Vats and Knudson (2021) for demonstrating multivariate Gelman-Rubin convergence diagnostics in Bayesian logistic regression.

Usage

titanic_sub

Format

A data frame with 100 rows and 6 variables:

Survived

Binary indicator of survival (1 = survived, 0 = died)

Pclass

Passenger class (1 = 1st, 2 = 2nd, 3 = 3rd)

Sex

Sex of passenger (1 = female, 0 = male)

Age

Age in years

SibSp

Number of siblings/spouses aboard

Fare

Passenger fare

Source

Vats, D. and Knudson, C. (2021). Revisiting the Gelman–Rubin Diagnostic. Statistical Science, 36(4), 518–529. doi:10.1214/20-STS812.

mirror server hosted at Truenetwork, Russian Federation.