The Impact of Loss Estimation on Gibbs Measures
David T. Frazier Jeremias Knoblauch Christopher Drovandi



	#1 1
	
	
	
	The Impact of Loss Estimation on Gibbs Measures
    David T. Frazier^*, Jeremias Knoblauch^*, and Christopher Drovandi
    
======================================================================

1.8 
	
    
    ^*Both authors contributed equally.
    In recent years, the shortcomings of Bayes posteriors as inferential devices has received increased attention.
    
    A popular strategy for fixing them has been to instead target a Gibbs measure based on losses that connect a parameter of interest to observed data.
    
    While existing theory for such inference procedures relies on  these losses to be analytically available, in many situations these losses must be stochastically estimated using pseudo-observations.
    
    The current paper fills this research gap, and derives the first asymptotic theory for Gibbs measures based on estimated losses.
    
    Our findings reveal that the number of pseudo-observations required to accurately approximate the exact Gibbs measure depends on the rates at which the bias and variance of the estimated loss converge to zero.   
    
    These results are particularly consequential for the emerging field of generalised Bayesian inference, for  estimated intractable likelihoods, and for biased pseudo-marginal approaches. 
    
    We apply our results to three Gibbs measures that have been proposed to deal with intractable likelihoods and model misspecification by <cit.>, <cit.>, and <cit.>.
    
	






PART: 
 



§ INTRODUCTION



Bayesian 
inference has long been the gold standard for principled statistical methodology reliant on the inclusion of prior knowledge.

Beyond that, the Bayes posterior provides a natural way for quantifying uncertainty about the data-generating process under study <cit.>.

However, orthodox Bayesianism also suffers from various pathologies in modern data-rich environments, complex models, and machine learning tasks <cit.>.



For example, under model misspecification, the Bayes posterior concentrates towards sub-optimal parameter values and  provides miscalibrated parameter inference <cit.>.

Worse still, it produces brittle inferences in the presence of outliers and heterogeneous observations <cit.>.

An entirely different disadvantage emerges  in models with intractable likelihood functions <cit.>.

For example, since traditional methods for Bayesian computation rely on  analytically available likelihood functions, they cannot be deployed for simulation-based models.

To circumnavigate this, methods like approximate Bayesian computation (ABC) obtain posterior distributions by approximating the likelihood function itself <cit.>.

A similar problem arises for likelihoods that are only known up to an intractable normalisation constant: once again, standard Markov chain Monte Carlo (MCMC) methods cannot be deployed without first approximating the likelihood function itself <cit.>.


While model misspecification and intractable likelihoods are unrelated problems, various solutions proposed to both problems take the same form: directly changing the inference target.

More specifically, rather than focusing on the Bayes posterior, a host of contemporary proposals instead base parameter inferences on 

alternative Gibbs measures <cit.>.

While algorithmically straightforward, this  leads to a new and hitherto unaddressed theoretical challenge. 

Specifically, many of the desired Gibbs measures must be constructed from losses that cannot be computed without stochastic estimation based on simulations from a statistical model for the data.

Notable examples include  β-divergence losses (see <Ref>) and the maximum mean discrepancy (see <Ref>), both of which are  popular in the literature on generalised Bayes <cit.>.

Much like in ABC, these losses lead to a discrepancy between the nominal target measure and the de-facto target from which an MCMC algorithm  produces samples.


This poses several practical problems: firstly, it is completely unclear what level of computational effort is required to make the de-facto target sufficiently close to the nominally targeted Gibbs measure.

In fact, our lack of understanding about the approximation quality raises the question under which conditions this new class of algorithms and Gibbs measures constitutes a preferable alternative to Bayes posteriors.

The existing literature provides little beyond heuristic arguments.

The current paper addresses this research gap, and provides the first suite of theoretical guarantees for algorithms that sample from Gibbs measures based on estimated losses.



§ MOTIVATION AND SETUP


Throughout, we consider the finite-dimensional data set =(y_1,…,y_n)^⊤ with y_i∈𝒴 for all i=1,…,n, and where we assume y_1:n∼ P_0.

While the true data-generating measure P_0 will generally depend on n, we suppress this for notational convenience.

Using our data, we model P_0 through the class of parametric models 

rCl
	{P_θ: θ∈Θ} ⊂𝒫(𝒴^n) for Θ⊆ℝ^d,
    

where we write p_θ to denote the  density of P_θ, and  once again suppress dependence of P_θ and p_θ  on n.

In a Bayesian approach to statistical modelling, we express our prior beliefs about the true state of the world probabilistically through the prior measure Π∈𝒫(Θ) with density π. 

Given p_θ and π, the Bayesian approach for quantifying uncertainty about the value of θ that best describes the observed data proceeds via Bayes' Rule, which yields the Bayes posterior measure with density

rCl
    π(θ|)=p_θ()π(θ)/∫_Θp_θ()π(θ)θ.
    

Conceptually, the Bayes posterior updates our prior beliefs about the world with data, and represents an internally coherent belief update.

Crucially,  our beliefs about the world are expressed exclusively in terms of θ—a direct consequence of the fact that the  Bayesian approach forces us to explicitly posit a model for how the data  was generated.




 §.§ New Posterior Belief Distributions


It is easy to generalise the Bayes posterior: letting ω>0 be a scalar and _n:Θ×𝒴^n →ℝ a loss function expressing the information about the parameter θ contained in , then so long as ∫_Θexp{-ω·_n(θ, y_1:n)}π(θ)θ<∞, we can define the Gibbs measure

rCl
    π(θ|_n) := exp{-ω·_n(θ, y_1:n)}π(θ)/∫_Θexp{-ω·_n(θ, y_1:n)}π(θ)θ,
    


which recovers the Bayes posterior π(θ|) when _n(θ, ) = - log p_θ() and ω=1.

Gibbs measures make appearances in many areas of Bayesian methodology, including ABC  <cit.>, Provably Approximately Correct (PAC) Bayes <cit.>,  and generalised Bayesian methods <cit.>.



  §.§.§ Approximate Bayesian Computation (ABC)



The motivation of ABC algorithms is the approximation of Bayes posteriors for simulator models with intractable likelihoods, where p_θ is not available analytically, but 
can be sampled with low computational overhead using simulated data U_1:miid∼ P_θ

for each θ∈Θ. See <cit.> for a handbook of these methods, and <cit.> for a recent review. 

Given a vector of summary statistics S(), a distance d, and a positive tolerance ε, ABC methods can be interpreted as sampling algorithms targeting a Gibbs measure given by (<ref>) with _n^ε(θ, ) = - log𝔼_U_1:miid∼ P_θ[ 1_{d{S(), S(U_1:m)}≤ε}· p_θ(U_1:m)] and ω=1.

This measure is used as an approximation to the Bayes posterior π(θ|) with its accuracy depending on the choice of S,d, and ε.

Exemplifying the more general problem under study in the current paper, the expectation in _n^ε is unavailable.

ABC approximates this expectation by drawing pseudo-observations u_1:miid∼ P_θ to construct an estimator for _n^ε.


An alternative to ABC in the setting of simulation-based inference is Bayesian synthetic likelihood <cit.>. 

The corresponding target posterior is a Gibbs measure as in (<ref>) with loss 

rCl
    _n(θ,)=-1/2lndet{Σ(θ)}-1/2{S_n()-b(θ)}^⊤Σ(θ)^-1{S_n()-b(θ)}.
    

Here, 
Σ(θ)=_U_1:miid∼ P_θ{S(U_1:m)} and b(θ)=𝔼_U_1:miid∼ P_θ{S(U_1:m)}. 

As with ABC, the expectations defining _n are generally intractable, so that Bayesian synthetic likelihood approximates them with m samples from P_θ. 




  §.§.§ Provably Approximately Correct (PAC) Bayes


Provably approximately correct (PAC) Bayesian results have direct links with Bayesian inference, but have distinct origins and motivations <cit.>.

