Introduction to Moreau-Yosida MCMC Importance Sampling (MYIS)

Shikhar Tyagi, Arvind Pandey, Bhupendra Singh, Vrijesh Tripathi

Overview

The MYIS package provides a distribution-independent framework for Moreau-Yosida Markov Chain Monte Carlo Importance Sampling (MY-IS) based on the theoretical paradigm established by Shukla, Vats, and Chi (2025).

In modern statistical modeling, Bayesian target posteriors \(\pi(\theta \mid \mathbf{x}) \propto \exp(-\psi(\theta))\) often exhibit non-differentiable priors (e.g., Lasso, fused Lasso, nuclear norm penalties) or light tails (e.g. Poisson random effects), precluding standard gradient-based MCMC algorithms such as MALA or HMC.

MYIS addresses this challenge by approximating non-smooth potentials \(\psi(\theta)\) with their smooth Moreau-Yosida envelopes: \[\psi_\lambda(\theta) = \inf_{\eta} \left\{ \psi(\eta) + \frac{1}{2\lambda} \|\eta - \theta\|^2 \right\}\]

Gradient-based samplers (\(\pi_\lambda\)-MALA, \(\pi_\lambda\)-HMC, or \(\pi_\lambda\)-RWM) draw samples from the smooth importance density \(\pi_\lambda(\theta) \propto \exp(-\psi_\lambda(\theta))\), while self-normalized importance weights: \[w_\lambda(\theta) = \exp\big(-(\psi(\theta) - \psi_\lambda(\theta))\big) \le 1\] re-weight samples to yield consistent, finite-variance estimators \(\hat{\theta}_n^{\text{MY}}\) and Bayesian marginal quantiles (Chen and Shao, 1999).

Basic Usage Example

Below is a simple demonstration estimating the rate parameter of an Exponential distribution from complete observations:

set.seed(123)
true_rate <- 1.8
sample_data <- rexp(100, rate = true_rate)

# User provides custom PDF function
my_pdf <- function(x, theta) {
  dexp(x, rate = theta[1])
}

# Run Moreau-Yosida MCMC Importance Sampling
fit <- my_is_estimate(
  pdf = my_pdf,
  data = sample_data,
  initial_theta = c(1.0),
  par_lower = 0.001,
  sampler = "mala",
  n_samples = 200,
  burnin = 50
)

# Print Summary
print(fit)
#> 
#> =========================================================
#>   Moreau-Yosida MCMC Importance Sampling (MY-IS) Fit     
#> =========================================================
#> 
#> Call:
#> my_is_estimate(pdf = my_pdf, data = sample_data, initial_theta = c(1), 
#>     par_lower = 0.001, sampler = "mala", n_samples = 200, burnin = 50)
#> 
#> Censoring Scheme: complete 
#> MCMC Sampler:     MALA 
#> Smoothing Lambda: 0.1 
#> Acceptance Rate:  97.5 %
#> Effective SS:     120.9 (Ratio ne/n = 0.6047 )
#> 
#> Parameter Estimates & Standard Errors:
#>         Estimate Std. Error    2.5%   97.5%
#> theta_1  1.72635    0.01681 1.42127 2.07502
#> =========================================================

Inference Under Right Censoring

MYIS supports arbitrary censoring schemes including complete, right, left, interval, Type-I, Type-II, progressive Type-II, and truncation.

set.seed(456)
n <- 80
x_obs <- rexp(n, rate = 1.2)
c_times <- rexp(n, rate = 0.8)
t_obs <- pmin(x_obs, c_times)
delta <- as.numeric(x_obs <= c_times)

fit_censored <- my_is_estimate(
  pdf = my_pdf,
  data = t_obs,
  initial_theta = c(1.0),
  censoring = "right",
  censoring_params = list(delta = delta),
  par_lower = 0.001,
  sampler = "mala",
  n_samples = 200,
  burnin = 50
)

summary(fit_censored)
#> 
#> =========================================================
#>   Moreau-Yosida MCMC Importance Sampling (MY-IS) Summary 
#> =========================================================
#> 
#> Censoring Scheme: right 
#> Sampler:          MALA 
#> Lambda:           0.2 
#> Acceptance Rate:  98.5 %
#> Effective SS:     86.8 (ne/n = 0.4338 )
#> 
#> Summary Table:
#>         Estimate Std. Error    2.5%   97.5%
#> theta_1  1.33216    0.02408 0.95092 1.74447
#> =========================================================

Diagnostic Plots

MYIS includes diagnostic plotting functions to inspect trace plots, autocorrelation, posterior density, and weight distribution:

plot(fit_censored, type = "all")

References