PaperlayerView exact revision on arXiv

Source-linked AI summary

Efficient implementation of Markov chain Monte Carlo when using an unbiased likelihood estimator

Arnaud Doucet, Michael Pitt, George Deligiannidis, Robert Kohn

arXiv:1210.1871v4stat.ME

TL;DR

Unbiased likelihood estimators make Metropolis–Hastings feasible when the likelihood is intractable, but increasing Monte Carlo samples trades higher computational cost for lower asymptotic variance. The paper analyzes this trade-off for general proposal and target densities, derives Gaussian-noise guidelines, and shows that the recommended log-likelihood-noise standard deviation depends on baseline chain efficiency.

  • Problem

    Choosing the number of Monte Carlo samples requires balancing lower asymptotic variance against the increasing cost of constructing an unbiased likelihood estimator.

  • Method

    The paper bounds relative computing time for general proposal and target densities, relaxing the Gaussian noise assumption and deriving quantitative guidance under Gaussian noise.

  • Results

    The recommended log-likelihood-noise standard deviation is around 1.0 for efficient exact-likelihood chains, 1.7 for inefficient chains, and 1.2 when efficiency is unknown.

  • Takeaways & Limitations

    Selecting N through the log-likelihood-noise standard deviation provides practical guidance when exact-likelihood chain efficiency is unknown or the exact likelihood cannot be implemented.

  • Takeaways & Limitations

    The quantitative results rely on a Gaussian noise assumption, while the parameter-independent noise assumption is acknowledged not to hold in practical scenarios.

Abstract

from arXiv · show

When an unbiased estimator of the likelihood is used within a Metropolis--Hastings chain, it is necessary to trade off the number of Monte Carlo samples used to construct this estimator against the asymptotic variances of averages computed under this chain. Many Monte Carlo samples will typically result in Metropolis--Hastings averages with lower asymptotic variances than the corresponding Metropolis--Hastings averages using fewer samples. However, the computing time required to construct the likelihood estimator increases with the number of Monte Carlo samples. Under the assumption that the distribution of the additive noise introduced by the log-likelihood estimator is Gaussian with variance inversely proportional to the number of Monte Carlo samples and independent of the parameter value at which it is evaluated, we provide guidelines on the number of samples to select. We demonstrate our results by considering a stochastic volatility model applied to stock index returns.

1 Introduction

The pseudo-marginal algorithm uses unbiased likelihood estimators within Metropolis–Hastings, but choosing the number of Monte Carlo samples requires balancing estimator cost against asymptotic variance. This paper develops sample-size guidelines under broader proposal and target settings.

  • Background: The pseudo-marginal algorithm uses unbiased likelihood estimators within Metropolis–Hastings, commonly constructed with importance sampling or particle filters.Its theoretical properties have been established, and the method has attracted substantial interest in Bayesian statistics.
  • Core problem: The main practical decision is choosing N, the number of Monte Carlo samples or particles used to estimate the likelihood.For fixed N, the pseudo-marginal transition kernel leaves the posterior distribution invariant.
  • Core problem: More Monte Carlo samples usually lower the asymptotic variance of pseudo-marginal averages, but constructing the likelihood estimator becomes more expensive as N increases.Computing time is typically proportional to N times the average’s asymptotic variance.
  • Contribution: The paper minimizes explicit upper bounds on computing time for general proposal and target densities while relaxing the Gaussian noise assumption.Quantitative guidelines are then derived under a Gaussian noise assumption.
  • Contribution: The recommended standard deviation of log-likelihood noise is around 1.0 for efficient exact-likelihood Metropolis–Hastings, 1.7 for inefficient chains, and 1.2 when efficiency is unknown.These values provide practical guidance for selecting N when the exact-likelihood algorithm cannot be implemented.

2 Metropolis–Hastings method using an estimated likelihood

The estimated-likelihood Metropolis–Hastings scheme augments the parameter with likelihood-estimation noise and targets a joint density whose marginal is the desired posterior. Under a parameter-independent noise assumption, its transition kernel and acceptance probability can be analyzed explicitly.

  • Exact-likelihood baseline: The exact-likelihood Metropolis–Hastings chain uses proposal density q and acceptance probability based on the posterior and proposal-density ratio.The chain cannot be simulated when the likelihood is intractable.
  • Estimated likelihood: When the likelihood is intractable, a non-negative unbiased likelihood estimator based on auxiliary variables enables a Metropolis–Hastings construction.Unbiasedness ensures that the resulting joint density has the correct posterior marginal.
  • Pseudo-marginal construction: The pseudo-marginal algorithm targets the augmented density with proposal density q(θ, ·)m(·), retaining the desired posterior as its marginal.In practice, the recorded state uses θ and the log estimated likelihood rather than the auxiliary variables themselves.
  • Noise parameterization: The additive noise is defined as Z = log estimated likelihood − log likelihood and parameterizes the augmented target through its conditional density g(z | θ).This reparameterization replaces the auxiliary-variable representation with log-likelihood-estimation noise.
  • Transition kernel: With the noise parameterization, the pseudo-marginal transition kernel proposes a new parameter and noise value, then accepts using their noise difference and the exact-likelihood proposal ratio.The acceptance probability is αQ{(θ,z),(ϑ,w)} = min{1, exp(w − z) rex(θ,ϑ)}.
  • Noise assumption: Assumption 1 makes the noise density independent of θ, enabling detailed performance analysis of the pseudo-marginal algorithm.The authors note that this simplifying assumption is not satisfied in practical scenarios, although it may be reasonable where the noise varies little in high-probability regions.

