PaperlayerView exact revision on arXiv

Source-linked AI summary

GetDist: a Python package for analysing Monte Carlo samples

Antony Lewis

arXiv:1910.13970v2astro-ph.IMastro-ph.COphysics.data-an

TL;DR

Monte Carlo density estimation must remain accurate with weighted or correlated samples and hard prior boundaries. GetDist addresses this with boundary-corrected, bias-reduced KDE and automatic bandwidth selection, together with tools for analysis and visualization. The package is broadly useful beyond cosmology, while its current implementation assumes stationarity and does not fully handle dynamic-sampling correlations.

  • Problem

    Accurate density estimation from Monte Carlo samples must account for smoothing, hard prior boundaries, correlations, and weights.

  • Method

    GetDist combines linear boundary kernels, multiplicative bias correction, correlation- and weight-aware effective sample sizes, and automatic bandwidth selection within an analysis and visualization package.

  • Results

    GetDist provides robust density estimation and publication-quality analysis tools for cosmological and broader low-dimensional sample-density applications.

  • Takeaways & Limitations

    The package supports fast exploration and reproducible scripted analysis of constrained Monte Carlo parameter densities and figures.

  • Takeaways & Limitations

    The current implementation assumes stationarity and treats samples from nested and other dynamic sampling methods as independent.

Abstract

from arXiv · show

Monte Carlo techniques, including MCMC and other methods, are widely used in Bayesian inference to generate sets of samples from a parameter space of interest. The Python GetDist package provides tools for analysing these samples and calculating marginalized one- and two-dimensional densities using Kernel Density Estimation (KDE). Many Monte Carlo methods produce correlated and/or weighted samples, for example produced by MCMC, nested, or importance sampling, and there can be hard boundary priors. GetDist's baseline method consists of applying a linear boundary kernel, and then using multiplicative bias correction. The smoothing bandwidth is selected automatically following Botev et al., based on a mixture of heuristics and optimization results using the expected scaling with an effective number of samples (defined here to account for both MCMC correlations and weights). Two-dimensional KDE uses an automatically-determined elliptical Gaussian kernel for correlated distributions. The package includes tools for producing a variety of publication-quality figures using a simple named-parameter interface, as well as a graphical user interface that can be used for interactive exploration. It can also calculate convergence diagnostics, produce tables of limits, and output in LaTeX, and is publicly available.

I. INTRODUCTION

GetDist addresses the challenges of estimating and visualizing densities from correlated, weighted Monte Carlo samples, including samples subject to hard prior boundaries. It combines specialized KDE methods with analysis, diagnostics, visualization, and interactive tools.

  • Motivation: Monte Carlo methods underpin Bayesian inference, but accurate density estimation must handle smoothing, boundaries, correlations, and weights.These challenges motivate methods beyond simple histograms and basic independent-sample KDE.
  • Package scope: GetDist provides optimized KDE, convergence diagnostics, statistical analysis, and publication-quality visualization through a Python package.The package also offers an interactive graphical interface alongside programmatic control.
  • Applications: GetDist supports correlated samples and hard prior boundaries while producing marginalized one- and two-dimensional densities and constraints.The Planck example illustrates marginalization, correlations, and boundary-prior handling in a high-dimensional parameter space.
  • Methodological contributions: Its key methods include bandwidth selection accounting for correlations and boundary effects, higher-order bias correction, and effective sample-size estimators for correlated weighted samples.These methods target the sampling properties encountered in MCMC, importance, and nested sampling.
  • Scope: The methods were developed primarily for cosmological parameter estimation but are applicable to broader low-dimensional density, constraint, and plotting problems.The paper reports use of GetDist plots and results across hundreds of published papers.

II. WEIGHTED SAMPLES

GetDist treats weighted samples as a general representation of Monte Carlo output and distinguishes their effects from correlations in MCMC chains. The resulting estimators support unbiased weighted means under specified sampling conditions.

  • Weighted estimators: Each sample point X_i is assigned a weight w_i, which may equal unity for unweighted samples.Weighted sums then estimate means of functions under the target distribution.
  • Visualization example: Figure 1 shows GetDist constructing marginalized 1D plots and 2D 68% and 95% density contours from all samples using KDE.The example includes a hard z_re > 6.5 prior and a tight H0–Ωm correlation.
  • Weighted estimators: The weighted mean estimator is unbiased when the sampling distribution and weights satisfy the paper’s stated expectation condition.The condition covers importance sampling with appropriate non-stochastic weights and integer-weighted MCMC rejection counts.
  • Sampling origins: Importance sampling and nested sampling naturally produce weighted samples, while MCMC can produce repeated points through rejected proposals.MCMC samples are also correlated, so weighting alone does not characterize their Monte Carlo noise.

