PaperlayerView exact revision on arXiv

Source-linked AI summary

Pareto Smoothed Importance Sampling

Aki Vehtari, Daniel Simpson, Andrew Gelman, Yuling Yao, Jonah Gabry

arXiv:1507.02646v9stat.COstat.MEstat.ML

TL;DR

Importance sampling can become unreliable when proposal-to-target ratios have heavy tails or vary greatly, especially in high dimensions. The paper introduces PSIS, replacing extreme weights with generalized-Pareto expected quantiles, and reports reduced mean square error, improved Monte Carlo error estimates, diagnostics, and convergence-rate results. Its broader scope includes a Pareto k-hat diagnostic for Monte Carlo estimators, with reliability bounded by tail behavior.

  • Problem

    Heavy-tailed or highly variable importance ratios can make importance sampling estimates unstable, sometimes with infinite variance, while finite variance alone may be inadequate in high dimensions.

  • Method

    PSIS stabilizes importance sampling by replacing the largest weights with expected order-statistic quantiles from a fitted generalized Pareto distribution.

  • Results

    The paper reports reduced mean square error, improved Monte Carlo standard error estimates, diagnostics for estimate reliability, and empirical convergence rates closely following known theoretical results.

  • Takeaways & Limitations

    The Pareto k-hat diagnostic provides a finite-sample convergence-rate assessment applicable to any Monte Carlo estimator, while PSIS stabilizes importance sampling estimates.

  • Takeaways & Limitations

    When k-hat > 1, the mean is expected not to exist and any estimate for the mean is invalid.

Abstract

from arXiv · show

Importance weighting is a general way to adjust Monte Carlo integration to account for draws from the wrong distribution, but the resulting estimate can be highly variable when the importance ratios have a heavy right tail. This routinely occurs when there are aspects of the target distribution that are not well captured by the approximating distribution, in which case more stable estimates can be obtained by modifying extreme importance ratios. We present a new method for stabilizing importance weights using a generalized Pareto distribution fit to the upper tail of the distribution of the simulated importance ratios. The method, which empirically performs better than existing methods for stabilizing importance sampling estimates, includes stabilized effective sample size estimates, Monte Carlo error estimates, and convergence diagnostics. The presented Pareto $\hat{k}$ finite sample convergence rate diagnostic is useful for any Monte Carlo estimator.

1 Introduction

Importance sampling adapts Monte Carlo expectations when direct draws from the target are unavailable, but heavy-tailed ratios can make estimates unstable. The paper proposes PSIS and Pareto k-hat diagnostics to stabilize estimates and assess convergence across settings.

  • Motivation: Importance sampling uses an auxiliary distribution g(theta) that is easier to sample from than the target distribution p(theta), which may be known only up to proportionality.The method reweights draws from g(theta) to estimate expectations under p(theta).
  • Motivation: Poor proposal distributions can produce heavy right-tailed importance ratios, unstable estimates, and sometimes infinite variance.The problem is especially acute when the proposal poorly captures the target distribution.
  • Motivation: In high dimensions, importance sampling can fail even with finite-variance ratios because a few draws may dominate the estimator.Ratios can vary by several orders of magnitude, making finite variance insufficient for reliable behavior.
  • Contributions: The paper proposes Pareto smoothed importance sampling to stabilize importance sampling estimates.PSIS is presented as a method for modifying problematic importance weights.
  • Contributions: The paper shows that PSIS retains consistency and finite variance for well-behaved self-normalized importance sampling estimators.These are described as usual properties of a well-behaved self-normalized estimator.
  • Contributions: The paper proposes a numerical Pareto k-hat diagnostic for finite-sample convergence rates and demonstrates PSIS and k-hat in low- and high-dimensional examples.The diagnostic is intended to support convergence assessment alongside stabilization.

2 Stabilizing Importance Sampling Estimates by Modifying the Ratios