3 Main results

The paper bounds the inefficiency and relative computing time of pseudo-marginal chains using a bounding kernel whose acceptance probability factorizes into exact-chain and noise components. These bounds yield computable guidance for choosing the log-likelihood noise scale while clarifying when the bounds approximate the actual chain.

  • 3.3 Inefficiency of the bounding chain: The bounding kernel Q* separates parameter and noise effects through a product-form acceptance probability, enabling an exact expression for its inefficiency.Peskun ordering gives IF(h, Q) ≤ IF(h, Q*), so Q* provides an upper bound for the pseudo-marginal chain.
  • 3.2 Inefficiency of Metropolis–Hastings type chains: The paper introduces a novel spectral representation of integrated autocorrelation time, which underpins the main inefficiency bounds.The representation is derived using the jump-chain formulation and spectral properties of an auxiliary kernel.
  • 3.4 Bounds on the relative inefficiency of the pseudo-marginal chain: Relative inefficiency RIF(h, Π) compares a kernel’s inefficiency with that of the exact-likelihood kernel Qex, and the paper bounds RIF(h, Q) using only exact-chain inefficiency in tractable corollaries.The resulting bounds avoid dependence on the full autocorrelation sequence appearing in the more intricate expression for RIF(h, Q*).
  • 3.5 Optimizing the computing time under a Gaussian assumption: Relative computing time is defined as RCT(h, Q; σ) = RIF(h, Q; σ) / σ^2, reflecting the inverse relationship between σ^2 and likelihood-estimator sample size.Because RIF is intractable, the paper minimizes upper bounds uRCTi(h; σ) instead.
  • 3.5 Optimizing the computing time under a Gaussian assumption: As exact-chain inefficiency increases, the computing-time bounds flatten and their minimizing σopt increases, making performance less sensitive to σ.For fixed σ, the bounds decrease with corresponding exact-chain inefficiency, suggesting a smaller pseudo-marginal penalty for less efficient exact chains.
  • 3.5 Optimizing the computing time under a Gaussian assumption: The two extreme scenarios give σopt = 0.92 for a perfect proposal and σopt = 1.68 for a very inefficient proposal.Using 1.68 instead of 0.92 raises relative computing time from 5.36 to 12.73 in the perfect-proposal scenario, whereas the reverse choice raises it from 1.51 to 2.29 in the inefficient scenario.
  • 3.5 Optimizing the computing time under a Gaussian assumption: For small-step random-walk proposals, large exact-chain inefficiency suggests RCT(h, Q*; σ) ≈ RCT(h, Q; σ).The approximation follows because the exact-chain acceptance probability tends toward unity as the step size decreases.
  • 3.5 Optimizing the computing time under a Gaussian assumption: The numerical results rely on the Gaussian noise assumption, although the relative-inefficiency bounds can also be calculated for other noise distributions subject to their conditions.The paper cautions that lower bounds apply to Q* but not generally to Q, while some upper bounds become exact for Q* as exact-chain inefficiency grows.

4 Application

The application evaluates the log-likelihood error assumptions and computing-time guidance in a two-factor stochastic-volatility model for S&P 500 returns. For T = 300 and T = 2700, Gaussian approximations are accurate, and the empirical computing-time optimum occurs at σ between 1.2 and 1.5.

  • Model and implementation: The model uses T daily S&P 500 returns and Bayesian inference for a 9-dimensional parameter vector with a bootstrap particle-filter likelihood estimator.The pseudo-marginal algorithm uses a multivariate Student-t random-walk proposal and N particles.
  • Log-likelihood error: For T = 40, the true log-likelihood-error histogram disagrees with the assumed Gaussian density, with lighter right tails and heavier left tails.The resulting discrepancy is smaller for the stationary target density, and the acceptance rate is higher than the Gaussian assumption suggests.
  • Log-likelihood error: For T = 300 and T = 2700, the assumed Gaussian densities closely match the corresponding histograms, while the central-limit approximation becomes very good.The posterior and stationary proposal distributions are sufficiently concentrated that the variance of the error changes little with θ.
  • Log-likelihood error: Figure 2 compares error histograms and theoretical densities across T = 40, 300, and 2700, using N = 4, 80, and 700, respectively.The columns show distributions associated with gN(z | θ) and πN(z | θ), evaluated at the posterior mean.
  • Pseudo-marginal performance: For most parameters, the optimal σ is between 1.2 and 1.5, corresponding to N = 40 and 60.The empirical results agree with the upper-bound analysis, and the two marginal-acceptance-probability curves in Figure 3 are similar.
  • Pseudo-marginal performance: Figure 3 reports inefficiencies and relative computing times against σ, averaging inefficiencies over the 9 parameter components, with marginal acceptance probability in the right panel.Relative computing time is formed from inefficiency and σ, using N = 2000 as the reference for the exact-likelihood inefficiency.