III. KERNEL DENSITY ESTIMATION (KDE)

KDE estimates a continuous density by smoothing sampled points with kernels, offering a systematic alternative to histogram-based estimates. GetDist implements weighted, Gaussian-kernel KDE with bandwidth-controlled smoothness and FFT-based evaluation for large samples.

  • KDE framework: KDE places a smooth kernel at each sample point and sums the contributions to estimate a continuous probability density.This avoids the bin-width and bin-placement sensitivity of histograms.
  • KDE framework: GetDist uses slightly truncated zero-centered Gaussian kernels whose width h controls the smoothness of the estimated density.The framework also permits a covariance matrix for multivariate kernels.
  • Computation: FFT-based convolution enables fast evaluation of finely binned KDE densities for large sample sets and low-dimensional marginals.The implementation is designed for the relatively large effective sample sizes common in converged MCMC analysis.
  • Weighted KDE: Weighted KDE replaces equal sample contributions with weights w_i and uses weighted bin sums in the discretized implementation.This extends the estimator to importance and nested-sampling outputs.
  • KDE framework: As h approaches zero, the kernel estimator converges to the underlying density when the kernel tends to a delta function.Bandwidth selection therefore controls the practical bias–variance and smoothness trade-off.

A. KDE bias and linear boundary kernels

Kernel smoothing is biased near hard prior boundaries because samples are absent on the excluded side. GetDist uses linear boundary kernels and multiplicative normalization to reduce this bias, including for non-axis-aligned boundaries.

  • Boundary effects: Boundary priors bias KDE because smoothing extends into regions where no samples are allowed.The resulting density is suppressed near the boundary.
  • Illustrations: Figure 2 shows that boundary correction is essential for truncated densities, while higher-order corrections can become less stable.In the example, first-order correction outperforms zeroth-order correction, but this is not guaranteed.
  • Boundary effects: Normalizing by the allowed-region convolution removes the leading boundary bias but leaves linear bias when the boundary density has a gradient.
  • Linear boundary kernels: Linear boundary kernels use non-symmetric kernels near boundaries to remove the remaining first-order gradient bias.GetDist restricts the construction to linear kernels because higher-order multidimensional forms become cumbersome.
  • Linear boundary kernels: The residual bias is O(h2) even at the boundary, and the correction can incorporate arbitrary prior boundaries not aligned with parameter axes.
  • Practical correction: The linear boundary estimator is not guaranteed to remain positive, so GetDist imposes positivity and renormalizes the density.

B. Statistical and total error

KDE error reflects a trade-off between smoothing bias and statistical noise. The optimal bandwidth depends on effective sample size, dimension, kernel properties, and density curvature, with the stated scaling derived for independent samples.

  • Error decomposition: Broader kernels reduce statistical noise but increase bias, whereas narrower kernels do the opposite.
  • Error decomposition: For small h, the statistical term dominates the error expansion.
  • Error decomposition: The total mean integrated error combines bias and statistical terms, motivating the Asymptotic Mean Integrated Squared Error criterion.
  • Bandwidth selection: The asymptotically optimal bandwidth scales with effective sample number as N_eff^-1/(4+d).The scaling depends on the dimensionality d.
  • Bandwidth selection: Greater density curvature requires smaller bandwidth because broader smoothing removes more structure.This result is specific to the simple linear estimator assuming independent sample points.

C. Multiplicative bias correction

Boundary-corrected kernels remove leading bias but can still broaden peaks and overestimate errors. Multiplicative bias correction flattens the density before convolution, achieving higher-order accuracy away from boundaries while retaining smooth estimates for smaller samples.

  • Motivation: Boundary kernels leave systematic peak broadening that can produce overestimated errors unless h is much smaller than the density width.
  • Method: Multiplicative bias correction uses an estimate of the density shape to make the convolved quantity nearly flat.
  • Results: The corrected estimator has bias O(h4) away from boundaries when the density is sufficiently smooth.
  • Trade-offs: Multiplicative correction produces smooth density estimates with relatively small sample sizes, but that smoothness can mask sampling uncertainty.

D. Correlated samples