The paper develops PSIS by modeling the upper tail of importance ratios and replacing extremes with Pareto-based expected order statistics, improving stability while retaining consistency and finite variance under broader conditions.

  • 2 Stabilizing Importance Sampling Estimates by Modifying the Ratios: Truncation is simulation-consistent with finite variance when E(|r_s|) < ∞, whereas unmodified importance sampling requires the stronger condition E(|r_s|^2) < ∞.The simpler truncation can nevertheless be too severe and produce larger-than-necessary finite-sample bias.
  • 2.2 Our Proposal: Pareto Smoothed Importance Sampling: PSIS adapts the modification of extreme importance ratios to the problem instead of applying Ionides’s universal truncation rule.The method replaces the largest ratios above a threshold with well-spaced values consistent with the estimated importance-ratio tail.
  • 2.1 Modeling the Tail of the Importance Ratios: The generalized Pareto distribution models the upper tail of the M largest importance ratios, with the shape parameter k indicating the number of finite fractional moments.The tail size uses M = floor(min(0.2S, 3√S)); the minimum tail proportion reduces small-sample bias, while the square-root rule improves high-S stability.
  • 2.2 Our Proposal: Pareto Smoothed Importance Sampling: Expected order statistics from the fitted generalized Pareto distribution provide low-bias, lower-variance replacements for the largest ratios.The inverse-CDF transformation approximates expected order statistics, becoming more accurate as M increases.
  • 2.2 Our Proposal: Pareto Smoothed Importance Sampling: For k ∈ (0.5, 1), replacing the largest ratios with expected order statistics yields finite variance and an error distribution converging to normal.The algorithm also warns when k-hat exceeds min(1 − 1/log10(S), 0.7), indicating likely instability or high bias; for S > 2000, the threshold is 0.7.
  • 2.2 Our Proposal: Pareto Smoothed Importance Sampling: PSIS is consistently better than ordinary importance sampling in the reported RMSE comparisons, but truncated importance sampling is slightly better in one intermediate-sample case.That exception occurs for θ = 2, corresponding to k-hat ≈ 0.5; the authors report it as the only example where TIS outperformed PSIS.

3 Using ˆk as a Diagnostic

The paper uses the estimated Pareto shape parameter ˆk to diagnose finite-sample convergence, effective sample size, and practical reliability of PSIS. A threshold near 0.7 marks a transition where bias, slow convergence, and impractically large errors become dominant, including for bounded ratios in high dimensions.

  • 3.1 PSIS Is Reliable When ˆk < 0.7: PSIS is generally reliable for moderate sample sizes when ˆk < 0.7, whereas estimates with ˆk > 0.7 often converge too slowly for practical use.The paper also provides a sample-size-specific threshold for small samples.
  • 3.3 ˆk Is a Good Diagnostic in Finite Samples and in High Dimensions: The ˆk diagnostic can reveal poor finite-sample behavior even when importance ratios are bounded and have finite variance, including failures caused by increasing dimension.In high-dimensional examples, ˆk captures effective-sample-size collapse and deteriorating convergence before asymptotic CLT behavior is reached.
  • 3.2.2 Pareto Means and IS L1 Deviation: For PSIS L1 deviation, required sample size grows as k increases beyond 0.5, while PSIS requires fewer samples than IS and becomes infeasible when k > 0.7.The comparison is supported by theoretical and empirical results for the two estimators.
  • 3.2.6 PSIS RMSE Convergence Rate Given ˆk: PSIS RMSE follows approximately S−1/2 in the CLT regime and S^(k−1) in the generalized-CLT regime, with a smooth transition near k = 0.5.The convergence efficiency begins declining before k = 0.5 and matches the theoretical prediction at k = 0.5.
  • 3.2.3 Truncated Pareto Means and PSIS RMSE: For k < 0.7, variance dominates PSIS RMSE and variance-based MCSE works well; above 0.7, bias dominates and the variance-based estimate begins to fail.The threshold reflects both increasing bias and difficulty diagnosing accuracy, not an abrupt change in finite-sample scaling.
  • 3.2.6 PSIS RMSE Convergence Rate Given ˆk: At k = 0.7, approximately 10 times more samples are needed to halve PSIS RMSE, and theoretical guidance becomes approximate when the ratio distribution is not exactly Pareto.The bulk and unobserved tail can influence convergence rates in non-Pareto ratio distributions.

4 Convergence of PSIS