In particular, despite their name and the re-appearance of a `prior' measure Π, they do not rely on Bayes' Theorem.

Rather, they are justified through generalisation error bounds.

To illustrate this, we adapt Theorem 2.1 in <cit.> which itself is a version of the bound in <cit.>. 

Taking 𝒫(Θ) as the set of probability measures on Θ,
assuming that  are sampled independently from P_0, and supposing 0 ≤_n(θ, ) ≤ C for all ∈𝒴^n and θ∈Θ, it holds with probability at least 1-ε that for all Q ∈𝒫(Θ),

rCl
    ∫𝔼_Y_1:n iid∼P_0[ω·_n(θ, Y_1:n)] Q(dθ) 
        ≤    
    ∫ω·_n(θ, Y_1:n)Q(dθ)
    +
    KL(Q  Π)
    +
    ω^2 C^2/8n
    -
    logε,
    

where KL denotes the Kullback-Leibler divergence.

Minimizing the right hand side of this bound with respect to Q ∈𝒫(Θ), one  
recovers (<ref>).

A proof can be found in Theorem 1 of <cit.>, which follows by the arguments of <cit.>.

For other significant PAC-Bayes bounds that  recover Gibbs measures, see also <cit.>, <cit.>, or <cit.>.






  §.§.§ Generalised Bayes


To provide reliable and well-calibrated uncertainty quantification, the Bayes posterior requires that the model is well-specified.

If this condition is violated, π(θ|) lacks robustness, and produces miscalibrated representations of parameter uncertainty.

A recent line of work  has therefore argued for generalised Bayesian beliefs that are the outcome of an optimisation problem over 𝒫(Θ) <cit.>, recovering Gibbs measures as 

rCl
    π(θ|_n)
        =    
    _q ∈(Θ){∫_Θω·_n(θ, y_1:n)q()+KL(q||π)}.
      


When _n is additive over the data so that _n(θ, ) = ∑_i=1^n ℓ_i(θ, y_i) for some sequence of losses ℓ_i: Θ×𝒴→ℝ, such Gibbs measures coincide with general Bayesian updating  as advocated for by <cit.>.

The implication of (<ref>) is significant: as the Bayes posterior is the special case associated to _n(θ)=-log p_θ(), some of the drawbacks of relying on Bayes' Theorem can be addressed by choosing different losses _n.


In particular, noting the brittleness of the standard Bayes posterior π(θ| y_1:n) under model misspecification, 
various generalised posteriors are based on losses _n that guarantee robustness.

Most of these are built on discrepancy measures between P_θ and P_0 such as scoring rules <cit.>, α-, β-, and γ-divergences <cit.>, or the maximum mean discrepancy  <cit.>.

More recent work has also combined this desire for robustness with improved computational properties: for example, Stein discrepancies <cit.>
and generalisations of score matching
<cit.> provide conjugate posteriors for exponential family likelihoods with intractable normalisers.





 §.§ Computational Challenges



While many methodologically interesting Gibbs measures are straightforward to compute, several expose a computational hurdle whose impact on asymptotic behaviour has consistently been left unaddressed in existing contributions  <cit.>.

In particular, many methods in the literature rely on losses _n that are intractable, and  have to be estimated through a proxy _m,n that relies on m samples u_1:m with u_j ∼ P_θ for j=1,2,… m. 

Yet, prior theory operates directly on the idealised target posterior π(θ|_n)—rather than the approximate version based on the estimated loss _m,n—several examples of which are given in the following section.

Not accounting for the randomness in _m,n means that no existing results tell us how to choose m in practice, and that existing asymptotic results often do not apply to the Gibbs measure we compute in practice.


In the current paper, we address this issue head-on, and provide comprehensive results on asymptotic behaviours in both m and n.

The results are the first of their kind,  and essential in understanding exactly how large m should be to ensure that _m,n leads to reliable approximations of π(θ|_n).

While our findings are most immediately relevant for the emerging field of generalised Bayesian posteriors, they also inform a host of other settings.

For example, they provide the first set of generic asymptotic results for MCMC-based inference based on biased likelihood estimates, extending existing results on the behavior of pseudo-marginal methods to the case where the likelihood estimator is  biased.





 §.§ Applicability of Results



To provide a clear exposition to the  research gap our results are closing, we  provide three examples of posteriors that suffer from computationally intractable losses.

The first and second are generalised Bayes posteriors based on the β-divergence and the maximum mean discrepancy (MMD), respectively.

These two generalised posteriors are popular in the literature, and have not only been proposed for general methodology <cit.>, but have also  found specific use for filtering problems <cit.>,  changepoint detection <cit.>, differential privacy <cit.>,  Bayesian deep learning <cit.>, and simulation-based inference <cit.>.

While asymptotic results were derived by <cit.> and 
<cit.>, they are based on the exact loss _n, rather than the posterior calculated from the estimate _m,n that must be used in practice.

The third is a Bayes posterior whose intractable likelihood model is approximated using kernel density estimation (KDE) introduced by <cit.>, and which has similarities with more contemporary—albeit much harder to analyse—black box density estimation methods  <cit.>.

The use of KDE for intractable models can deliver more informative  inferences than classical ABC based on summary statistics <cit.>, but the behaviour of the resulting posteriors has not been studied previously. 





  §.§.§ β-divergence





    Suggested for parameter estimation by <cit.>, the β-divergence is also often called the density power divergence.
    
    As the second name suggests, it is a discrepancy that relies on the existence of densities, and  owes its name to the fact that these densities are raised to a power β>0.
    
    
    
    For a given β, this divergence is written as
    
    rCl
        _β(p_0,p_θ)=∫{p_θ^1+β(y)-(1+1/β) p_0(y) p_θ(y)^β+1/β p_0^1+β(y)} y.
        

    
    If we have access to a sample y_1:n of n data points for which y_i iid∼ p_0, and if we wish to minimise this expression over θ, we can ignore the last term and obtain the loss function
    
    rCl
        ^β_n(θ, y_1:n)= n∫p_θ^1+β(u)u
        -(1+1/β) ∑_i=1^n p_θ(y_i)^β.
        
    
    Minimising this loss yields an estimator for θ that converges to the true data-generating parameter as n→∞ if the model is well-specified, and which exhibits robustness under 
    under model misspecification <cit.>.
    
    Here, β determines the degree of robustness: the loss becomes more robust as β→∞, but less robust as β→ 0.
    
    As β→ 0, its limit eventually recovers the non-robust negative log likelihood -log p_θ(y_1:n).
    
    This makes _n^β ideally suited for the computation of robust generalised Bayes posteriors, as it  recovers the standard Bayes posterior for β→ 0.
    
    However, outside a small subset of exponential families, the integral ∫ p_θ^1+β(u) u needs to be approximated.
    
    As a result, one generally has to draw samples u_j iid∼ p_θ for j=1,2,… m, and then use the approximation
    
     rCl
        ^β_m,n(θ, y_1:n)=  
        n/m∑_j=1^mp_θ^β(u_j)
        -(1+1/β) ∑_i=1^n p_θ(y_i)^β.
        
    
    
    This means that the actual target posterior depends on ^β_m,n rather than ^β_n.
    
    This is not taken into account in existing theory, and prior work has opted to study the idealised target π(θ|^β_n)  instead <cit.>.




  §.§.§ Maximum Mean Discrepancy




The maximum mean discrepancy (MMD) is a popular distance for testing and parameter estimation <cit.>. 

There are generally two reasons for using the MMD as a loss in the context of generalised Bayes posteriors.

The first is robustness, which can be formally shown to hold whenever the kernel k:𝒴×𝒴→ℝ_+ it is based on is bounded <cit.>.

The second relates to computation: unlike most other losses, the MMD does not require an analytically tractable likelihood function. 

Instead, all we need is access to a generative model for P_θ from which we can draw samples, which makes the MMD particularly useful in simulation-based inference <cit.>.


For a positive-definite and symmetric kernel function k:𝒴×𝒴→ℝ_≥ 0 associated to the reproducing kernel Hilbert space ℋ with norm ·_ℋ, the corresponding MMD between the empirical measure P_n = 1/n∑_i=1^nδ_y_i and P_θ is given by

rCl
    MMD_k(P_n, P_θ)
        =    
    sup_f ∈ℋ: f_ℋ≤1|
        𝔼_U∼P_θ[ f(U) ]
        -
        1/n∑_i=1^n f(y_i)
    |.
    

This is an integral probability metric <cit.>, and appears intractable at first glance.

Fortunately, one can solve the underlying optimisation problem analytically. 

After squaring, this becomes 

rCl
    MMD_k(P_n, P_θ)^2
        =    
    1/n^2∑_i=1^n∑_i' = 1^n
        k(y_i, y_i')
    -
    21/n∑_i=1^n𝔼_U∼P_θ[ k(U, y_i) ]
    +
    𝔼_U ∼P_θ, U' ∼P_θ[ 
        k(U, U')
    ].
    

As the first term does not depend on θ, a natural loss to construct from the above is  

rCl
    ^k_n(θ, y_1:n)
        =    
    -
    21/n∑_i=1^n𝔼_U∼P_θ[ k(U, y_i) ]
    +
    𝔼_U ∼P_θ, U' ∼P_θ[ 
        k(U, U')
    ],
    

which is also often referred to as the kernel score <cit.>,
and which relies on a number of intractable expectations.  

Just as for the expectations featuring in the β-divergence, one can approximate them by taking u_j iid∼ P_θ for j=1,2,… m, which results in the estimator

rCl
    ^k_m,n(θ, y_1:n)
        =    
    -
    21/nm∑_i=1^n∑_j=1^m  k(u_j, y_i) 
    +
    1/m^2∑_j=1^m∑_j'=1^m
        k(u_j, u_j').
    


Similarly to the case of generalised Bayes posteriors based on β-divergences, MMD-based posteriors are computed using ^k_m,n, while existing asymptotic theory is applicable only to the idealised version π(θ|^k_n) <cit.>.




  §.§.§ Approximating Intractable Likelihoods via Kernel Density Estimation




Our third and last example relates to a Bayes posterior whose likelihood function p_θ is not analytically available, but where for each θ, it is easy to generate samples from the corresponding model  P_θ.

In such situations, working with the Bayes posterior π(θ|) demands estimation of the normalised likelihood.

One way of doing so is via kernel density estimation (KDE) with a suitable isotropic kernel function k:𝒴×𝒴→ℝ_+, and bandwidth h≥0. 

Given simulated data u_jiid∼ P_θ for j=1,2,… m, we can estimate the intractable likelihood using the density estimate:






rCl
    p̂_θ(y)
        =     1/mh∑_j=1^m k_h(y - u_j) ,  k_h(x)=k(x/h),
    

where we have followed conventional practice and overloaded notation as k(y,y')=k(y-y'), for y,y'∈𝒴.

From the above, we estimate the negative log likelihood _n^h(θ, ) = -∑_i=1^nlog p_θ(y_i) as

rCl
    ^h_m,n(θ, y_1:n)
        =    
    -∑_i=1^n logp̂_θ(y_i).
    

<cit.> and <cit.> use  ^h_m,n(θ, y_1:n) within MCMC to conduct inference on θ even thought the actual density is intractable.[The above estimator differs slightly from the one suggested by <cit.>, which is based on simulating m=M· n data points for some M∈ℕ. The two approaches are equivalent when the simulated data are iid. ] 


More recently, alternative density estimation methods have been proposed that replace the KDE-based estimator with approaches that instead estimate the likelihood using neural networks <cit.>. 

While the reliance on black box models makes it impossible to explicitly verify our  theory in this setting, the results we derive for the KDE-based estimator should be understood as the best case scenario we could hope for with these sequential neural likelihood approaches.




  §.§.§ Wider Applicability



The three estimators presented above are constructed using three distinct principles (averages, U-statistics, and kernel density estimators), and our theory applies to all  of them.

While we will confine ourselves to these three examples throughout the remainder of the paper, it is important to note our setup applies to a much larger class of posteriors which rely on approximating some idealised loss _n with an estimator _m,n.

This includes posteriors based on Hellinger distances <cit.> and  γ-divergences <cit.>, sliced distances <cit.>, as well as the so-called ρ-Bayes approach <cit.>.

In the context of ABC, it also applies to the  Wasserstein-ABC posterior <cit.>, 
approaches based on general integral probability metrics <cit.>, as well as other distance-based variants of ABC <cit.>.

For all of these cases, previous work has exclusively studied theoretical properties of the idealised posterior π(θ|_n), not to the computationally feasible version based on _m,n.

In addition, the setting also includes posteriors computed using MCMC based on biased likelihood estimates
<cit.>.





 §.§ Motivating Example



Before providing our theoretical results, we motivate their necessity using the losses described in the previous section and demonstrate the practical ramifications of using _m,n instead of _n in posterior computations.

To this end, we simulate n=100 observations from a Student's t location model with fixed unit variance and 10 degrees of freedom.

Inference is performed only on the location parameter, which is set to θ=1.

We then construct generalised Bayes posteriors based on the losses introduced in <Ref> for different m, and compare their accuracy.

The results are collected in <Ref>, and demonstrate that even large values of m can lead to poor approximations.



Beyond that, the results show that across all three  loss functions, the posteriors become more concentrated around the true value θ=1 as m increases.

For the MMD and KDE losses, we see a marked increase in estimation accuracy when using more draws, while this effect is less pronounced for the β-divergence loss. 

Put together, this simple example illustrates that for accurate approximations of the target Gibbs measure, we must carefully choose the number of draws m used in our loss estimates. 

In particular, in the MMD and KDE cases it is clear that location and variability of the posteriors change as m increases.

This suggests that the choice of m will itself also influence concentration of the posterior, an empirical finding that we will verify theoretically in the next section. 





§ A NEW ASYMPTOTIC THEORY



To set the stage, we first note the clear parallels with pseudo-marginal MCMC: exactly as in that setting, using _m,n instead of _n to construct our posterior implies that we do not actually target π(θ|_n).

Instead, when we construct a sampling algorithm based on _m,n, we obtain an unbiased estimator for the implicit target posterior 

rCl
    π(θ|_m,n)         ∝      π(θ) ·_U_1:m iid∼ P_θ[
        exp{-ω·_m,n(θ, y_1:n)}
    ].
    

Here, _U_1:miid∼ P_θ[ exp{-ω·_m,n(θ, y_1:n)}] is the expected 
exponentiated loss exp{-ω·_m,n(θ, y_1:n)}.

The expectation is taken over the distribution of  the samples u_1:m  used in constructing _m,n(θ, y_1:n).

In the case of the β-divergence and MMD  losses of Sections <ref> and <ref> for instance, these samples are used to approximate an intractable integral.

The expectation over U_1:miid∼ P_θ is crucial when defining the implicit target: it reflects the fact that using sampling methods like MCMC to approximately sample from Gibbs measures π(θ|_n) would recompute _m,n based on a new sample u_1:miid∼ P_θ at each iteration of the sampling algorithm—rather than re-using a fixed sample at each iteration.

This resulting target posterior π(θ|_m,n) is at best a biased estimator for π(θ|_n), even if the loss estimator is unbiased so that _U_1:miid∼ P_θ[  _m,n(θ, y_1:n) ] = _n(θ, y_1:n).

In particular, whenever _m,n(θ, y_1:n) is non-constant in u_1:m, Jensen's Inequality implies 

rCl
	_U_1:m iid∼ P_θ[
        exp{-ω_m,n(θ, y_1:n)}
    ]
        ≠   
    
    exp{-ω·_U_1:m iid∼ P_θ[ _m,n(θ, y_1:n) 
    ]}.
    

This inequality raises two questions:

    
  1. For a fixed n, how different are π_(|_m,n) and π(|_n) as a function of m?
    
  2. As n→∞, how fast does m need to increase for  π_(|_m,n) to behave like π(|_n)?


In the remainder, we address this by deriving the rate at which π(θ|_m,n) converges to π(θ|_n) (<Ref>), deriving the posterior concentration of π(θ|_m,n) as a function of both n and m (<Ref>, <Ref>), 
and by establishing the size of m required to ensure that a Bernstein-von-Mises type result for π(θ|_n) transfers  to  π(θ|_m,n)  (<Ref>). 




 §.§ Preliminaries


To state our assumptions and results, we make use of several common notational shorthands.

For two sequences of scalars {a_n}_n∈ℕ and {b_n}_n∈ℕ, we take a_n ≲ b_n to mean that for all n large enough and a constant C>0, we have a_n ≤ Cb_n. 

Moreover, we take a_n ≍ b_n to mean that both a_n ≲ b_n and b_n ≲ a_n.


Similarly, for given sequences of random variables {X_n}_n∈ℕ and scalars {a_n}_n∈ℕ, we take X_n ≲ a_n to mean that the sequence {X_n/a_n}_n∈ℕ is  stochastically bounded by a constant C >0 almost surely, so that ℙ( lim_n→∞ X_n/a_n ≤ C ) = 1. 

Similarly, we write X_n≍ a_n if a_n≲ X_n and X_n≲ a_n.

It will be clear from context whether we use this notation for asymptotics in  m only, or for joint asymptotics in n with m=m(n).

In addition, we also use  X_n = O_p(a_n) to mean that there is a C>0 so that  lim_n→∞( X_n / a_n ≤ C ) = 1, and X_n = o_p(a_n) to mean that for all ε>0,  we have lim_n→∞ℙ(  X_n/a_n < ε) = 1.




Further, for Θ⊆ℝ^d, and a function ℓ:Θ→ℝ, we denote ∇_θℓ(θ_0) as the vector of first partial derivatives of ℓ evaluated at θ_0, and ∇_θ^2 ℓ(θ_0) as the d× d matrix of its second partial derivatives.




We will also use x to denote the usual Euclidean norm for any vector x ∈ℝ^d, and min{z, y } (max{z, y }) for the minimum (maximum) between two scalars z,y ∈ℝ.





 §.§ Approximation Quality (fixed n, large m)


To ensure that π(θ|_m,n) is a good approximation for π(θ|_n) as m→∞, a
minimal condition to impose is that the loss estimate _m,n is pointwise consistent for _n.

While it is easy to demonstrate that   |π(θ|_m,n) - π(θ|_n)| vanishes for all θ∈Θ as m→∞ under this assumption, such a result would have little practical use for two main reasons:

it tells us nothing about the speed of convergence, and it does not necessarily even imply  convergence in distribution of π(θ|_m,n) to π(θ|_n).



We address this by deriving the rate at which π(θ|_m,n) converges to π(θ|_n) in  total variation as n is held fixed and m→∞ (cf. <Ref>).

To achieve this, we require some regularity conditions for the estimation quality of _m,n  stated  below.





For constants κ_1, κ_2 > 0 and a sequence of functions {σ_n^2(θ)}_n∈ℕ with σ_n(θ): Θ→ℝ_+  allowed to depend on y_1:n, it holds for any fixed  that 


    (i)  _U_1:miid∼ P_θ[_m,n(θ, y_1:n)] - _n(θ, ) ≲σ^2_n(θ)m^-κ_1 as m→∞; 
    
    
    (ii) _U_1:miid∼ P_θ[ {_m,n(θ, y_1:n)-_U_1:miid∼ P_θ[_m,n(θ, y_1:n)]}^2] ≲σ^2_n(θ)  m^-κ_2.
    






The interpretation of <Ref>  is  straightforward:  <ref> allows _m,n to be a biased estimator for _n, but requires that the pointwise bias goes to zero like m^-κ_1 as  m →∞ for some κ_1.

Similarly, <ref> ensures that the estimation error of _m,n(θ, y_1:n) goes to zero as m^-κ_2 for some κ_2 and as  m →∞.

On a technical level, these requirements  allow the derivation of a pointwise asymptotic bound of the moment-generating function of the random variable -_m,n(θ, y_1:n) in U_1:m with the constant exp{-_n(θ, y_1:n)}, so that _U_1:miid∼ P_θ[ exp{-_m,n(θ, y_1:n)}] ≲exp{-_n(θ, y_1:n)}+σ^2_m,n(θ)m^-κ_2. 

This explains why <Ref> looks superficially similar to classic conditions  required for concentration inequalities for sub-Gaussian random variables  <cit.>.

Note however that <Ref> is much weaker: classic sub-Gaussianity requires the left-hand side of <ref> to be zero for all m, and that  <ref> holds for all finite m (rather than asymptotically as m→∞).

In this sense, <Ref> is akin to an asymptotic sub-Gaussianity assumption, with σ_n^2(θ) serving as a form of asymptotic variance proxy (see <Ref>).





Having upper bounded bias and variance in terms of σ^2_n(θ), it is  straightforward to define interpretable conditions that allow us to strenghten pointwise to posterior-averaged convergence.

In particular, if certain moments  of π(θ) and π(θ|_n) exist, then the posterior-averaged biases and variances will vanish uniformly over Θ at rates κ_1 and κ_2, respectively.



	For any n≥1 and some p ≥ 1, 
	∫_Θ^pπ_(|_n) < ∞, 
    ∫_Θθ^pσ^2_n(θ)π(θ)θ< ∞,
    ∫_Θσ^2_n(θ)θ^pπ(θ|_n)θ< ∞, 
    0<∫_Θexp{-ω·_n(θ)}π(θ)θ<∞,
    ∫_Θθ^pπ(θ)θ < ∞.


Under <Ref>, and  the moment conditions in  <Ref>, we can control the difference term |π(θ|_m,n) - π(θ|_n)| uniformly over Θ, for any n≥ 1, and as m→∞.




theoremmasymptotics

	If Assumptions <ref> and <ref> are satisfied, then for all ξ∈ [0,2], as  m→∞  
	
	rCl
		∫_Θθ^ξ|π(θ|_m,n)-π(θ|_n)|θ    ≲   m^-min{κ_1,κ_2}
		
	
	which also implies that
	
	rCl
		∫_Θ|π(θ|_m,n)-π(θ|_n)|θ    ≲   m^-min{κ_1,κ_2}, 

		∫_Θθπ(θ|_m,n)θ-∫_Θθπ(θ|_n)θ     ≲   m^-min{κ_1,κ_2}.
		



The above result holds conditionally for any sample . 

While Theorem <ref> has three parts, the last two equations are direct consequences of the first, and correspond to the special cases of convergence in total variation (ξ=0) and convergence in mean total variation (ξ=1).

Taken together,  <Ref> tells us that under mild regularity conditions, inferences based on π(θ|_m,n) can be made arbitrarily close to those based on the infeasible target π(θ|_n) by increasing m.

Beyond that,  <Ref> also tells us that the rate at which m increases the approximation quality of π(θ|_m,n) for π(θ|_n) depends on the slower of the two rates κ_1 and κ_2 that govern how quickly bias and variance of _m,n(θ, ) vanish.

In other words, it tells us that a biased estimator _m,n(θ,) for _m(θ,) will not adversely affect the approximation quality of π(θ|_m,n) as long as the order κ_1 of the bias is upper-bounded by that of the estimator's variance.





 §.§ Joint Asymptotics and Concentration (large n, large m)



While the previous result demonstrates that the behaviour of π_(|_m,n) is desirable for any fixed n as m→∞, it tells us nothing about how π_(|_m,n) behaves as m and n increase.

This is an important question since the convergence results in <Ref> only hold conditionally for a fixed , and therefore obscure dependence on a constant factor, say C(), that generally depends on both n and . 

Depending on the nature of C(), ensuring that π_(|_m,n)
remains a good approximation may become computationally infeasible: for instance, if C() ≍exp(n), the required number of samples m  would have to increase exponentially in n.


To investigate these interactions between n and m, we take m=m(n) in the remainder, and study π_(|_m,n) as n→∞.

To do this, we have to impose additional regularity conditions to control the behaviour of _n as n→∞.

To this end, define the limiting object

rCl
    (θ)     =    
    lim_n→∞𝔼_Y_1:n ∼P_0[ 
        n^-1
        
        _n(θ, Y_1:n)
    ].
    

We will assume that (θ) exists for all θ∈Θ,  and define its minimiser as

rCl
    θ_0     =    
    _θ∈θ (θ).
    

We impose the below assumptions to guarantee posterior concentration of  π(θ|_n).



For any δ > 0, (θ_0)< inf_θ-θ_0≥δ(θ).

There exists a function γ:Θ→ℝ_+ such that ∫_Θγ^2(θ)π(θ)θ < ∞, and for some N ∈ℕ and all n>N: 


    (i) _Y_1:n∼ P_0[{n^-1_n(θ, Y_1:n)-(θ)}^2]≤ Cγ^2(θ)/n;
    
    (ii)  
    
    ∫_Θ(θ)π(θ|_n)θ < ∞, and ∫_Θγ^2(θ)π(θ|_n)θ < ∞.
    





<Ref> part <ref> is a standard regularity condition for posterior concentration of π(θ|_n) (see, e.g., Condition 1 of  for a similar requirement), while part <ref> and the prior moment condition on γ^2 requires further discussion.

Paralleling how σ_n^2(θ) in <Ref> controls the pointwise variance in _m,n(θ, ) due to U_1:miid∼ P_θ and is used for the asymptotics in m (for fixed n),  γ^2(θ) in <Ref> controls the pointwise variance in _n(θ; Y_1:n) due to the variability in the observed data Y_1:n∼ P_0, and will drive asymptotics in n (with m=m(n)).

While part <ref> and the prior moment condition on γ^2(θ) are not standard and may be hard to verify in general, they are relatively mild conditions.

Whenever γ^2(θ) can be derived explicitly as for our examples in <Ref>, the prior  and posterior moment conditions can also be verified explicitly and numerically, respectively. 

For the case where Θ is compact or sup_θ∈Θ(θ) < ∞, both conditions hold automatically.


To extend our asymptotics to the case where n→∞ and m=m(n), we also require minor regularity conditions for the prior, as well as some refined bounds for the moment conditions in <Ref> that allow n →∞ and m = m(n). 

These requirements are jointly summarised in our last set of suppositions.




Suppose that the prior π is continuous on Θ, and that π(θ_0)>0. 

Further, assume that m is chosen as m=m(n), and that for  constants α_1, α_2, η_1, η_2 > 0, 


    (i) for ℬ_n={θ∈Θ:|(θ)-(θ_0)|≤α_1√(n)}, we have ∫_ℬ_nπ(θ)θ≥exp{-α_1√(n)};
    
    (ii)
    m^-κ_1∫_Θσ_n^2(θ)π(θ|_n)θ≍ m^-η_1 and m^-κ_2∫_Θσ_n^2(θ)π(θ)θ≍ m^-η_2 as n→∞;
    


where σ_n^2(θ) is the same variance proxy term as in <Ref>.



Almost all conditions in <Ref> are standard for posterior concentration.

The exception is <ref>, which extends some of the moment conditions in <Ref> for m=m(n) and n→∞.

On a technical level, <ref> ensures that the bounds  obtained by integrating out the bias and variance bounds in <Ref> do not diverge as n→∞ with m=m(n).





  §.§.§ Posterior Concentration


Together, our assumptions allow the derivation of novel concentration rates for the posterior π(θ|_m,n). 

They constitute the first results applicable to the setting we address, and provide direct information about the rate at which m=m(n) should scale with n as n →∞.

Our main result shows that if we want π(θ|_m,n) to place most of its posterior mass in regions of Θ that correspond to small (θ), we must choose m ≍ n^0.5 / min{η_1, η_2}. 



theoremconcentrationone

	If Assumptions <ref>-<ref>, are satisfied, then for any M_n→∞, it holds that
	
    P_0(∫_Θ|(θ)-(θ_0)|π(θ|_m,n)θ  > M_n/min{√(n),m^min{η_1, η_2}})⟶ 0.




The above shows that posterior concentration  occurs at rate min{√(n),m^min{η_1, η_2 }}.

In contrast, concentration of π(θ|_n)  occurs at rate √(n).

The implication is clear: for  π(θ|_m,n) to concentrate at the same rate as π(θ|_n), we require m≍ n^0.5 / min{η_1, η_2}.


Under a local  Hölder continuity condition for  in a neighbourhood of θ_0, concentration onto regions of low loss in <Ref> extends to open balls around θ_0. Such conditions are often encountered in the study of Gibbs posteriors, see for instance <cit.>, <cit.>, and <cit.>.


corollaryconcentrationtwo
    
	Under Assumptions <ref>-<ref>, and if on some open neighborhood of θ_0,  θ-θ_0≤ C|(θ)-(θ_0)|^α for some C>0 and α>0, then for any M_n →∞ and the sequence  S_n = {θ∈Θ:θ-θ_0>C(M_n/min{√(n),m^min{η_1,η_2}})^α},
    it holds  P_0-almost surely that
    rCl
        ∫_S_n π(θ|_m,n)≲M_n^-1.
        



While local  Hölder continuity condition is a mild requirement, it relates to the limiting object , and is thus generally hard to verify.

In most settings of practical interest however, it is unproblematic to simply assert  this condition for some α > 0, which implies that π(θ|_m,n) concentrates around θ_0 at rate min{n^α/2, m^αmin{η_1, η_2}}. 



  §.§.§ Asymptotic Normality and Posterior Mean


<Ref> implies that the contraction rate of π(θ|_n) is transferred to π(θ|_m,n) as long as m ≍ n^0.5·min{η_1, η_2}.

A similar result also holds for the posterior mean and  asymptotic posterior shape. 

To state this, let 𝖧=[ω·∇_θ^2(θ_0)]^-1, θ̂_n:=_θ∈Θ_n(θ), and 𝒩(θ;μ,Σ)  the density of a Gaussian distribution at θ with  mean μ and variance Σ.


corollarynewcorollary

    
    If Assumptions <ref>-<ref> are satisfied, and if ∫_Θ|π(θ|_n)-𝒩(θ;θ̂_n,n^-1𝖧^-1)|θ=o_p(1), then for θ̅_n=∫_Θθπ(θ|_m,n)θ, we have
    
    
        (i) ∫_Θ|π(θ|_m,n)-𝒩(θ;θ̂_n,n^-1𝖧^-1)|θ=o_p(1) if m^-min{η_1, η_2}→ 0 as n→∞;
        
        (ii)
        √(n)(θ̂_n-θ̅_n)=o_p(1) if √(n)m^-min{η_1, η_2}→ 0 as n→∞.
        
    
    


Part <ref> of the above result is especially useful for calibration of generalised Bayes: such techniques often rely on asymptotic normality of π(θ|_n) to scale _m,n by a value of ω in (<ref>) that  achieves approximate frequentist coverage. 

Part <ref> implies that if we choose m large enough, then the posterior mean of π(θ|_m,n)  asymptotically behaves like the computationally infeasible risk minimizer θ̂_n. 

In order for this behavior to emerge however, we must account for the simulation-induced bias of estimators constructed from π(θ|_m,n).

In particular, ensuring that this bias does not pollute the asymptotic distribution of θ̅_n requires choosing m ≫ n^0.5 / min{η_1, η_2}.





In the last part of this paper, we show that our theory is sufficiently flexible to cover a wide range of settings. 

This includes the MMD-based loss _m,n^k (see <Ref>), 
for which one can apply our results with κ_1 = κ_2 = η_1=η_2 = 1, which implies that the computationally feasible scaling m ≍√(n)  is sufficient to achieve the standard √(n)-concentration rate for π(θ|_m,n^k).

In contrast, a guarantee of √(n)-concentration with the KDE-based loss _m,n^h (see <Ref>) makes κ_1, κ_2, η_1, η_2 depend on the kernel bandwidth, and ultimately requires choosing m ≍ n^2.25 or larger.

Similarly, for the losses _m,n^β induced by β-divergences (see <Ref>) we have min{η_1, η_2} = 0.2, so that m≍ n^2.5  achieves √(n)-concentration.





§ APPLICATIONS




We apply our theory to Gibbs measures based on losses introduced in <Ref>: the β-divergence, the MMD, and a KDE-based likelihood estimator.

We focus on verifying  the most crucial conditions, namely Assumptions <ref>, <ref> part <ref>, and <ref> part <ref>.

We summarise our findings in Lemmas <ref> and <ref>, and conclude that _m,n^k and _m,n^h  require m≍√(n) and m≍ n^2.25 respectively for √(n)-concentration. 

We empirically verify this on simulated data using  the g-and-k model and a copula model.  

For the β-divergence, we find min{η_1, η_2} = 0.2, so that m≍ n^2.5 is required for standard √(n)-contraction rates (see <Ref>).







 §.§ Maximum Mean Discrepancy (MMD)


If the kernel k used to construct _n^k and _m,n^k is bounded and characteristic so that the corresponding MMD satisfies the identity of discernibles, all key assumptions hold.

This is unrestrictive, and true for most kernels of interest; including Gaussian and Matérn kernels.

lemmammd

	If the kernel k:𝒴^2 →ℝ is characteristic, and if sup_y, y'k(y,y') ≤ K, then for _m,n^k in (<ref>),  <Ref> holds for σ^2(θ) = 16K^4 + 2K, and κ_1 = κ_2 = 1.
    
    If it additionally  holds that  is independently and identically distributed, then we also have that <Ref> part <ref> holds with γ^2(θ) = 16K^4 + 2K and C=1.







Under the conditions of the above result, _m,n^k, _n^k, and the corresponding limiting loss ^k are all uniformly bounded by K, and both σ_n^2(θ) and γ^2(θ) are constant functions. Thus, <Ref> is actually verifiable by elementary arguments.

Further,  <Ref> part <ref> yields η_1 = κ_1 = 1 and η_2 = κ_2 = 1.

Hence, <Ref> implies that the difference between π(θ|_m,n^k) and π(θ|_n^k) vanishes at rate m^-1 for any fixed  n.

Lastly,
<Ref> implies that m ≍√(n) is sufficient for achieving the standard √(n) concentration rate.




  §.§.§ MMD-based Copula Model Inference


To verify the practical relevance of our theory, we provide numerical demonstrations using MMD-based inference in a copula model based on the loss _m,n^k proposed in <cit.>, details of which are provided in  <Ref>.



We compute generalised posteriors based on  _m,n^k, and compare the accuracy of posterior moments.

In particular, we generate  n=100 observations from a Student's t copula model with 10 degrees of freedom and dependence parameter θ = -0.1.

Inference is performed on θ only, holding the degrees of freedom fixed at 10.

Across 100 replications of this procedure, we compare mean and variance of the distribution π(θ|_m,n^k) for different choices of m in <Ref>. 

Our results show that posterior means become more accurate as m increases. Similarly, posterior variances tend to shrink as m increases. 

As our theory predicts,  improvements are substantial for smaller m, but negligible once  m≥√(n) = 10. Indeed, there is little difference in the boxplots for posterior means and variances between m=10 and m=20.






 §.§ Kernel Density Estimation (KDE) for Intractable Likelihood


When the likelihood  is intractable so that p_θ() cannot be evaluated exactly but we can sample  U_1:miid∼P_θ for each θ∈Θ, <cit.> propose to estimate the likelihood via KDE.

This amounts to computing π(θ|_m,n^h) based on _m,n^h as introduced in <Ref> to approximate the exact Bayes posterior π(θ|) = π(θ|_n) with _n(θ, ) = -log p_θ().

To show that our theory applies to this setting, we first establish conditions under which <Ref> holds.


lemmakde


Assume that for an isotropic kernel k
and a bandwidth h>0, we have


    (i)
        C_j(θ):=C_jsup_y∈𝒴|∂ ^j/∂ y^jp_θ(y)|<∞ for C_j<∞ and j∈{0,1,2,3};
    (ii)
        for every θ∈Θ, there exist δ(θ)>0 so that inf_y∈𝒴p_θ(y)≥δ(θ);
    (iii)
       ∫ k(v)v=1, k(v)=k(-v), ∫ k^2(v) v<∞, and ∫ v^2 k(v) v<∞.


Then, as h→0, and m→∞ with mh→∞, it holds for any finite n that 

rCl
_U_1:m iid∼ P_θ[logp̂_θ()]-logp_θ()
   ≤   max{C_2(θ), C_3(θ)}/δ(θ)nh^2; 


_U_1:m iid∼ P_θ[ {logp̂_θ()-logp_θ()}^2]
   ≤    n/mh{C_2(θ)+C_0(θ)h}/δ(θ)^2+nh^4max{C_1(θ), C_2(θ)}^2/δ(θ)^2.




Additional assumptions on p_θ and the data-generating process are needed to verify <Ref> <ref>.

For instance, if observations are sampled independently and are identically distributed according to P_0, we can simplify <Ref> by noting that

rCl
    _Y_1:niid∼ P_0[{n^-1logp_θ(Y_1:n)-𝖫(θ)}^2]= _Y_1:niid∼ P_0{n^-1logp_θ(Y_1:n)}=_Y_1∼P_0{logp_θ(Y_1)}/n.
    

Hence, we can take γ^2(θ) = _Y_1∼ P_0{log p_θ(Y_1)}, so that <Ref> <ref> simplifies to requiring that the Bayes posterior-averaged variance and expectation of the log likelihood is finite so that ∫_Y_1∼ P_0{log p_θ(Y_1)}π(θ|)dθ <∞ and ∫𝔼_Y_1∼ P_0{log p_θ(Y_1)}π(θ|)dθ <∞.

This generally depends on p_θ and P_0, but will generally hold if we have a second moment for log p_θ(Y_1) in Y_1 ∼ P_0, a condition commonly required for proving concentration of generalised posteriors <cit.>.



Given this, <Ref><ref> implies that π(θ|_m,n^h) behaves as π(θ|) so long as nh^2 → 0, and n/(mh)→0.

More specifically, it implies that for m=m(n)→∞ and as n→∞, 

rCl
    ∫_Θ|
    π(θ|_m,n^h)
    -
    π(θ|)
    |θ≲max{n(mh)^-1, nh^2}.
    


Here, nh^2→0 ensures that the bias of _m,n^h(θ, ) = logp̂_θ() vanishes, and that h^2≪ n^-1/2. 

To  control the variance, we also need n/(mh)→0. 

For m=n^a and h=n^-b, a,b>0, the bias condition implies b>1/2, and the variance condition requires a>1+b. 

As a result, we require m ≍ n^2.25 for π(θ|_m,n^h) to concentrate at rate √(n)
(see Appendix <ref> for details).





	



  §.§.§ KDE-based Inference for G-and-K-Model 


The scaling required for π(θ|_m,n^h) to concentrate at the √(n)-rate calls the computational feasibility of KDE-based likelihood estimators into question.

We demonstrate this by conducting inference on the g-parameter of the intractable g-and-k distribution (see <Ref> for additional details).

Throughout our experiments, we compare against an almost exact  benchmark posterior obtained using the computationally prohibitive likelihood approximation strategy in <cit.>.

<Ref> plots the results for a single run with different choices for m and n. 

The results reflect our theory: regardless of n, the approximation quality of  π(θ|_m,n^h) is poor for smaller choices of m.  

In particular, π(θ|_m,n^h) is biased and has higher variance than the target posterior.

While bias and variance of π(θ|_m,n^h) decrease as m increases, the difference with the target is notable even for large m.

Perhaps counter-intuitively, the results also illustrate that for small n,  m has to be increased by a larger multiple than for larger n to achieve the same performance. 

This points to a limitation of our theoretical guarantees: they are asymptotic in nature and thus ignore small-sample effects.[

  To illustrate this, inspect (<ref>). 

  Choosing the MSE-optimal bandwidth rate of h≍ m^-2/5, we obtain κ_1 = 3/5, so that the non-asymptotic version of the upper bound in (<ref>) for n=25 and m=100· n becomes Cmax{ 0.23,  0.05} for some C.
  
  Since the righthand side is also trivially upper bounded by 2, this is vacuous if C > 2/0.23 > 9. 
  
  While this upper bound will eventually become meaningful as n and m=m(n) diverge, ensuring that this term is small for an arbitrary C may require choosing a very large value of m when n is also small.
 ]





To supplement the qualitative findings of <Ref>, we replicate these computations 100 times for different choices of n and m, and summarise the results in <Ref> for the model parameter g, and in <Ref> for the remaining parameters.

To keep comparisons fair across different n and m, the MCMC algorithm used for posterior computations uses a total of 10 million simulations of the full dataset from the model for each setting.[This means that for any m, we use 10· n million  simulations from the g-and-k distribution.]

The findings validate our theory, and show that m has to be chosen much larger than n for reliable inference.
























§ DISCUSSION



Frequently, a target Gibbs measure π(θ|_n) is computationally infeasible, and has to be approximated through the measure π(θ|_m,n), which substitutes an intractable loss _n for the tractable approximation _m,n.

While the use of π(θ|_m,n) is already common in contemporary practice, our results provide the first generic and comprehensive asymptotic theory for such objects.

Our key finding is that their approximation quality depends crucially on the quality of the estimator _m,n. 

In particular, our theory shows that if η_1 and η_2 denote the rates at which posterior-averaged biases and variances of _m,n vanish, the approximation error of π(θ|_m,n) vanishes at rate m^-min{η_1, η_2} for any fixed n; see Assumption <ref> in Section <ref> for the definition of η_1,η_2 and further details.

To ensure that the approximation quality of π(θ|_m,n) stays constant as n increases, our findings further indicate that we need to scale the number of simulations as m ≍ n^0.5 / min{η_1, η_2}.

Note that all results apply regardless of the choice for ω in (<ref>)—a crucial tuning parameter often referred to as the learning rate in the literature on generalised Bayes.


While we focused on working out the technical details for three illustrative examples, the theory derived in the previous section has much wider applicability.

For generalised Bayesian methods, it applies immediately for the setting of losses constructed from α- and γ-divergences <cit.>, total variation  <cit.>, Wasserstein distances <cit.>, and sliced distances <cit.>.


Further, the application of our theory to kernel density estimation (KDE) suggests that various newer black-box density estimation methods, like sequential neural likelihood <cit.>, for inference with intractable likelihoods may also require a similar level of draws in order to accurately approximate an intractable posterior. To the best of our knowledge, no formal results exist on the accuracy of such likelihood estimators.


Lastly, our results establish the first general asymptotic theory for posteriors constructed from MCMC methods based on biased likelihood  estimators, and thus extend existing results on the large-sample behavior of pseudo-marginal MCMC methods <cit.>.

Surprisingly, the results show that a biased likelihood estimator may be preferable to an unbiased alternative if computational effort is a major concern.

In particular, our theoretical results suggest that instead of using unbiased likelihood estimators with slowly decaying variances, we should consider using biased estimators whose variance decays more quickly.

Given that existing methods for pseudo-marginal methods do not allow for bias, our contributions could help to liberate the field from the restrictive requirement of unbiased likelihood estimators.





	1.0 

 
	chicago
	
	
	
	

    tocsectionAppendix 
    

PART: 
Appendix 
    
    

	

§ FURTHER DETAILS




 §.§ <Ref>




Importantly, while the moment bounds <ref> and <ref> in <Ref> are pointwise for each θ∈Θ, they are tied together via σ^2_n(θ).

Hence, σ^2_n(θ) can be thought of as an asymptotic variance proxy.

To illustrate this, define the sequence of functions {σ^2_m,n}_n ∈ℕ, m∈ℕ  with σ^2_m,n:Θ→ℝ_+ via the identities

rCl
    _U_1:m iid∼ P_θ [ {_m,n(θ, y_1:n)-_U_1:m iid∼ P_θ[_m,n(θ, y_1:n)]}^2 ] 
    =
    σ^2_m,n(θ)m^-κ_2.
    


We can now see that for any fixed n,  σ^2_n(θ) is an asymptotic upper bound over {σ^2_m,n(θ)}_m∈ℕ.

To make this bound as tight as possible, we could choose σ^2_n(θ) so that  σ^2_n(θ)m^-κ_2≍σ^2_m,n(θ)m^-κ_2 as m→∞.

This also implies an upper bound on the bias term in <ref>: decomposing the left hand side of (<ref>) into a squared bias and variance term yields 

rCl
    (
        _U_1:m iid∼ P_θ [ _m,n(θ, y_1:n) ] - _n(θ, )
    )^2
    +
    _U_1:m iid∼ P_θ [ {_m,n(θ, y_1:n)-_n(θ, )}^2 ],
    

and since the second term in this decomposition is positive, an application of Jensen's inequality reveals that σ_n(θ)m^-κ_1 with κ_1 = 0.5 ·κ_2 is a valid  asymptotic upper bound on the bias term.[This is true since we can always increase σ_n(θ) if necessary to ensure σ_n(θ) > 1.]

In practice however, κ_1 = 0.5 ·κ_2 is typically a much slower rate  than what reasonable estimators achieve.

As our convergence results depend on the minimum between κ_1 and κ_2, we thus state our assumptions  in terms of both rates.



 §.§ Numerical Demonstrations with Copula Model



Motivated by the impact of outliers on likelihood-based copula inference, <cit.> propose a point estimator for the dependence parameter in copula models by minimising the function θ↦_n^k(θ, ).


While more robust to misspecification than the negative log likelihood, this function has no closed form, and estimation must instead be based on the approximate version _m,n^k. 

Following <cit.>, we can also use _m,n^k for a generalised Bayesian version of this estimation problem. 

Here, our theory predicts that as long as m≍√(n), the resulting posterior  π(θ|_m,n^k) should be close to π(θ|_n^k).

In practice, this means that we expect  π(θ|_m,n^k) to provide an excellent approximation to π(θ|_n^k),  even if  m  is quite small.



To verify the predictions of our theory, we consider a simple bivariate copula model. 

Suppose that we observe the bivariate random variables Y_i=(Y_i,1,Y_i,2)^⊤iid∼ F from the joint cumulative probability function (CDF) with marginal CDFs F_j.  Then, Sklar's Theorem states that F can be modeled by some copula C:[0,1]^2→ℝ_+, where for  y=(y_1,y_2),

rCl
    F(y_1,y_2)=C(F_1(y_1),F_2(y_2)) = C(u_1,u_2),   where  u_j=F_j(y_j)  for  j=1,2.
 

Based on this, a common approach for modelling such bivariate observed data y_1:n is via some parametric copula function: C_θ:[0,1]^2×Θ→ℝ_+, where the parameter θ measures the dependence between y_i,1 and y_i,2 for each i=1,2,… n.

In practice, the marginals F_j are of course unknown, so that one must estimate them. 

A common approach is to do this via the empirical marginal CDFs F_n,j(t) = 1/n∑_i=1^n{y_i,j≤ t}.

One can then construct a corresponding set of transformed observations û_i=(û_1,i,û_2,i)=(F_n,1(y_1,i),F_n,2(y_1,i))^⊤, where the hat notation  signifies that û_j,i is an estimate of u_j,i=F_j(y_j,i).


Given {û_i}_i=1^n, pseudo-likelihood estimation of θ can be carried out via the copula probabilisty density function c_θ
<cit.>.

While pseudo-likelihood inference is extremely common in copula modelling, the method is known to be non-robust to outliers and contamination bias, due to the behavior of copula densities near the extremes of the support. 

To make copula-based inference more robust, <cit.> thus propose to conduct inference on θ by minimising the MMD between the empirical distribution of {û_i}_i=1^n (often called the empirical copula process) and the copula model.

For some kernel function k:[0,1]^2×[0,1]^2→ℝ_+, this gives rise to the loss

rCl
    _n^k(θ, )=
    𝔼_U ∼c_θ, U' ∼c_θ[ 
        k(U, U')
    ]
    -
    2/n∑_i=1^n 
    𝔼_U ∼c_θ[ 
        k(U, û_i)
    ].
    
    

Unfortunately, the integrals featuring in _n^k are intractable even for the simplest copulas.

However, it is straightforward to approximate them via Monte Carlo: 


defining 𝒰_2 as the distribution of a bivariate uniform distribution on [0,1]^2, we can simulate m sets of n independently and identically distributed draws from the copula model via the pushforward u_j^θ = (u_1,j^θ,u_2,j^θ)^⊤iid∼C_θ∘𝒰_2 for j=1,… m, and construct 




rCl
    _m,n^k(θ, )=1/m^2∑_j'=1,j=1^mk(u_j^θ,u_j'^θ)-
    2/mn∑_j=1^m∑_i=1^nk(u_j^θ,û_i).
    



    

 §.§ Discussion on m and h in the KDE


While Lemma <ref> shows that Theorem <ref> is satisfied by the KDE likelihood estimator for any fixed n, it is not directly clear what rates on h and m are needed to control the bias and variance of the estimated likelihood as n→∞. 

To this end, we parameterise m≍ n^a with a>0, and choose the bandwidth h≍ n^-b, b>0.

Given these choices, let us explore how we can control the terms nh^2 and n(mh)^-1 as n→∞. 

Clearly, we must choose a and b such that mh=n^a-b→∞, and nh^2≍ n^1-2b→0, which requires that a>b>1/2. 

However, if we additionally want to ensure that the desired √(n)-convergence is satisfied, <Ref> implies that we must choose these terms so that 


    nh^2≪ n^-1/2,  and  n/(mh)≪ n^-1/2.
 

Using the parameterisations m=n^a and h=n^b,  we can ensure this is the case by analysing the values of a and b for which it holds that


    1-2b=-1/2 and 1+b-a=-1/2.


The first equation implies b=3/4, and plugging this into the second  yields a=2.25. 

Hence, in order to control the bias of the approximation and ensure that the posterior concentrates at the canonical √(n) parametric rate, we must choose the bandwidth like h≪ n^-3/4, and simulate m≫ n^2.25 data points at each MCMC iteration. 

Taken together, this showcases that to obtain accurate results using the KDE likelihood estimator, one must simulate a very large number of data points—resulting in a rather large computational cost. 




 §.§ Numerical Demonstrations with G-and-K Model



To numerically illustrate that choosing m to be sufficiently large is crucial, we use the  g-and-k distribution <cit.>.

This is a commonly used example for likelihood-free inference <cit.>. 

It has no closed form density, and is instead defined in terms of its quantiles 


    Q{z(p);θ}   = a + b[1 + c1-exp{-gz(p)}/1 + exp{-gz(p)}] {1 + z(p)^2}^kz(p),


where p denotes the g-and-k distribution's quantile of interest,  z(p) the quantile of a standard normal, and  θ = (a,b,c,g,k)^⊤ the model parameters. 

Throughout, we follow common practice and fix c = 0.8 (see <cit.> for a justification).  

While the likelihood function for the g-and-k distribution is not available in closed form, it can be approximated numerically to arbitrary precision using the strategy outlined in <cit.>.

Though far more cumbersome than approximations based on _m,n^h, this provides near-exact computations of the Bayes posterior and can be used as the benchmark against which we compare π(θ|_m,n^h).


The parameter g is the most challenging parameter for inference in the g-and-k model.

Accordingly, in the main text we limit attention to inference on g, and defer results for the remaining parameters to <Ref>. 

Throughout, we generate n data points from the g-and-k model, and then compare π(θ|_m,n^h) against near-exact computations of the Bayes posterior following the likelihood approximation strategy in <cit.>.

Figure <ref> plots these posteriors for g for different m and n. The results reflect our theory: regardless of n, the approximation quality of  π(θ|_m,n^h) is very poor for smaller choices of m.  

In particular, π(θ|_m,n^h) is biased and has higher variance than the target Bayes posterior.

While bias and variance of π(θ|_m,n^h) decrease as m increases, the difference with the target is notable even for large m.



Though instructive, <Ref>  only shows posterior inferences based on a single data set.

To ascertain that the inflated variance and the bias are a problem across all possible data sets, we compare  posteriors computations across 100 datasets and different choices of m and n.



MCMC is used for posterior computations in each case so that 10 million model simulations of the full dataset are used in total for each setting.


Thus, regardless of m and n the number of MCMC iterations is adjusted to ensure a roughly constant computational burden: 
for example, for m=50 (m=100)  we run MCMC for 200,000  (100,000) iterations.

<Ref> summarises the results of this analysis for g, and <Ref> contains results on the remaining parameters.

These simulations validate the findings of <Ref> and of <Ref> by showing that m has to be chosen much larger than n for reliable posterior inferences.




 §.§ Additional Numerical Results: G-and K

    
This section gives results on the KDE-based posterior for the a,b,k parameters in the g-and-k model (see Tables <ref>, <ref>, <ref>, <ref>). 

The results show that the other parameters are less impacted by the choice of m than g.

This is due to the fact that their impact on the density is not as   nonlinear as the impact of g, which controls the skewness of the density. 

That being said, our theory still applies, and the numerical results show that higher samples sizes still require larger choices for m to obtain more accurate inferences. However, the results are not as extreme as those for the g parameter.

	

	
	


	

	
	


	


§ TECHNICAL LEMMAS

	We now state and prove several lemmas used to help prove the main results of the paper. The first result uses 
	Assumptions  <ref>-<ref> to establish a concentration bound for π^(θ|_n). We note that this result is not for the posterior π(θ|_m,n) but the idealised version π^(θ|_n).  To further simplify our notation in the appendix we write _n(θ) for  _n(θ,).





 §.§ <Ref>



    
	Under Assumptions <ref>-<ref>, for any M_n→∞ as n→∞, 
	
    P^_0( ∫_Θ |(θ)-(θ_0)|π(θ|_n)θ>M_n/n^1/2)≲1/M_n.



The result follows a similar set of arguments to those in <cit.>, <cit.>, as well as Theorem 1 in <cit.>.  Critically, however, unlike these results we do not assume uniform boundedness of expectations for _n(θ). 
To simplify certain arguments that follow, we assume without loss of generality that (θ_0)=0. Note that this is without loss of generality as we can always redefine _n(θ) (and (θ)) as _n(θ):=_n(θ)-(θ_0) (and (θ):=(θ)-(θ_0)) without changing the resulting posterior distributions we study. 




Firstly, since 0=(θ_0)≤(θ) for all θ∈Θ, 
it follows that 

    ∫_Θ|(θ)-(θ_0)|π(θ|_n)θ=|∫_Θ{(θ)-(θ_0)}π(θ|_n)θ|=|∫_Θ(θ)π(θ|_n)θ-(θ_0)|.
Hence, for any δ>0, 


    P_0( ∫_Θ |(θ)-(θ_0)|π(θ|_n)θ>δ)	
    =P_0( |∫_Θ(θ)π(θ|_n)θ-(θ_0)|>δ).






Now, note that 

∫_Θ(θ)π(θ|_n)θ≤   	∫_Θn^-1_n(θ)π(θ|_n)θ+∫_Θ|n^-1_n(θ)-(θ)|π(θ|_n)θ.
	From  <Ref>, with probability at least 1-δ, 
|n^-1_n(θ)-(θ)|≤ Cγ(θ)/√(n)δ by Markov's inequality. By  <Ref>, ∫_Θγ(θ)π(θ|_n)θ=C_1<∞, so that we have 


	∫_Θ(θ)π(θ|_n)θ≤    	∫_Θn^-1_n(θ)π(θ|_n)θ+C_1/√(n)δ.


Adding n^-1KL(ρπ)≥0, and noting that π(θ|_n) solves (<ref>) and satisfies ∫γ(θ)π(θ|_n)θ<C for some constant C>0 (<Ref>), we have 

rCl
	∫_Θ(θ)π(θ|_n)θ    ≤    	∫_Θn^-1_n(θ)π(θ|_n)θ+n^-1KL{π(·|)π}+C_1/√(n)δ
 

     ≤   inf_ρ∈𝒫(Θ):∫γ(θ)ρ(θ)<C{∫_Θn^-1_n(θ)ρ(θ)+n^-1KL(ρπ)}+C_1/√(n)δ .    


By Assumption <ref>,  it holds with probability at least 1-δ for the set of all ρ∈𝒫(Θ) for which ∫γ(θ)ρ(θ)<C simultaneously that


    ∫_Θ n^-1_n(θ)ρ(θ)≤∫_Θ(θ)ρ(θ)+1/√(n)δ∫_Θγ(θ)ρ(θ)≤∫_Θ(θ)ρ(θ)+C_2/√(n)δ.


Using this, we combine (<ref>) and (<ref>)  to conclude that 




    ∫_Θ(θ)π(θ|_n)θ≤inf_ρ∈𝒫(Θ):∫γ(θ)ρ(θ)<C{∫_Θ(θ)ρ(θ)+n^-1KL(ρπ)}+C_1+C_2/√(n)δ.


Clearly, (θ)≤(θ_0)+|(θ)-(θ_0)|, so that the RHS of the above is bounded above by  

    ∫_Θ(θ)π(θ|_n)θ≤inf_ρ∈𝒫(Θ):∫γ(θ)ρ(θ)<C{∫_Θ{(θ)-(θ_0)}ρ(θ)+n^-1KL(ρπ)}+C_1+C_2/√(n)δ .



Let 
    ρ_n(θ):=
    π(θ)/Π(ℬ_n)    if θ∈ℬ_n,
    0    else,
 

 denote the prior measure normalized by the prior magnitude of the set ℬ_n:={θ∈Θ:{(θ)-(θ_0)}≤α_1/√(n)}, and note that, by Assumption <ref>, 
Π(ℬ_n)≥exp{-α_2√(n)}. Further, we note tha ρ_n(θ) is in the stated set, since ∫_Θγ(θ)^2π(θ)θ<C by Assumption <ref>. It then follows that 

    ∫_Θ{(θ)-(θ_0)}ρ_n(θ)≤α_1/√(n),   n^-1KL{ρ_n||π}=-logΠ(ℬ_n)≤α_2/√(n).

The above terms deliver an upper bound for the first term in equation (<ref>), so that for 0<δ<1,  

    ∫_Θ{(θ)-(θ_0)}π(θ|_n)θ≤α_1+α_2+C_1+C_2/√(n)δ,
with probability at least 1-δ. From the above, and (<ref>), we have that 

P_0( ∫_Θ|(θ)-(θ_0)|π(θ|_n)θ>α_1+α_2+C_1+C_2/√(n)δ)≤δ.
Taking δ={(2C+α_1+α_2)}/M_n then yields 

	P_0( ∫_Θ|(θ)-(θ_0)|π(θ|_n)θ>M_n/√(n),)≤(2C+α_1+α_2)/M_n≲M_n^-1,
as desired. 
		

	


 §.§ <Ref>


	We use the following general result to bound the difference between the infeasible posterior π(θ|_n), in equation (<ref>), and π(θ|_m,n). As far as we are aware, the following result is similar to, but distinct from, arguments in the literature that characterize sub-Gaussian behavior of random variables; see, e.g., Theorem 2.1 in <cit.>.
	
	
    
		Let Z be a positive, scalar-valued random variable with mean μ that satisfies [(Z-μ)^2]≤ B b^2 for some B,b>0. Then, for any 0≤λ≤ 1/b, 
		
    [exp(-λ Z)]≤exp(-λμ)+B(λ b)^2.

	
	
	
		For any Z≥0, a Taylor expansion of exp(-λ Z)  around Z=μ yields 
		 
			exp(-λZ)     = exp(-λμ)∑_n=0^∞(-1)^nλ^n/n!(Z-μ)^n=τ(Z)+ℰ(Z),
		where τ(Z) is the truncation of the series
		
    τ(Z):=exp(-λμ){1-λ(Z-μ)}

		and
		
    ℰ(Z):=exp(-λζ)λ^2/2!(Z-μ)^2
is the Lagrange remainder term with ζ a constant between Z and μ. 
		
		Boundedness of exp(-x) over x≥0, implies that for any λ≥0 and Z≥0
		
    ℰ(Z)≤λ^2(Z-μ)^2.

		Consequently,
		
			exp(-λZ) = τ(Z)+ℰ(Z)≤exp(-λμ){1-λ(Z-μ)}+λ^2(Z-μ)^2.
		By hypothesis, [(Z-μ)^2]≤ B b^2<∞ for some B,b>0. Take expectations of both sides to obtain  
		
			[exp(-λZ)]    ≤exp(-λμ){1-λ(Z-μ)}+λ^2(Z-μ)^2
   ≤exp(-λμ)+B(λb)^2.
		
	
	
	

§ PROOFS OF MAIN RESULTS



Recall that _m,n(θ,) depends on a sequence of random variables U_1:miid∼P_θ, where the dependence on U_1:m was previously subsumed via the dependence on the index m. 	Define g_n^ω(θ):= exp{-ω·_n(θ)} and recall the generalized posterior:
	
    π^(θ|_n)=exp{-ω·_n(θ)}π(θ)/∫_Θexp{-ω·_n(θ)}π(θ)θ≡g_n^ω(θ)π(θ)/∫_Θ g_n^ω(θ)π(θ)θ,
	which is the (generalized) posterior that would result if we were able to analytically calculate _n(θ). In addition, let us simplify our notation and subsume the dependence of _m,n(θ,) on  via the index n and write _m,n(θ) for _m,n(θ,).



 §.§ <Ref>


*
 


We first demonstrate that, uniformly over Θ,
		
			𝔼_U_1:m[exp{-ω_m,n(θ)}]= exp{-ω_n(θ)}[1+C{σ_n^2(θ)/m^κ}]+C{σ^2_n(θ)/m^κ_2},
		where the dependence of _U_1:m on θ and  is suppressed for notational simplicity, and where we remind the reader that U_1:miid∼P_θ. 
  
		From <Ref>,
        
    _U_1:m{_m,n(θ)- _n(θ)}≲σ_n^2(θ)/m^κ_1, and   _U_1:m[{_m,n(θ)- _U_1:m[_m,n(θ)]}^2]≲σ^2_n(θ)/m^κ_2,
 where, throughout the remainder of this proof, the notation ≲ is in relation to sequences in m for fixed n (as the result is true for fixed n). 
        
        Take Z=_m,n(θ), μ= _U_1:m_m,n(θ), λ=ω, and apply
		Lemma <ref> to obtain the following upper bound for _U_1:m[exp{-ω_m,n(θ)}]: with b=m^-0.5 ·κ_2, m≥1, and B=Cσ^2_n(θ), for some C>0 such that _U_1:m[{_m,n(θ)- _n(θ)}^2]≤ Bb^2,
    
       
		
			𝔼_U_1:m[exp{-ω_m,n(θ)}]   ≤exp{-ω_U_1:m[_m,n(θ)]}+B(λb)^2
   ≤exp[-ω_n(θ)+C{ωσ_n^2(θ)/m^κ_1}]+C{σ^2_n(θ)/m^κ_2}
   ≲exp{-ω_n(θ)}[1+{ωσ_n^2(θ)/m^κ_1}]+σ^2_n(θ)/m^κ_2.
		
  
   
  Where the first inequality uses <Ref> part <ref>, and the second holds due to the fact that
  exp{ωσ_n^2(θ)/m^κ_1}=1+exp{ωσ_n^2(θ)/m^κ_1}-1≲ 1+{ωσ_n^2(θ)/m^κ_1}.
		
		Define g^ω_n(θ)=𝔼_U_1:m[exp{-ω_m, n(θ)}]. From equation  (<ref>) and the definitions of g_n^ω(θ) and g^ω_n(θ),
  
		
			|g_n^ω(θ)-g^ω_n(θ)|≤g^ω_n(θ) C{σ^2_n(θ)/m^κ_1}+C{σ_n^2(θ)/m^κ_2}.
		By Assumption <ref>, for any n≥1, ∫_Θ g^ω_n(θ)π(θ)θ<∞, so that, applying (<ref>), 
		
			|∫_Θ g_n^ω(θ)π(θ)θ-∫_Θ g^ω_n(θ)π(θ)θ|   ≤∫_Θ|g_n^ω(θ)-g^ω_n(θ)|π(θ)θ
   ≲1/m^κ_1∫_Θ σ_n^2(θ)g^ω_n(θ)π(θ)θ+1/m^κ_2∫_Θσ^2_n(θ)π(θ)θ
   = 1/m^κ_1∫_Θ g^ω_n(θ)π(θ)θ∫_Θ σ_n^2(θ)π(θ|_n)θ
   +1/m^κ_2∫_Θσ^2_n(θ)π(θ)θwhere the second line follows from equation (<ref>),  and the equality in the third line from
		
        rCl
            ∫σ_n^2(θ) g_n^ω(θ) π(θ) θ    =    
            ∫g_n^ω(θ) π(θ)dθ·( 
                
                    ∫σ^2_n(θ) g_n^ω(θ) π(θ)θ
                    ∫g_n^ω(θ) π(θ)θ
            )
            

                =    
            ∫g_n^ω(θ) π(θ)θ·∫σ^2_n(θ)  π(θ|_n)θ.
            
               
        
        
		  Notice that <Ref> ensures that ∫_Θσ^2_n(θ)θ^2π(θ)θ and ∫_Θσ^2_n(θ)θ^2π(θ|_n)θ are finite for any n. 
  
        Since it also holds for any ξ∈ [0, 2] that ∫_Θσ^2_n(θ)θ^ξπ(θ)θ≤∫_Θσ^2_n(θ)θ^2π(θ)θ and ∫_Θσ^2_n(θ)θ^ξπ(θ|_n)θ≤∫_Θσ^2_n(θ)θ^2π(θ|_n)θ, we may apply the exact same arguments as before to conclude that 
        
		
			|∫_Θ θ^ξg_n^ω(θ)π(θ)θ-∫_Θ θ^ξg^ω_n(θ)π(θ)θ|   ≲1/m^κ_1∫_Θθ^ξσ_n^2(θ)π(θ|_n)θ
   +1/m^κ_2∫_Θθσ^2_n(θ)π(θ)θ.
		
  
        The above display can be re-expressed as
        
		
			|∫_Θ θ^ξg_n^ω(θ)π(θ)θ-∫_Θ θ^ξg^ω_n(θ)π(θ)θ|≲m^-min{κ_1, κ_2}.
		

  


Aplying equation (<ref>) twice with ξ=0, once for the numerator and once for the denominator, yields
		
			|∫_Θ g_n^ω(θ)π(θ)θ-∫_Θ g^ω_n(θ)π(θ)θ|/∫_Θ g^ω_n(θ)π(θ)θ≲m^-min{κ_1, κ_2}.
		 










	Write  π(θ|_m,n)-π_(θ|) as 
		
			π(θ|_m,n)-π(θ|_n)=    g_n^ω(θ)π(θ) /∫_Θ g_n^ω(θ)π(θ)θ-g^ω_n(θ)π(θ)/∫_Θ g^ω_n(θ)π(θ)θ 

			=   { g_n^ω(θ)-g^ω_n(θ)}π(θ)/∫_Θ g^ω_n(θ)π(θ)θ ∫_Θ g^ω_n(θ)π(θ)θ/∫_Θ g_n^ω(θ)π(θ)θ 

			   -g_n^ω(θ) π(θ)(1/∫_Θ g_n^ω(θ) π(θ)θ-1/∫_Θ g_n^ω(θ) π(θ)θ),
		 and apply the triangle inequality to obtain
		
			|π(θ|_m,n)-π(θ|_n)|    ≤|g_n^ω(θ)-g^ω_n(θ)| π(θ)/∫_Θ g_n^ω(θ)π(θ)θ
   +|∫_Θ g_n^ω(θ)π(θ)θ-∫_Θ g^ω_n(θ)π(θ)θ|/∫_Θ g_n^ω(θ)π(θ)θ π(θ|_n).
		 
        
        Multiplying by θ^ξ, integrating both sides, and using the exact same line of arguments as above, leads to
		
			   ∫_Θθ^ξ|π(θ|_m,n)-π(θ|_n)| θ≤1/∫_Θ g_n^ω(θ)π(θ)θ ∫_Θθ^ξ|g_n^ω(θ)-g_n^ω(θ)|π(θ) θ
   +|∫_Θg_n^ω(θ)π(θ)θ-∫_Θg_n^ω(θ)π(θ)θ|/∫_Θg_n^ω(θ)π(θ)θ∫_Θθ^ξπ(θ|_n)θ.
   
  By (<ref>)-(<ref>), we have that 
  
    1/∫_Θg_n^ω(θ)π(θ)θ∫_Θθ^ξ|g_n^ω(θ)-g_n^ω(θ)|π(θ)≲m^-min{κ_1, κ_2}/∫_Θg_n^ω(θ)π(θ)θ≲ m^-min{κ_1, κ_2}
where the second inequality follows since 0<∫_Θ g_n^ω(θ)π(θ)θ<∞ by Assumption <ref>. 
  For the second term in (<ref>), first note that by Assumption <ref>, ∫_Θθ^ξπ(θ|_n)θ<∞, hence, applying again (<ref>) shows that this term is also bounded by a multiple of m^-min{κ_1, κ_2}. 		This completes the proof.
	




 §.§ <Ref>


*
 

We apply <Ref> and <Ref> to derive a concentration bound for π(θ|_m,n). This result relies on asymptotics in (m,n) and so we express m=m(n), but suppress the dependence on n for notational simplicity. As in the proof of <Ref>, we assume without loss of generality that (θ_0)=0, since we can always redefine _n(θ) (and (θ)) as _n(θ):=_n(θ)-(θ_0) (and (θ)=(θ)-(θ_0)) without altering the posterior distribution we study. 


Define κ = min{κ_1, κ_2}. Using similar arguments to those of Lemma <ref>,

    ∫_Θ |(θ)-(θ_0)|π(θ|_m,n)θ=|∫_Θ(θ)π(θ|_m,n)θ-(θ_0)|.
Similarly, with probability at least 1-δ, we have 

∫_Θ(θ)π(θ|_m,n)θ≤    ∫_Θ(θ)π(θ|_m,n)θ+∫_Θ|(θ)-n^-1_n(θ)|π(θ|_m,n)θ
≤    ∫_Θn^-1_n(θ)π(θ|_m,n)θ+1/√(n)δ∫_Θγ(θ)π(θ|_m,n)θ.

We can decompose the second term into two components that both converge to zero. 

The first of these vanishes by Theorem <ref>, and the second by Assumption <ref> via

	1/√(n)∫_Θγ(θ)π(θ|_m,n)θ≤    1/√(n)∫_Θγ(θ)|π(θ|_m,n)-π(θ|_n)|θ+1/√(n)∫_Θγ(θ)π(θ|_n)θ
≤    C/√(n)m^κδ+C/√(n)δ,	
where the second line utilizes a moment bound for γ(θ) in Assumption <ref><ref>. 
Focusing on ∫_Θ n^-1_n(θ)π(θ|_m,n)θ, we obtain the upper bound

∫_Θn^-1_n(θ)π(θ|_m,n)θ   =∫_Θ(θ)π(θ|_m,n)θ+∫_Θ{n^-1_n(θ)-(θ)}π(θ|_m,n)θ
   ≤∫_Θ(θ){π(θ|_m,n)-π(θ|_n)}θ+∫_Θ(θ)π(θ|_n)θ
   +∫_Θ{n^-1_n(θ)-(θ)}π(θ|_m,n)θ.


Consider the last term in (<ref>). Similarly as in (<ref>), by <Ref>, Assumptions <ref> and <ref>, and with probability at least 1-δ, we have the bound

∫_Θ{n^-1_n(θ)-(θ)}π(θ|_m,n)θ   ≤1/√(n)δ∫_Θγ(θ)π(θ|_m,n)θ
   ≤1/√(n)δ∫_Θγ(θ){π(θ|_m,n)-π(θ|_n)}θ
   +1/√(n)∫_Θγ(θ)π(θ|_n)θ
   ≤C/√(δ^2m^2κn)+C/δ√(n). 

Consider the second term in (<ref>), and note that, again with probability at least 1-δ, 


∫_Θ(θ)π(θ|_n)θ   =∫_Θ{(θ)-n^-1_n(θ)}π(θ|_n)θ+∫_Θn^-1_n(θ)π(θ|_n)θ
   ≤1/√(n)δ∫_Θγ(θ)π(θ|_n)θ+∫_Θn^-1_n(θ)π(θ|_n)θ.


Using similar arguments as those in (<ref>) and <Ref>, for some C>0, with probability at least 1-δ


∫_Θn^-1_n(θ)π(θ|_n)θ   ≤C/√(n)δ.


Recall that, from Assumption <ref>, for m=m(n) and n→∞, m^-κ_2∫_Θσ^2_n(θ)(θ)π(θ)θ≍ m^-η_1 and m^-κ_1∫_Θσ^2_n(θ)(θ)π(θ|_n)θ≍ m^-η_2 for some η_1,η_2>0. Under this assumption, using similar arguments as those that led to (<ref>) in the second half of the proof of Theorem <ref> we obtain  the bound

rCl
          
	∫_Θ (θ)|π(θ|_m,n)-π(θ|_n)|θ

	    ≲    m^-min{η_1,η_2}+m^-min{η_1,η_2} ·∫_Θ g_n^ω(y_1: n |θ) π(θ) d θ/∫_Θ g_n^ω(y_1: n |θ) π(θ) d θ{∫_Θ(θ) π(θ|_n) d θ} 

	    ≍      m^-min{η_1,η_2}.



Collecting the terms in (<ref>)-(<ref>), we have that, with probability at least 1-δ, for some C>0,

	∫_Θ(θ)π(θ|_m,n)θ≲(θ_0)+1/δ(C/√(nm^2κ)+3C/√(n)+C/m^min{η_1,η_2}).
 
For 0<δ<1, we then have that 

	P^_0{ ∫_Θ|(θ)-(θ_0)|π(θ|_m,n)θ>1/δ(C/√(nm^2κ)+C/√(n)+C/m^min{η_1,η_2})}   ≤δ.
	Noting that √(n)≥√(n m^2κ), and taking δ=C/M_n, 

    P^_0{∫_Θ |(θ)-(θ_0)|π(θ|_m,n)θ > M_n/min{√(n),m^min{η_1,η_2}}}≤ C/M_n≲ 1/M_n.





 §.§ <Ref>


*

 
The result follows by appropriately modifying the argument used in, e.g., Corollary 2 of <cit.>. Define m_η,n:=min{m^min{η_1,η_2},√(n)}. Firstly, recall that, for A⊂Θ, 

    Π(θ∈ A|_m,n)=∫_Aπ(θ|_m,n)θ.

By Markov's inequality, and Theorem <ref>, we have that


    Π{|(θ)-(θ_0)|≥ M_n/m_η,n|_m,n}≤[∫_Θ |(θ)-(θ_0)|π(θ|_m,n)θ]/M_n/m_η,n≲ M_n^-1.

Let δ>0 be such that {θ∈Θ:θ-θ_0≤δ}⊂ V:={θ∈Θ:θ-θ_0≤ C|(θ)-(θ_0)|^α}. By Assumption <ref>, for any δ>0 such that θ-θ_0>δ, there exists ϵ such that 0<ϵ≤|(θ)-(θ_0)|. Let ϵ_n={M_n/m_η,n}→0 as n→∞, and let n be large enough such that ϵ_n:={M_n/m_η,n}<ϵ, which implies that 

    {θ∈Θ :|(θ)-(θ_0)|≤ϵ_n}⊆{θ∈Θ:θ-θ_0≤δ}⊂ V.

Hence,  from equation (<ref>) in the proof of Theorem <ref>, for any M_n→∞,  

1-M_n^-1≳Π{|(θ)-(θ_0)|≤M_n/m_η,n|_m,n}=    Π{C|(θ)-(θ_0)|^α≤C{M_n/m_η,n}^α|_m,n}
≥    Π{θ-θ_0≤C{M_n/m_η,n}^α|}.





 §.§ <Ref>


*



    By the triangle inequality 

∫_Θ|π(θ|_m,n)-𝒩(θ;θ̂_n,𝖧^-1/n)|θ   ≤∫_Θ|π(θ|_m,n)-π(θ|_n)|θ
   +∫_Θ|π(θ|_n)-𝒩(θ;θ̂_n,𝖧^-1/n)|θ.
Taking (θ)=1 in the proof of <Ref>, see equation (<ref>), we can show that ∫_Θ |π(θ|_m,n)-π(θ|_n)|θ=O_p(m^-min{η_1,η_2}), which implies that 
    ∫_Θ |π(θ|_m,n)-𝒩(θ;θ̂_n,𝖧^-1/n)|θ≤ O_p(m^-min{η_1,η_2})+o_p(1). Hence, for m^-min{η_1,η_2}→0, the first stated result follows. 

To obtain the second result, decompose θ̅_n as 

    θ̅_n=∫θπ(θ|_m,n)θ= ∫θ{π(θ|_m,n)-π(θ|_n)}θ+∫θπ(θ|_n)θ.

Applying Theorem <ref>, we have

θ̅_n=   O_p(m^-min{ η_1, η_2 })+∫θπ(θ|_n)θ
=   O_p(m^-min{ η_1, η_2 })+∫θ{π(θ|_n)-N(θ;θ̂_n,n^-1𝖧^-1)}θ+∫θN(θ;θ̂_n,n^-1𝖧^-1)θ
=   θ̂_n+O_p(m^-min{ η_1, η_2 })+∫θ{π(θ|_n)-N(θ;θ̂_n,n^-1𝖧^-1)}θ.

Applying the change of variable t=√(n)(θ-θ̂_n),

θ̅_n-θ̂_n=O_p(m^-min{ η_1, η_2 })+1/√(n)∫t {π(t|_n)-N(t;0,𝖧^-1)}t

Since the total variation distance is invariant to changes of location and scale, from the hypothesis of Corollary <ref>, we have that ∫ t {π(t|_n)-N(t;0,𝖧^-1)} t=o_p(1), and we can re-arrange the above equation as 

    √(n)(θ̅_n-θ̂_n)=O_p(√(n)m^-min{η_1, η_2 })+o_p(1).
So long as √(n)m^-min{η_1, η_2 }→0 as n→∞, the stated result follows. 







§ APPLICATIONS




 §.§ β-divergence




Recalling the definitions of _n^β and _m,n^β in <Ref>, it is straightforward to verify <Ref>.

In particular, it is immediately clear that _U_1:miid∼ P_θ[ _m,n^β(θ, ) ] = _n^β(θ, ), so that part <ref> holds automatically for any κ_1 > 0 and σ_n^2(θ) = 0.

For part  <ref>, note that

rCl
    _U_1:m iid∼ P_θ[ {_m,n^β(θ, ) - _n^β(θ, ) }^2 ]
        =    
    n^2/m
    _U_j ∼P_θ[
        p_θ^2β(U_j)
    ]
     = 
     n^2/m
     ∫p_θ^1+2β(u)du,
     

so that <Ref> part <ref> holds for σ_n^2(θ) = n^2 ∫ p_θ^1+2β(u) u and κ_2 = 1.


To verify <Ref>, notice that if the data is independently and identically distributed, part <ref> is easy to establish: since  𝔼_Y_1:niid∼P_0[ _n^β(θ; Y_1:n)] = (θ) = n · d_β(p_0, p_θ), the expression on the left of <ref> is equivalent to Var_Y_1:niid∼P_0(n^-1_n(θ, Y_1:n)). Further,   

rCl
    Var_Y_1:n iid∼P_0(n^-1_n(θ, Y_1:n))
        =    
    ( 1+ 1/β)^2
    n^-1
    Var_Y_i ∼P_0(
        p_θ(Y_i)
    ),
    

so that <Ref> part <ref>  holds for C =  ( 1+ 1/β)^2 
and γ(θ) = √(Var_Y_i∼ P_0( p_θ(Y_i) )).

Notice that the moment conditions on γ^2 in <Ref> now become equivalent to the mild assumption that Var_Y_i∼ P_0( p_θ(Y_i) ) has finite prior and posterior expectation.



For <Ref>  part <ref>, recall that σ_n^2(θ) = n^2 ∫ p_θ^1+2β(u) u and κ_2 = 1.

Therefore, provided that ∫(∫ p_θ^1+2β(u) u)π(θ)θ <∞,  we require m^-1 · n^2 ≍ m^-η_2.

Taking m ≍ n^ξ, this leads to the requirement 2-ξ  = - ξη_2.

To make posterior computation simultaneously as accurate and inexpensive as possible, we wish to pick the smallest ξ that still achieves the fastest possible posterior concentration rate of π(θ|_n,m^β). 

Since m ≍ n^ξ, <Ref> tells us that such a ξ must satisfy n^ξη_2≍ n^0.5, which leads to the relation ξη_2 = 0.5.

Plugging this into the first relation, we thus find that ξ = 2.5, which yields η_2 = 0.2.





 §.§ <Ref>


 *
 


We provide the argument for
<Ref> only, as the arguments for <Ref> part <ref> follow along the exact same lines.

First, decomposing the mean square error into squared bias and variance, we find

rCl
           𝔼_U[ 
        (_m,n(θ, ) - _n(θ, ) )^2
    ]
    

         =    
     {
     𝔼_U[ 
        _m,n(θ, )
    ] - _n(θ, )
    }^2
    +
    𝔼_U[ 
        {_m,n(θ, ) - 𝔼_U[_m,n(θ, )]}^2
    ].
    

Since both terms in the last line are non-negative, this means that to complete our proof, it suffices to show that

rCl
    𝔼_U[ 
        (_m,n(θ, ) - _n(θ, ) )^2
    ]
        ≤   
    16K^4/m.
    


To show that this bound indeed holds, note that by definition, _m,n^k(θ, ) = ^2(P_m,θ, P_n) and _n^k(θ, ) = ^2(P_θ, P_n), where P_n = 1/n∑_i=1^nδ_y_i and P_m,θ = 1/m∑_j=1^m δ_u_j for u_1:miid∼ P_θ.

Using this and |a^2 - b^2| = (a+b)|a-b|, it is  straightforward to show that

rCl
          
    |_m,n^k(θ, ) - 
    _m^k(θ, )
    |
    

         =    
     ((P_θ, _n) +
		(P_m,θ, _n))
		|(P_θ, P_n) -
		(P_m,θ, P_n)|
		

		    ≤    
		4K ·(P_θ, P_m,θ ), 

where we have used sup_y,y'k(y,y') ≤ K for the first term, and the triangle inequality for the second term in the last line.

Notice that the triangle inequality is applicable because k was assumed to be characteristic, so that the MMD is a proper metric on 𝒫().

We now complete the proof by squaring the above on both sides, taking the expectation, and then using that _U[ ^2(P_θ, P_m,θ )] ≤K^2/m (cf. Lemma 5 of <cit.>).





 §.§ <Ref>


*


	First, for each y∈, we have 
	
    logp_θ(y)-log p_θ(y)=1/p_θ(y){p_θ(y)-p_θ(y)}-1/21/ζ^2{p_θ(y)-p_θ(y)}^2,
 for some real number ζ between p_θ(y) and p_θ(y). Note that, since from the definition of ζ, we have that ζ>0, and  
	
		logp_θ(y)-logp_θ(y)≤   1/p_θ(y){p_θ(y_i)-p_θ(y)}
≤    δ(θ)^-1{p_θ(y_i)-p_θ(y)}, 	
	
	where we recall that δ(θ)<inf_y∈𝒴p_θ(y). Since (<ref>) holds for all y∈𝒴, we can use this, and the fact that 
	
    p̂_θ(y_1:n)=∏_i=1^n{1/mh∑_l=1^mk_h(y_i-u_l)}.
 to obtain the stated results. To this end, we first derive the mean and variance of p_θ(y_i) wrt u_i=(u_1,…,u_m), where we recall that u_i⊥ u_j for all j i. 
	


  §.§.§ Term _U_i{p_θ(y_i)}
 From the definition of p_θ(y_i), and since u_iiid∼P_θ, 
	
		_U{1/mh∑_i=1^mk(u_i-y/h)}	-p_θ(y)=1/h∫p(z|θ)k(u-y/h)u-p_θ(y)
	Use the change of variables u-y=hv, and a second-order Taylor series expansion to obtain
	
		   1/h∫p(u|θ)k(u-y/h)u
=   ∫p(y+hv|θ)k(v)v
=   ∫{p_θ(y)+∂p_θ(y)/∂y hv +2^-1∂^2 p_θ(y)/∂y^2 (hv)^2+6^-1∂^3 p(y|θ)/∂y^3 (hv)^3}k(v)v,
	 for some intermediate value y. Considering the last term, we have that, from our assumption on  ∂^3 p_θ(y)/∂ y^3, 
	
    6^-1∫∂^2 p(y|θ)/∂ y^2 (hv)^3 k(v) v≤h^3/6sup_y∈𝒴|∂^2 p(y|θ)/∂ y^2|∫ v^3k(v) v= C_3(θ).
Noting that ∫ v k(v) v=0 by assumption, and that ∫ v^2 k(v)<∞, we can apply the above into equation (<ref>) to obtain
	
		_U{1/mh∑_i=1^mk(u_i-y/h)}	-p_θ(y)   ≤h^3C_3(θ)+h^2∫v^2 k(v)v/4|∂^2 p_θ(y)/∂y^2|.
	Writing C_2(θ)=∫ v^2 k(v) v/6sup_y∈𝒴|∂^2 p_θ(y)/∂ y^2|, we have the individual upper bound, 
	
		_U{1/Lh∑_j=1^Lk(z_j-y/h)}	-p_θ(y)   ≤{C_3(θ) ∨C_2(θ)} h^2,
	since h^2≫ h^3 for all m large enough, and which is valid for any y∈𝒴.
	


  §.§.§  Term var_U_i{p_θ(y)}
	
	
 
	A similar change of variable and a first-order Taylor expansion yields 
	
		var_U_i{p_θ(y)}   =1/mh^2[∫p(y+hv|θ)k(v)^2hv-{h∫p(y+hv|θ)k(v)v}^2]

		   =1/mh^2[∫{p_θ(y)+∂p(y|θ)/∂y(hv)}k(v)^2hv-{h∫p(y+hv|θ)k(v)v}^2]

		   = 1/mhp_θ(y)∫k^2(v)v+C_1(θ)1/m-1/m{∫p(y+hv|θ)k(v)v}^2

		   ≤p_θ(y)κ_2/mh+C_1(θ)/m+m^-1p_θ(y)

		   ≤ p_θ(y){1+κ_2}/mh+m^-1C_0(θ)
   ≤{C_0(θ)+C_1(θ)h}/(mh),
	where C_1(θ)=sup_y∈𝒴|∂ p_θ(y)/∂ y|κ_2, κ_2=∫ k^2(v)  v, C_0={1+κ_2}sup_y∈𝒴p_θ(y) 
	and since 
	{∫ p(y+hv|θ)k(v) v}^2≤ p_θ(y)+hC_1(θ)∫ vk(v) v= p_θ(y).

	
	Now, taking the expectation on both sides of (<ref>), and applying (<ref>), we have, uniformly for y∈𝒴,  
	
		_U_i{	logp_θ(y)-logp_θ(y)}≤    δ(θ)^-1_U_i{ p_θ(y)-p_θ(y)}
≤    δ(θ)^-1{C_2(θ)∨C_3(θ)}h^2
	
	In addition, 
	
		_U_i{	logp_θ(y)-logp_θ(y)}^2≤    δ(θ)^-2_U_i{ p_θ(y)-_U_ip_θ(y)}^2
+   δ(θ)^-2[{_U_ip_θ(y)-p_θ(y)}]^2
	Applying equations (<ref>) and (<ref>) to the above then yields 
	
		_U_i{	logp_θ(y)-logp_θ(y)}^2≤   1/mh{C_0(θ)+C_1(θ)h}/δ(θ)^2+{C_2(θ)∨C_3(θ)}^2h^4/δ(θ)^2.
		
	
	
	Overall Mean: Since each U_i is independent, and since (<ref>) is independent of y, the first result follows.  
	

	
	Overall Variance: Since each U_i is independent, 
	
		_U{logp_h(|θ,U)-logp(|θ)}^2=   ∑_i=1^n_U_i{logp_θ(y_i)-_U_ilogp_θ(y_i)}^2
+   ∑_i=1^n[
		_U_i{logp_θ(y_i)-logp(y_i|θ)}]^2
	
	For each i, we can now apply the bounds in (<ref>) and (<ref>) to obtain  the upper bound
	
		_U{logp_h(|θ,U)-logp(|θ)}^2≤    n{C_0(θ)+C_1(θ)h}/δ(θ)^2mh+δ(θ)^-2{C_2(θ)∨C_3(θ)}^2h^4n.
	


	
	