Correlations and weights both affect the statistical uncertainty of KDE estimates, although correlations have limited impact on the optimal bandwidth. GetDist therefore uses an effective sample number tailored to the parameter subspace and local correlation structure.

  • Correlation effects: MCMC correlations increase density-estimation uncertainty for a fixed number of samples, while weights also contribute to sampling noise.
  • Bandwidth selection: Correlations do not strongly change the optimal bandwidth because their main effect is on statistical uncertainty rather than bandwidth selection.
  • Effective sample size: The KDE variance depends on the number of sample-point sequences lying within bandwidth h, as measured by a local K∗K filter.
  • Scope limitation: For nested and other dynamic sampling outputs, GetDist currently treats samples as independent, a limitation identified for future improvement.
  • Effective sample size: The effective sample number is defined for the selected parameter subspace and incorporates weights, duplicated samples, and short-range MCMC correlations.
  • Implementation: The effective KDE sample number depends on h, so GetDist evaluates it at a fiducial bandwidth near h ≈ 0.2σ.

E. Choice of kernel bandwidth

GetDist selects kernel bandwidths by balancing smoothing bias and sampling variance, using plug-in methods, effective sample scaling, and fallbacks for difficult distributions. Tests compare these choices across one- and two-dimensional Gaussian mixtures and show improved error over basic estimators in most cases.

  • Bandwidth motivation: Bandwidth selection balances smoothing bias against statistical variance, with the basic one-dimensional scaling h ∝ N^-1/5 and convergence ∝ N^-4/5.A bandwidth that is too broad washes out features, while one that is too narrow exposes sampling noise.
  • Bandwidth estimation: Plug-in bandwidth methods estimate unknown density derivatives and can perform especially well for multimodal distributions.Simple scale rules are retained mainly as fallbacks and for setting fiducial scales.
  • Bias correction: Multiplicative bias correction changes the scaling to h ∝ N^-1/9 with overall convergence ∝ N^-8/9 away from boundaries.The higher-order bias produces an error form proportional to Ah^8 + B/(Nh).
  • Bandwidth comparison: For effective sample sizes of order 1000, the selected smoothing width is about twice as broad as the basic estimator and has lower statistical noise.The effective sample count accounts for the sampling information relevant to the bandwidth choice.
  • Two-dimensional tests: The two-dimensional tests compare 68% and 95% true-density contours with estimated contours and evaluate normalized average integrated squared error.The two-dimensional bandwidth may be estimated from separate one-dimensional results, although GetDist also offers a two-dimensional expression.
  • Empirical tests: Across 500 simulations of 10000 samples, the default multiplicative-bias-corrected linear-boundary estimator usually has lower error than the basic Parzen estimator.Higher-order correction can be better in some cases but is more sensitive to bandwidth optimization; the default is suboptimal for Gaussian shapes whose leading bias vanishes.

Multivariate bandwidth matrix

GetDist constructs two-dimensional KDE bandwidths that adapt to sample correlations, boundaries, and bias-correction scaling. Its strategy combines covariance-based rotations, fixed-point or fallback rules, and numerical safeguards, but has limitations for complex correlations and boundaries.

  • Kernel shape: The two-dimensional kernel is generally elliptical, matching correlation directions to smooth degeneracies without unnecessarily broadening well-constrained directions.A single position-independent kernel works well for one clearly correlated direction but may be suboptimal for more complex distributions.
  • Limitations: The implementation may be inaccurate for strongly elliptical kernels, and boundary treatment cannot easily approximate odd derivative terms through the DCT-based approach.When fixed-point optimization fails, the code uses a plugin bandwidth estimate; isotropic derivative evaluation can remain inaccurate for strongly elliptical optima.
  • Limitations: For distributions with boundaries in both x and y, the inability to preserve both boundaries during rotation constrains the bandwidth strategy.The general workflow includes separate cases for one boundary, no boundaries, and boundaries in both directions.
  • Bandwidth construction: GetDist uses sample-covariance Cholesky rotations when boundaries permit, preserving any boundary aligned with an unchanged coordinate.The transformed samples are scaled and binned before bandwidth estimation, then the bandwidth matrix is rotated back.
  • Boundary handling: With boundaries in both coordinates, the method avoids rotation and evaluates a diagonal bandwidth rule on untransformed samples, using a covariance-based Gaussian fallback for very high correlations.Boundary-aware fixed-point solutions can otherwise become substantially too large.
  • Bias correction: Bandwidths scale asymptotically as h ∝ N^-1/6 for quadratic bias and h ∝ N^-1/10 for quartic bias, with multiplicative correction applied afterward.The empirical 1.1 factor gives significant gains over the basic Parzen estimator across most tested cases, although it slightly underestimates Gaussian bandwidths.