Appendix 1

Appendix 1 develops technical results that support the paper’s bounds on pseudo-marginal inefficiencies and relative computing times. The arguments use positivity, reversibility, spectral representations, and comparisons between exact, pseudo-marginal, and bounding kernels.

  • Kernel comparisons: The appendix proves inequalities linking inefficiencies under the exact-likelihood kernel, the pseudo-marginal kernel, and the bounding kernel Q∗.The comparison follows because Q∗ has lower acceptance probabilities, so Peskun ordering gives variance under Q∗ at least as large as under Q.
  • Spectral representation: The appendix derives a spectral representation for the inefficiency of the bounding kernel Q∗.This representation is used alongside the earlier inequalities to analyze the effect of likelihood-estimator noise.
  • Positivity and spectral arguments: Positivity of the auxiliary jump kernel implies that its inefficiency is at least 1, enabling several subsequent bounds.The appendix also uses reversibility and ergodicity to justify the relevant spectral arguments.
  • Bounds: The proofs establish upper bounds for relative inefficiency and relative computing time using spectral quantities and moments of the noise distribution.These bounds involve terms such as IF(h/̺ ex, eQex), πz(1/̺ z), and πz(̺ z).

Appendix 2

Appendix 2 uses lower and upper computing-time bounds to bracket the optimal log-likelihood-noise standard deviation. The resulting interval narrows in relative computing time as inefficiency increases, while the lower endpoint for σ increases.

  • Bounding the optimum: The appendix computes an interval containing σopt by comparing the lower bound lRCT1 with the minima of the upper bounds uRCT3 and uRCT4.Table 1 reports this interval and the corresponding lower- and upper-bound minima for different values of IF(h/̺ ex, eQex).
  • Bounding the optimum: As IF(h/̺ ex, eQex) increases, both endpoints of the interval containing optimal relative computing time decrease, while the lower endpoint of the σ interval increases.The appendix also supplies technical proofs, an additional example, and implementation information.

B Proof of Proposition 2

The proof of Proposition 2 transfers asymptotic-variance control between an ergodic Markov chain and its jump-chain representation. It establishes the required function-space, operator, and nonstationary-start results through Hilbert-space and martingale arguments.

  • Operator properties: Ergodicity and reversibility make the relevant operators positive, self-adjoint, injective, and equipped with densely defined inverse square roots.These properties support the spectral and operator arguments in the proof.
  • Functional-analytic setup: The proof constructs Hilbert spaces whose elements are functions with finite inefficiencies under the original and jump-chain kernels.These spaces are shown to be Hilbert spaces, and multiplication by 1/̺ is extended as an operator between them.
  • Operator extension: The proof uses density and closed-operator arguments to extend the transformation from a dense subspace and identify its limit with h/̺.The convergence argument uses completeness of the target Hilbert space and absolute continuity of the jump-chain invariant measure.
  • Functional-analytic setup: When IF(h, P) is finite, the transformed function h/̺ belongs to the corresponding jump-chain space and has finite inefficiency under the jump kernel.This is the key operator result used in the proposition.
  • Asymptotic argument: A martingale approximation and asymptotic negligibility argument establish the central-limit behavior for the jump-chain representation.The result remains valid when the chain starts from the invariant jump-chain distribution, and the asymptotic variance is unchanged.
  • Conclusion: The proposition concludes that the transformed jump-chain functional has finite asymptotic variance under the stated ergodicity and integrability conditions.The proof combines the martingale limit, negligibility, and the finite-inefficiency result for h/̺.

C Proofs of other technical results in the main paper

The proofs establish irreducibility, ergodicity, small-set structure, and geometric drift properties for the relevant pseudo-marginal kernels. These properties support the stated finiteness and ergodicity results.

  • C = (−∞, z0] is established as a small set using a probability measure concentrated on [0, z0] ⊂ C.
  • The drift function V(z) = 1/̺z(z) yields geometric ergodicity of eQz, implying the required finiteness inequality.
  • The n-step density of Q* is strictly positive on the relevant product set, establishing its irreducibility and ergodicity.
  • The extended exact-likelihood kernel eQex is ergodic by arguments paralleling those for Qex.