The paper establishes that idealized Pareto smoothed importance sampling is simulation-consistent and has finite variance under standard conditions, while noting that its asymptotic theory may require infeasible sample sizes in practice.

  • PSIS is asymptotically simulation-consistent and has finite variance under standard conditions.
  • The asymptotic analysis treats k and σ as fixed in an idealized PSIS variant, although their estimates are consistent when M →∞ and M/S →0.
  • Reaching the asymptotic regime can require infeasible sample sizes even in a simple example, and asymptotic normality for PSIS remains an open question.
  • As S →∞ with M/S →0, the effect of smoothing eventually vanishes and ordinary importance-sampling consistency is recovered.
  • The PSIS estimator can be written as a Winsorized truncated estimator plus an additional bias-correction term.
  • Under the theorem’s mild conditions, PSIS converges in L1, its variance goes to zero, and it is asymptotically unbiased.

5 Practical Examples

The paper demonstrates PSIS in approximate Bayesian inference and leave-one-out cross-validation, showing improved estimates, practical speedups, and Pareto k-hat diagnostics for identifying unreliable cases.

  • Approximate Bayesian inference: PSIS improves logistic Gaussian process density estimation, achieving small KL-divergence with only S = 100 draws.The comparison used 1000 runs against MCMC and included split-normal approximation, IS, and PSIS.
  • Approximate Bayesian inference: The Pareto k-hat diagnostic indicates slightly higher variability and slower convergence for an unscaled normal approximation than for the split-normal approximation.With S = 10^4, the mean k-hat was 0.60 for the simple normal approximation.
  • LOO for Stack Loss Data: For stack-loss leave-one-out estimation, PSIS gives the smallest RMSE, while MCSE estimates are accurate for k-hat < 0.5 and optimistic for k-hat > 0.7.The evaluation compared IS, truncated IS, and PSIS for expected log predictive densities and leave-one-out predictive means.
  • LOO for 105 Protein Expression Data Sets: High k-hat values identify protein-expression leave-one-out cases where PSIS-LOO is unreliable; switching from Gaussian to Student-t models reduces the number of such values.The text attributes large k-hat values to substantial differences between proposal and target distributions, often associated with influential observations.
  • LOO for 105 Protein Expression Data Sets: PSIS-LOO+ combines PSIS-LOO with exact computation for cases with k-hat > 0.7, while K-fold cross-validation remains an alternative.The paper reports 352 such Gaussian-model cases and 53 Student-t-model cases in this example.

6 Discussion

The discussion presents PSIS as a stable importance-sampling variant with diagnostics, while noting its use across Monte Carlo and Bayesian workflows.

  • PSIS reduces importance-sampling mean square error by replacing the largest weights with expected generalized-Pareto order statistics.The method also improves Monte Carlo standard-error estimates and provides reliability diagnostics.
  • The Pareto ˆk diagnostic uses thresholds to distinguish accurate estimates, reliable MCSEs, expensive estimation, and invalid means.The stated thresholds include 1 − 1/log10(S), 0.7, and 1.
  • PSIS is presented as a reliable importance-sampling variant with a heuristic that signals when estimates become unreliable.The diagnostic is intended to make failure visible rather than silently returning an unreliable estimate.
  • Pareto ˆk and PSIS have been applied to cross-validation, Bayesian model combination, variable selection, sensitivity analysis, and variational approximations.The discussion lists implementations and applications including PSIS-LOO, stacking, projection predictive selection, and prior sensitivity analysis.
  • Pareto ˆk can diagnose general Monte Carlo estimates and support proposal comparison or adaptive importance-sampling decisions.The discussion also describes uses for stopping adaptation and selecting among weighting schemes.

Appendix A. Scaling of Distribution of Mean of Truncated Means

This appendix analyzes how PSIS error scales with sample size and tail shape, identifying CLT, generalized-CLT, and intermediate regimes.

  • PSIS replaces the M largest weights with expected Pareto order statistics and models their mean using truncated Pareto variables.The truncation is set to the largest expected order statistic, y = S^k.
  • At k = 0.5, the standard deviation scales approximately as (S/log(S))^-1/2.This matches the generalized central-limit scaling at the boundary.
  • S^-1/2 scaling applies when k is below approximately 0.5, while S^(k−1) scaling applies when k is above approximately 0.5 and below 1.The transition region narrows as S increases.
  • The scaling equations are approximate because they omit minor terms and PSIS-specific variance and bias differences, although experiments indicate the scaling order is unchanged.The appendix also notes that simple truncation can have larger bias than PSIS in some settings.