IV. CORRELATION LENGTHS AND SAMPLING ERROR ON PARAMETER MEANS

GetDist estimates Monte Carlo mean uncertainty and effective sample size while accounting for both sample weights and correlations. Its correlation-length calculations use autocorrelation structure and practical truncation to control sampling noise.

  • Sampling error: Weighted correlated samples require mean-error estimates that account jointly for importance weights and serial correlations.Weighted samples arise in importance and related sampling methods, while MCMC chains are typically highly correlated.
  • Effective sample size: Effective sample size is defined by comparing the variance of correlated weighted means with the variance from independent unit-weight samples.This extends the independent-sample relation σ²_X/N to an effective N_eff,X for correlated weighted samples.
  • Autocorrelation: The estimator uses weighted deviations whose equilibrium covariance depends on lag through an autocorrelation function.For d_i ≡ w_i(X_i − X̄), the expected pairwise product is represented by C_d(|i−j|).
  • Correlation lengths: Correlation lengths can be expressed in weight units for importance-sampled chains or sample units to indicate effective independence.The weight-unit correlation length scales with the arbitrary normalization of importance weights.
  • Practical estimation: In practice, the autocorrelation sum is truncated when correlation falls below a threshold such as 0.05 to limit sampling noise.The long-chain approximation assumes the relevant correlation length is much shorter than the chain length, although a finite estimator is also available.

V. CREDIBLE INTERVALS AND CONTOURS

GetDist derives credible intervals and two-dimensional regions from KDEs while explicitly accounting for active prior boundaries and distribution shape. It selects among equal-tailed, highest-density, and one-tailed reporting modes according to density behavior.

  • Interval construction: GetDist calculates credible intervals at multiple confidence levels, defaulting to 68%, 95%, and 99%.The multi-step procedure is designed to handle varied distribution shapes and parameter boundaries.
  • Range and priors: The KDE range is initialized from weighted quantiles, adjusted for nearby hard prior boundaries, and extended slightly where no active boundary is present.The default range_confidence of 0.001 retains 99.8% of the probability before boundary adjustments.
  • Boundary diagnosis: Boundary flags are set when density at an active prior edge exceeds the significant-density threshold, identifying posterior truncation by that prior.If both ends are prior-dominated, no interval limits are reported because the parameter is constrained by the priors rather than distinct tails.
  • Caveat: For multimodal distributions, the reported outermost limits can include lower-density regions between modes.The highest-density threshold remains useful for identifying enclosed regions, but the outermost interval need not exclude valleys between modes.
  • Tail limits and contours: When one boundary is active, the method reports a one-tailed limit; when neither is active, it reports a two-tailed interval using KDE and sample-based checks.Two-dimensional credible regions are obtained from KDE contours enclosing the requested probability, and may be disconnected for complex distributions.
  • Interval selection: For symmetric distributions, GetDist prefers sample-weight equal-tailed intervals, while asymmetric cases use highest-density intervals with equal endpoint density.This choice provides numerical stability for symmetric cases and avoids including low-probability regions that equal-tailed intervals can include for asymmetric distributions.

VI. DISCUSSION

The discussion identifies scope boundaries for GetDist’s KDE methods, including stationarity assumptions, axis-aligned prior treatment, complex multimodality, and task-specific bandwidth objectives. The package supports several sample formats and integrations, while remaining an open-source project.

  • Sampling assumptions: The current implementation assumes stationarity, which may not hold for all sampling methods such as nested sampling.Future work is identified for optimizing kernel estimators for non-stationary sampling distributions.
  • Prior boundaries: Bandwidth selection does not explicitly account for general prior boundaries beyond axis-aligned cases, potentially producing suboptimal choices for more general priors.The limitation concerns smoothing-scale optimization rather than the broader availability of boundary corrections.
  • Complex distributions: A global smoothing kernel may be far from optimal for highly multimodal likelihoods whose modes have different characteristic scales.The discussion suggests local adaptation, such as clustering methods, as a possible improvement.
  • Task scope: Kernel widths are optimized for density estimation, while other tasks such as tail confidence limits may require different optimization criteria.This creates a scope boundary for applying the current bandwidth objective beyond density estimation.
  • Interoperability: GetDist accepts NumPy arrays, COBAYA samples, ArviZ imports, and general text-based sample files, with tighter sampler integration left as a possible future improvement.The code base is maintained openly on GitHub and welcomes community contributions.
Loading 1910.13970v2…