D Statements and proofs of auxiliary technical results

The auxiliary results characterize how relative computing time depends on log-likelihood noise and exact-chain inefficiency. They identify optimal noise levels and show that the optimum shifts as exact-chain inefficiency increases.

  • When IF(h, Qex) = 1, relative computing time is minimized at σopt = 0.92.The associated relative inefficiency value is reported in the proposition but is truncated in the supplied passage.
  • As IF(h, Qex) tends to infinity, σopt increases to 1.02.
  • For IF(h, Qex) ≥ 1, the second derivative of relative computing time with respect to σ is positive.
  • The general optimum σopt = 1.02 is verified numerically by setting the derivative of relative computing time to zero.
  • Relative inefficiency and relative computing time decrease as IF(h, Qex) increases, holding σ constant.

E Asymptotic upper bound

Under Gaussian noise, the relative inefficiency and computing-time bounds have similar shapes across noise levels and attain minima at nearby noise standard deviations.

  • The exact relative inefficiency approaches the lower bound as l → 0 and approaches Ψ(σ) as l → ∞.
  • The relative computing-time curves are very similar in shape as functions of σ, regardless of l.
  • The lower-bound computing time is minimized at σ1 = 1.68, while the limiting inefficiency expression is minimized at σ2 = 2.00.
  • At σ1 = 1.68 and σ2 = 2.00, the corresponding values are lRCT2(σ1) = 1.51 and lRCT2(σ2) = 1.59.

F Simulation results

The simulations apply the pseudo-marginal algorithm to a linear Gaussian state-space model and add results for the stochastic volatility model. The state-space likelihood is also evaluated exactly through the Kalman filter.

  • The linear Gaussian state-space model is a first-order autoregression AR(1) observed with noise.
  • The parameters of interest are θ = (φ, μx, σx).
  • The Kalman filter calculates the likelihood exactly, so the true likelihood term is known when evaluating log-likelihood error.

F.1 Empirical results for the error of the log-likelihood estimator

The experiments assess whether the log-likelihood error follows the assumed Gaussian form and whether its distribution changes substantially across posterior parameter values. Agreement improves for longer time series, while the shortest series shows heavier-tailed or discrepant behavior.

  • For T = 40, the conditional error histogram differs slightly from the assumed Gaussian density and is heavier-tailed in the left tail.For T = 300 and T = 2700, the assumed Gaussian densities are very accurate.
  • For T = 300 and T = 2700, posterior-marginal error histograms closely match the assumed densities, supporting Assumptions 1 and 2.The agreement is reported for the error marginalized over θ distributed according to π(θ).
  • For T = 300, conditional moments vary only slightly across 100 posterior parameter values around expected values −0.5, 1, 0 and 3.The experiment estimates means, variances, and third and fourth central moments using 300 independent particle-filter runs per parameter value.
  • For the stochastic-volatility experiment, variability is greater because the true likelihood must be estimated, but it decreases as T rises to 2700.The resulting behavior is described as similar to the AR(1) plus noise experiment.

F.2 Empirical results for the pseudo-marginal algorithm

The pseudo-marginal experiments evaluate computing time as the likelihood-estimator noise varies under proposals with different efficiencies. Efficient proposals favor noise standard deviation near unity, whereas inefficient proposals favor larger noise and produce flatter computing-time curves.

  • The AR(1) plus noise study uses T = 300 data, a Kalman-filter likelihood for the exact chain, and particle-filter estimates for the pseudo-marginal chain.The proposal uses transformed parameters and an autoregressive Student-t random walk with ρ values from 0 to 0.9.
  • For N = 60 with σ(θ; N) = 0.92, the approximating asymptotic Gaussian densities closely match histograms of accepted and rejected log-likelihood errors.The histograms use marginal error values over posterior draws generated by the pseudo-marginal scheme.
  • For the independent proposal, computing time is minimized around N = 43 or 60, with corresponding σ values 1.11 and 0.92.These results support an optimal σ close to unity for an efficient proposal.
  • For the more persistent proposal with ρ = 0.9, the optimal N is around 31 for all three parameters, with σ(θ; N) = 1.34.The findings are consistent with a larger optimal noise level for an inefficient proposal.
  • As proposal inefficiency increases, the optimal σ increases and relative computing time becomes flatter as a function of σ.The bounds likewise become flatter, suggesting reduced sensitivity to the precise choice of σ.
  • Theoretical bounds are evaluated through numerical quadrature, autocorrelation summation, and Monte Carlo integration of the relevant terms.Autocorrelations descend rapidly, making the inefficiency estimates straightforward to compute.
Loading 1210.1871v4…