[Paper Review] Simulation-efficient shortest probability intervals

Ying Liu, Andrew Gelman+1 more|arXiv (Cornell University)|Feb 8, 2013

Statistical Methods and Inference8 references4 citations

TL;DR

This paper proposes Spin (Simulation-efficient Shortest Probability Intervals), a bootstrap-based quadratic programming method to compute more stable and accurate highest posterior density (HPD) intervals from MCMC simulations. By optimally weighting order statistics, Spin reduces Monte Carlo error compared to empirical shortest intervals, especially with low effective sample sizes, while maintaining consistency and outperforming central intervals in asymmetric posteriors.

ABSTRACT

Bayesian highest posterior density (HPD) intervals can be estimated directly from simulations via empirical shortest intervals. Unfortunately, these can be noisy (that is, have a high Monte Carlo error). We derive an optimal weighting strategy using bootstrap and quadratic programming to obtain a more compu- tationally stable HPD, or in general, shortest probability interval (Spin). We prove the consistency of our method. Simulation studies on a range of theoret- ical and real-data examples, some with symmetric and some with asymmetric posterior densities, show that intervals constructed using Spin have better cov- erage (relative to the posterior distribution) and lower Monte Carlo error than empirical shortest intervals. We implement the new method in an R package (SPIn) so it can be routinely used in post-processing of Bayesian simulations.

Motivation & Objective

• To address the high Monte Carlo error of empirical shortest probability intervals when posterior simulations have low effective sample sizes.
• To develop a computationally efficient and consistent method for estimating HPD intervals from MCMC output, particularly for asymmetric and bounded posterior distributions.
• To provide a practical alternative to central intervals that better captures the mode and highest-density regions in skewed posteriors.
• To improve the reliability of posterior interval summaries in routine Bayesian data analysis without increasing simulation cost.

Proposed method

• Uses empirical order statistics from posterior simulations as input for constructing shortest intervals.
• Applies bootstrap resampling to estimate the sampling variability of empirical shortest intervals.
• Employs quadratic programming to find optimal weights for order statistics that minimize estimation error.
• Derives a weighted average of order statistics to produce a more stable and accurate shortest interval estimate.
• Sets the bandwidth parameter b = √n by default, balancing stability in high- and low-density regions.
• Implements the method in an R package (SPIn) for routine post-processing of Bayesian simulations.

Experimental results

Research questions

• RQ1 Can a weighted combination of order statistics reduce Monte Carlo error in shortest probability interval estimation from MCMC simulations?
• RQ2 How does the performance of the proposed Spin method compare to empirical shortest intervals and central intervals in terms of coverage and mean squared error?
• RQ3 Does the method maintain consistency and stability across symmetric and asymmetric posterior distributions?
• RQ4 Can Spin effectively handle boundary-affected posteriors, such as variance parameters bounded at zero?
• RQ5 Is the computational cost of Spin justified in settings with high MCMC burn-in or low effective sample sizes?

Key findings

• Spin intervals exhibit significantly lower Monte Carlo error than empirical shortest intervals, especially with small effective sample sizes.
• In asymmetric posteriors, Spin intervals are shorter and better capture the highest-density region than central intervals, which can exclude high-posterior-density regions near boundaries.
• For symmetric posteriors, Spin performs comparably to central intervals in terms of root mean squared error and coverage probability.
• The method is simulation-consistent, converging to the true HPD interval as the number of simulations increases.
• The default bandwidth parameter b = √n provides a robust balance between stability and noise reduction across diverse distributional shapes.
• Empirical studies and real-data examples, including hierarchical models and overdispersion parameters, confirm that Spin provides more reliable posterior summaries than standard empirical HPD intervals.