Appendix B. Relative Convergence Rates

This appendix translates PSIS RMSE scaling into relative convergence rates and compares analytic approximations with empirical results.

  • When k < 0.5 − 0.5/log10(S), the relative convergence rate is α = 1 under CLT scaling.Here RMSE scales as S^-1/2.
  • For k > 0.5 + 0.5/log10(S), the relative convergence rate is α = 2 − 2k.This follows by equating (S^-1/2)^α with S^(k−1).
  • Piecewise and smooth approximations match empirical results across most reported examples, with one noted exception where observed first- and second-moment rates are better than predicted.The comparison includes results for S = 10^3, S = 10^4, and S = 10^6.
  • The appendix derives a smoother convergence-rate approximation from the truncated-Pareto mean for use in computer-assisted diagnostics.The smooth approximation complements the simpler piecewise expression.
  • At k = 0.5 and finite S, the convergence rate is less than 1, approaching 1 only as S tends to infinity.

Appendix C. Proof of Theorem 1

The appendix proves convergence and variance properties for Winsorized importance sampling and PSIS under tail regularity and boundedness assumptions.

  • Winsorized importance sampling converges in L1 and is asymptotically unbiased under the stated assumptions.The proof uses order-statistic representations and bounds for the bias terms.
  • The generalized Pareto tail-domain assumption is expressed through a von Mises regularity condition at infinity.This condition supports approximating the importance-ratio tail by a generalized Pareto distribution.
  • The Winsorized estimator's variance is bounded and, under refinement, decreases at a rate depending on M, S, and k.With M = O(S^1/2), the refined variance rate is S^(k/2−1) for k < 1.
  • PSIS converges in L1, has variance tending to zero, and is consistent and asymptotically unbiased.
  • With the von Mises condition and bounded h, PSIS has L1 convergence at most O(S^-1/2) and variance decreasing at least as O(S^(k/2−1)).Under these additional conditions, PSIS is √n-consistent.
  • The experiments assume unknown normalization terms and therefore use self-normalized importance sampling.

D.1 Exponential Target and Proposal

The exponential example provides analytically tractable importance-ratio tails and shows that PSIS stabilizes estimates while retaining favorable convergence behavior across moments.

  • When λ > 1, the importance-ratio distribution is Pareto type I with shape k = 1 − 1/λ.
  • At k ≈0.66, IS has high variability, whereas TIS and PSIS produce more stable normalization estimates.TIS truncates the largest weights; PSIS reduces variability without biasing the largest weights.
  • PSISa has the smallest deviation, while PSIS has similar bias and slightly larger deviation than PSISa; IS has the smallest bias but largest deviation.
  • The average sample size required for RMSE = 0.05 across k from 0.1 to 0.95 fits the theoretical result.
  • For first and second moments, h-specific kh and MCSE estimates work well, although all methods have impractical convergence rates at larger kh values.
  • Higher λ values lead to lower convergence rates and higher kh estimates, with typical rates near 0.87 at k = 0.5 and 0.6 at k = 0.7.

D.2 Univariate Normal and Student’s t

Univariate experiments vary target-proposal tail combinations and proposal displacement, showing that PSIS adapts across cases and improves stability when importance ratios are problematic.

  • The experiments compare four target-proposal pairs spanning matched tails, thicker-tailed targets, and proposals with slightly thinner tails.
  • As sample size increases, tail draws reveal increasingly distinct empirical weight shapes across the target-proposal pairs.
  • Finite-sample kh estimates can differ from results at much larger sample sizes, while bounded ratios may still appear unbounded empirically.
  • PSIS adapts well across the tested target-proposal pairs when kh ∈(0.5, 0.7), whereas differences are negligible at very low kh and all methods fail at very high values.
  • In the t20-normal example, PSIS is more stable, has smaller RMSE than IS, and has more accurate MCSE estimates than TIS.
  • For non-Pareto ratio distributions, convergence rates remain good despite overestimated kh values for thick-tailed proposals.

D.3 Multivariate Normal and Student’s t

