phmc: Generalized Proximal Hamiltonian Monte Carlo Sampler and...

In pHMC: Proximal Hamiltonian Monte Carlo for Non-Smooth Bayesian Inference

Generalized Proximal Hamiltonian Monte Carlo Sampler and Estimator

Description

Fits Bayesian models and estimates parameters using the Proximal Hamiltonian Monte Carlo (p-HMC) algorithm for non-differentiable target densities as proposed by Shukla, Vats, and Chi (2025).

Usage

phmc(
  fn,
  grad_f = NULL,
  g = NULL,
  prox_fn = "l1",
  start,
  data = NULL,
  n_draws = 2000,
  burnin = floor(n_draws/2),
  thin = 1,
  epsilon = 0.01,
  L = 10,
  lambda_g = 0.01,
  M = NULL,
  tune_lambda = FALSE,
  verbose = FALSE,
  ...
)

Arguments

fn	Function. The smooth potential component f(x) or complete potential U(x) = -\log \pi(x). Must accept parameter vector x as first argument.
grad_f	Function or NULL. Gradient of the smooth potential component f(x). If NULL, finite difference approximation is automatically computed.
g	Function or NULL. The non-smooth penalty component g(x).
prox_fn	Function or character. Proximal mapping operator for g(x), or name of built-in operator ("l1", "l2", "elastic_net", "nuclear_norm", "none"). Default is "l1".
start	Numeric vector or matrix. Initial parameter values.
data	Optional dataset passed as second argument to fn, grad_f, and g.
n_draws	Integer > 0. Total number of MCMC iterations to run (default 2000).
burnin	Integer >= 0. Number of initial draws to discard as burn-in (default floor(n_draws / 2)).
thin	Integer >= 1. Thinning interval (default 1).
epsilon	Numeric scalar > 0. Leapfrog step size (default 0.01).
L	Integer >= 1. Number of leapfrog steps per proposal (default 10).
lambda_g	Numeric scalar > 0. Moreau-Yosida regularization parameter (default 0.01).
M	Mass matrix or NULL (defaults to identity matrix).
tune_lambda	Logical. If TRUE, uses phmc_tune to automatically select optimal lambda_g.
verbose	Logical. If TRUE, prints sampling progress.
...	Additional arguments passed to fn, grad_f, g, and prox_fn.

Value

An object of class "phmc", which is a list containing:

draws	A numeric matrix of class "matrix" containing MCMC parameter samples after burn-in and thinning.
estimates	A summary matrix of class "matrix" with rows corresponding to parameters and columns containing posterior Mean, MAP, Median, StdErr, 2.5% and 97.5% credible interval bounds, ESS, and ESS per second.
accept_rate	A numeric scalar of class "numeric" giving the overall Metropolis-Hastings acceptance rate (between 0 and 1).
ess	A named numeric vector of class "numeric" containing Effective Sample Size estimates for each parameter.
ess_per_sec	A named numeric vector of class "numeric" containing Effective Sample Size per second for each parameter.
logLik	A numeric scalar of class "numeric" giving the log-likelihood value evaluated at the Maximum A Posteriori (MAP) parameter estimate.
AIC	A numeric scalar of class "numeric" giving the Akaike Information Criterion value for model assessment.
BIC	A numeric scalar of class "numeric" giving the Bayesian Information Criterion value for model assessment.
DIC	A numeric scalar of class "numeric" giving the Deviance Information Criterion value for model assessment.
elapsed_time	A numeric scalar of class "numeric" giving total sampler execution time in seconds.
lambda_g	A numeric scalar of class "numeric" specifying the Moreau-Yosida regularization parameter used during sampling.
call	An object of class "call" recording the matched function call.

References

Shukla A, Vats D, Chi EC (2025). “Proximal Hamiltonian Monte Carlo.” arXiv preprint, \Sexpr[results=rd]{tools:::Rd_expr_doi("10.48550/arXiv.2510.22252")}.

Examples

set.seed(42)
y_data <- rnorm(100, mean = 1, sd = 0.5)
f_smooth <- function(x, y) 0.5 * sum((y - x)^2)
grad_f_smooth <- function(x, y) -sum(y - x)
fit <- phmc(fn = f_smooth, grad_f = grad_f_smooth,
            prox_fn = "l1", start = 0.5,
            data = y_data, lambda_g = 0.01,
            n_draws = 500)
summary(fit)

pHMC documentation built on Aug. 21, 2026, 5:18 p.m.