Multivariate experiments show that increasing dimensionality amplifies weight variability and estimation difficulty, while PSIS generally outperforms IS and TIS in tested settings.

  • The study examines isotropic multivariate target-proposal pairs as dimensions increase from D = 1 to 64, including a thicker-tailed proposal case.
  • Increasing dimensions amplify even small distributional differences and increase importance-weight variability.
  • With D = 16, or D = 32 for the final case, PSIS has the smallest RMSE in all tested pairs; TIS converges more slowly and IS is more unstable.
  • In a representative t20-MVN example, PSIS has smaller RMSE than IS and TIS and more accurate MCSE estimates than TIS.
  • All methods eventually fail as dimensionality increases, and bounded finite-variance ratios can resemble infinite-variance behavior in finite samples.
  • TIS can show sudden jumps when one extreme weight raises its truncation level enough that other large weights remain untruncated.
  • Observed convergence rates relative to h-specific kh estimates are close to the theoretical result.

Appendix E. h-specific ˆkh

The appendix evaluates h-specific kh diagnostics and favors a two-tail approach because it is conservative and performs better across the tested scenarios.

  • The h-specific kh diagnostic is the maximum estimate from the smallest and largest values of h(θ)r(θ), representing the left and right tails.
  • Selecting the larger of two uncertain estimates creates positive bias when the estimates are close, which the authors accept as conservative.
  • A single diagnostic based on 1 + h(θ)^2 r(θ) can mix two tails and underestimate kh in finite samples.
  • Experiments found that the two-tail diagnostic for h(θ)r(θ) worked slightly better across different scenarios.
  • The authors recommend h-specific kh but not smoothing h(θ)r(θ), because separate smoothing adds normalization bias and can noticeably increase estimation error.

Appendix F. Marginal Distribution of k

The appendix extends the profile-posterior approach from point estimation to approximate marginal inference for k, while addressing uncertainty and small-sample variability. It also adapts tail-size selection to dependent MCMC draws using effective sample size.

  • Marginal posterior of k: Quadrature points are placed along the profile curve according to the prior, with normalized likelihood-based weights used to compute the profile posterior.The profile likelihood selects k along the curve that maximizes the likelihood for each θ.
  • Marginal posterior of k: The profile posterior can approximate the marginal posterior of k, extending Zhang and Stephens's point estimates of k and σ.The approximation is compared with the exact marginal posterior computed by adaptive quadrature.
  • Uncertainty assessment: Figure 27 compares profile and exact posterior marginals across tail sample sizes M = 20, 100, and 1000 and GPD shape values k ∈ (0.3, 0.5, 0.7, 0.9).Rows vary simulated-data size, while columns vary the generating GPD value of k.
  • Regularization: For small S, additional regularization shrinks the estimate k̂ toward 0.5, trading bias for lower variance and RMSE in PSIS estimates.With S = 100 and M = 20, regularization is strong; with S = 1000 and M = 95, it has approximately one-tenth the data weight.
  • MCMC adjustment: For dependent MCMC draws, the method increases the number of tail ratios using the estimated relative efficiency and split-chain effective sample size.The effective sample size used is related to, but not directly the effective sample size of, k̂.

Appendix I. Stan Program for Linear Regression on the Stack Loss Data

This appendix gives a Stan program for Bayesian linear regression on standardized predictors and response-scale variance. The program defines regression parameters, assigns priors, models the outcome normally, and records pointwise log likelihoods.

  • Data preparation: The program standardizes each predictor using its column mean and standard deviation before fitting the regression.It stores the standardized design matrix z alongside the predictor summaries.
  • Parameters and transformed quantities: The model defines an intercept, regression coefficients, positive residual variance, and a positive hierarchical scale parameter.The residual standard deviation is computed as sigma = sqrt(sigmasq), and the mean is mu = beta0 + z * beta.
  • Probabilistic model: The model uses normal, Cauchy, hierarchical normal, and inverse-gamma priors, with normally distributed observations centered on the regression mean.The response model is y ~ normal(mu, sigma).
  • Generated quantities: Generated quantities store each observation's pointwise normal log likelihood for downstream model assessment.The program computes log_lik[i] from y[i], mu[i], and sigma.
Loading 1507.02646v9…