arXiv:1601.07911v2 [math.ST] 22 Nov 2016

On asymptotic validity of naive inference with an
approximate likelihood
Helen Ogden
University of Southampton, UK
Abstract
Many statistical models have likelihoods which are intractable: it is impossible
or too expensive to compute the likelihood exactly. In such settings, a common
approach is to replace the likelihood with an approximation, and proceed with
inference as if the approximate likelihood were the exact likelihood. In this
paper, we describe conditions on the approximate likelihood which guarantee
that this naive inference with an approximate likelihood has the same first-order
asymptotic properties as inference with the exact likelihood. We investigate the
implications of these results for inference using a Laplace approximation to
the likelihood in a simple two-level latent variable model, and using reduced
dependence approximations to the likelihood in an Ising model on a lattice.
Keywords: Intractable likelihood, Ising model, Laplace approximation, Latent variable model

1

Introduction

For many models, it is impossible or infeasibly expensive to evaluate the likelihood function, typically because it involves a high-dimensional sum or integral.
In such cases, a common approach is to find an approximation L̃(.) to the likelihood L(.), and to use L̃(.) in place of L(.) to conduct inference about the
parameters of the model.
For instance, one could construct a point estimate of the parameters by
maximizing the approximate likelihood, and form confidence intervals based on
the curvature of the approximated log-likelihood about its maximum. From a
Bayesian perspective, an approximate posterior π̃(θ; y) ∝ L̃(θ; y)π(θ) could be
formed by substituting the approximate likelihood in place of the exact likelihood.
Such an approach is commonly used in practice. In latent variable models, where the likelihood is an integral over the latent variables, naive inference
using a Laplace approximation to the likelihood is used in both maximum likelihood (Pinheiro and Bates, 1995; Bates et al., 2015) and Bayesian (Rue et al.,
2009) settings. In Markov random field models, where the likelihood involves
an intractable normalizing constant, inference is often conducted by substituting an approximation to this normalizing constant in place of the exact

1

normalizing constant into the expression for the likelihood (Friel et al., 2009;
Tjelmeland and Austad, 2012).
In this paper, we provide conditions under which the naive approach is
asymptotically justified. Under these conditions, the approximate maximum
likelihood estimator is consistent and has the same asymptotic normal distribution as the exact maximum likelihood estimator, hypothesis tests based on
the approximate likelihood remain valid, and in Bayesian analysis the distance
between the approximate posterior and the true posterior shrinks to zero.
Douc et al. (2004) show that the approximate maximum likelihood estimator
will have the correct asymptotic normal distribution provided that the error in
the log-likelihood,
ǫn (θ) = log L̃n (θ) − log Ln (θ),
tends in probability to zero as n → ∞, uniformly in θ. We argue that this
measure is too strict in many practical examples, in which ǫn (θ) grows rapidly
with n and yet the inference remains asymptotically valid. Our conditions are
based instead on ∇θ ǫn (θ), the error in the approximation to the score function.
We provide two examples to demonstrate how the conditions may be used in
practice. The first is a simple two-level latent variable model, with mn repeated
observations for each of n items. We deduce the rate at which mn must grow
with n in order for the Laplace approximation to give asymptotically valid
inference. If mn grows with n at a slower rate, the estimator remains consistent,
but loses efficiency relative to the exact maximum likelihood estimator, and
naively constructed confidence intervals have lower than nominal coverage.
The second example is an Ising model on an m × m lattice, with the class of
reduced dependence approximations (Friel et al., 2009) used to approximate the
likelihood. For parameter values associated with weak dependence, we show that
the reduced dependence approximation may be used to obtain asymptotically
valid inference at cost polynomial in m, in contrast with the exponential cost
of computing the likelihood exactly.

2
2.1

Asymptotic validity of approximate likelihood
inference
Setup and notation

Consider a sequence of models indexed by n, with common parameter θ ∈ Θ ⊆
Rp . Write ℓn (θ; y) for the log-likelihood given observed data y under model
n, un (θ; y) = ∇θ ℓn (θ; y) for the corresponding score function and Jn (θ; y) =
−∇Tθ ∇θ ℓn (θ; y) for the observed information. We sometimes drop the data y
from the notation for convenience. Suppose that the data were generated from
the model for some θ0 ∈ Θ, and that as n → ∞, the amount of information
provided by the data about the parameter grows at some rate rn , such that
Jn (θ0 ) = O(rn ) in probability. This includes the case with n independent
replications as a special case, with rn = n, but we also wish to allow for more
complex settings. Write I(θ) for the Fisher information matrix, chosen such
that J¯n (θ) = rn−1 Jn (θ) → I(θ) in probability.
Let ℓ̃n (.; y) be an approximate log-likelihood, which in general may be any
function of the parameters θ, which will be used in place of ℓn (.; y). Write θ̂n

2

and θ̃n for the estimators maximizing ℓn (θ) and ℓ̃n (θ) respectively. Suppose that
ℓn (θ) and ℓ̃n (θ) are both three times differentiable, and write ũn (θ) = ∇θ ℓ̃n (θ)
and J˜n (θ) = −∇Tθ ∇θ ℓ̃n (θ) for the approximate score and information.
Write
ǫn (θ) = ℓ̃n (θ) − ℓn (θ)
for the pointwise error in the log-likelihood,
δn (θ) = k∇θ ǫn (θ)k = kũn (θ) − un (θ)k
for the absolute error in the score, and
γn (θ) = k∇Tθ ∇θ ǫn (θ)k = kJn (θ) − J˜n (θ)k
for the absolute error in the
we
P
P observed information matrix. For concreteness,
use the L1 norms kak = i |ai | for a vector a, and kAk = maxj { i |Aij |} for
a matrix A, although the same results hold for any choice of norms.
Write δn∞ (S) = supθ∈S δn (θ) for the uniform error in the score over any set
S ⊆ Θ, and let δn∞ = δn∞ (Θ). Similarly, define γn∞ (S) = supθ∈S γn (θ) and
γn∞ = γn∞ (Θ). For any θ0 ∈ Θ, write Bt (θ0 ) = {θ ∈ Θ : kθ − θ0 k ≤ t} for the
ball of radius t about θ0 .

2.2

Approximate maximum likelihood inference

First, we describe sufficient conditions to ensure that θ̃n is consistent. The
proofs of all results are given in the appendix.
We will assume some standard regularity conditions on the model. Writing ūn (θ) = rn−1 un (θ), and ū(θ) for the limit as n → ∞, we assume that
supθ∈Θ kūn (θ)
R − ū(θ)k → 0 in probability. We assume ū(.) is such that, for
any ǫ > 0, θ:d(θ,θ0)≥ǫ kū(θ)k > ū(θ0 ) = 0. These conditions are stronger than
necessary, and we expect the same result to hold in many other situations where
the exact maximum likelihood estimator is consistent.
Theorem 1. Suppose δn∞ = op (rn ) as n → ∞. Then θ̃n → θ0 in probability, as
n → ∞.
We now give conditions to ensure that θ̃n retains the same limiting distri−1/2
bution as θ̂n . Since θ̂n − θ0 is Op (rn ), this is equivalent to finding conditions
−1/2
under which θ̂n − θ̃n is op (rn ). The following lemma bounds the distance
between θ̃n and θ̂n in terms of the error in the score function near θ0 .
Lemma 1. Suppose that δn∞ = op (rn ), and that there exists t > 0 such that
δn∞ {Bt (θ0 )} = op (an ). Then θ̃n − θ̂n = op (an rn−1 ).
1/2

Applying Lemma 1 with an = rn
result.

leads directly to the asymptotic normality

Theorem 2. Suppose that δn∞ = op (rn ), and that there exists t > 0 such that
1/2
δn∞ {Bt (θ0 )} = op (rn ). Then
rn1/2 (θ̃n − θ0 ) → N (0, I(θ0 )−1 )
in distribution, as n → ∞.
3

It is also desirable for hypothesis tests constructed by using the approximate likelihood in place of the exact likelihood to have the correct asymptotic
distribution.
Consider testing the hypothesis H0 : θ ∈ ΘR , where ΘR ⊂ Θ and dim(ΘR ) =
q. Write θ̂nR for the restricted maximum likelihood estimator, and Λn = 2{ℓn (θ̂n )−
ℓn (θ̂nR )} for the likelihood ratio statistic. The approximate version of the likelihood ratio test statistic is Λ̃n = 2{ℓ̃n(θ̃n ) − ℓ̃n (θ̃nR )}, where θ̃nR is the restricted
approximate likelihood estimator.
Under the same conditions that were used to show that θ̃n has the correct
limiting distribution, plus a bound on the error in the information around θ0 ,
Λ̃n is asymptotically equivalent to Λn under H0 .
Theorem 3. Suppose that δn∞ = op (rn ), and that there exists t > 0 such that
1/2
δn∞ {Bt (θ0 )} = op (rn ) and γn∞ {Bt (θ0 )} = op (rn ). Then, under H0 , Λ̃n − Λn =
op (1).
The Wald and score test statistics, Wn and Sn , are asymptotically equivalent to the likelihood ratio test, so under H0 , all three statistics have limiting
distribution χ2p−q . Under the conditions of Theorem 3, I(θ0 ) is consistently
estimated by rn−1 J˜n (θ̃n ), so the approximate Wald and score test statistics W̃n
and S̃n are also asymptotically equivalent to Λn .

2.3

Approximate Bayesian inference

We now consider the approximate posterior
π̃(θ|y) ∝ L̃n (θ; y)π(θ),
where we suppose that the prior is such that log π(.) is three times differentiable. Under the same conditions that were used to show asymptotic correctness of maximum likelihood inference, the total variation distance between the
approximate and exact posteriors,
Z
1
π̃(θ|y) − π(θ|y) dθ,
dT V {π̃(θ|y), π(θ|y)} =
2 Θ
tends to zero as n → ∞.
Theorem 4. Suppose that δn∞ = op (rn ), and that there exists t > 0 such that
1/2
δn∞ {Bt (θ0 )} = op (rn ) and γn∞ {Bt (θ0 )} = op (rn ). Then
dT V {π̃(θ|y), π(θ|y)} = op (1)
as n → ∞.
If the Bernstein-von Mises Theorem holds for the exact posterior distribution, under the conditions of Theorem 4, it will also hold for the approximate
posterior distribution π̃(θ|y). In that case, credible regions formed from the
approximate posterior distribution will also be valid confidence sets.

4

2.4

Adjusted approximate likelihood inference
1/2

If Lemma 1 holds for some an > rn , then θ̃n will still be a consistent estimator,
but might not have the same limiting distribution as θ̂n . The approximate
likelihood may still be useful in practice, provided that the inference is adjusted
accordingly.
The sandwich information matrix (Godambe, 1960) is Gn (θ) = In (θ)Hn (θ)−1 In (θ)
where Hn (θ) = var{∇θ ũn (θ)} and In (θ) = E{J˜n (θ)}. Under suitable regularity
conditions,
−1
s1/2
)
n (θ̃n − θ0 ) → N (0, Ḡ(θ0 )
in distribution, where sn is the rate of convergence of θ̃n , chosen such that
Gn (θ0 ) = O(sn ), and Ḡ(θ) = limn→∞ s−1
n Gn (θ).
Composite likelihood estimators (Lindsay, 1988) also have this type of asymptotic behaviour, and many methods which have been proposed to adjust inference using a composite likelihood could also be used to adjust the inference with
an approximate likelihood, provided that θ̃n is a consistent estimator. For example, Varin et al. (2011) describe various methods to approximate the variance
of a composite likelihood estimator, and list some modification to the composite
likelihood ratio test statistic designed to ensure that the resulting test statistic
has approximately χ2p−q distribution. From a Bayesian perspective, the adjustments proposed by Pauli et al. (2011) and Ribatet et al. (2012) to posterior
distributions based on composite likelihoods could also be used in the context
of approximate likelihood inference.

3

Examples

3.1
3.1.1

Laplace approximation to the likelihood in a simple
latent variable model
A two-level model

Suppose Yi ∼ Binomial(m, pi ), where logit(pi ) = bi and bi ∼ N (0, θ2 ), for
i = 1, . . . , n. The likelihood is
Z Y
n
{logit−1 (bi )}yi {1 − logit−1 (bi )}m−yi φ(bi ; 0, θ)db
Ln (θ) =
=

Rn i=1
n
Y

L(i) (θ),

i=1

where
L(i) (θ) =

Z

∞

−∞

{logit−1 (bi )}yi {1 − logit−1 (bi )}mn −yi φ(bi ; 0, θ)dbi

and φ(.; µ, σ) is the N (µ, σ 2 ) density function.
If we take a Laplace approximation L̃n (θ) to the likelihood, it is intuitively
clear that m = mn will have to grow with n to give valid inference as n → ∞.
It is less obvious whether any choice of mn that grows with n will give valid inference, or whether mn needs to grow with n at some minimum rate. Rue et al.
(2009) suggest that any mn which grows with n will suffice, conjecturing that
5

“the error rate” is [number of latent variables / number of observations], although they note that this rate is not established rigorously. In this case, the
error rate refers to the error in approximating π(θ|y) with π̃(θ|y), found by using a Laplace approximation to the likelihood. The Integrated Nested Laplace
Approximations proposed by Rue et al. (2009) are based on this π̃(θ|y), with
further approximations used to approximate the marginal posterior distribution
of each component of θ, if p > 1. In this example θ is a scalar, so the Integrated
Nested Laplace Approximation to the posterior distribution is precisely π̃(θ|y).
3.1.2

Theoretical analysis

The factorization of the likelihood allows us to study the error for each item
ǫ(i) (θ) = ℓ(i) (θ) − ℓ(i) (θ) separately, and then combine the errors. In the supplementary materials, we show that for each fixed θ ∈ Θ,
δ(i) (θ) = k∇θ ǫ(i) (θ)k = Op (m−2
n ),
Pn
so δn (θ) ≤ i=1 δ(i) (θ) = Op (nm−2
n ).
However, the conditions are in terms of uniform rather than pointwise errors.
−1/2
Since δ(i) (θ) is maximized at a point θ∗ = Op (mn ), and decreasing for all
∗
θ>θ ,
∞
= δ(i) (θ∗ ) = Op (mn−1/2 ),
δ(i)
−1/2

so δn∞ = Op (nmn ).
The amount of information that the data provides about θ is bounded for
fixed n as mn → ∞, by the available information on θ given the value of each
bi . So if mn → ∞ as n → ∞, then rn = O(n), and since δn∞ = op (rn ), θ̃n will
be consistent.
Provided that θ0 6= 0, choose any fixed t ∈ (0, θ0 ). Then, for sufficiently
large n,
δn∞ {Bt (θ0 )} = δn (θ0 − t) = Op (nm−2
n )
and
γn∞ {Bt (θ0 )} = γn (θ0 − t) = Op (nm−2
n ).

If mn grows at a rate faster than n1/4 , then δn∞ {Bt (θ0 )} = op (n1/2 ) and
γn∞ {Bt (θ0 )} = op (n), so the Laplace approximation to the likelihood will give
first-order correct inference.
3.1.3

Numerical demonstration

To verify these results numerically, we simulate 10000 realizations from the
model with θ0 = 0.5, various values of n between 1000 and 10000, and mn =
min{1, 5 + 4(na − 1000a)}, for a = 0.2, 0.25 or 0.3. The three choices of mn are
shown in Figure 1a.
A very accurate approximation to the likelihood may be obtained by using
adaptive Gaussian quadrature with 20 quadrature points to approximate each
of the univariate integrals L(i) (θ). We use this approximation as a proxy for the
exact likelihood L(θ).
As n → ∞, we have rn = O(n), but for smaller sample sizes E(kJn (θ̂n )k)
still grows with mn . E(kJn (θ̂n )k) may be approximated by r̂n = Ê(kJn (θ̂n )k),
where Ê(.) the sample mean over the 10000 realizations. The functional form
6

1.2
1.0

30

0.6

0.8

20
5 10

2000 4000 6000 8000

2000 4000 6000 8000
(a) mn

(b) r̂n

δn (θ0 )

1.1

0.01

0.03

1.3

0.05

1.5

0.07

−1/2

2000 4000 6000 8000

(c) RMSE(θ̃n )

(d) RMSE(θ̃n )/RMSE(θ̂n )

0.65

0.20

0.75

0.30

0.40

0.85

2000 4000 6000 8000

2000 4000 6000 8000

2000 4000 6000 8000

(e) Coverage of 90% confidence intervals

(f) dT V {π̃(θ|y), π(θ|y)}

Figure 1: Comparison of exact and approximate likelihood inference, for a twolevel model with mn observations on each of n items, where mn = min{1, 5 +
4(na − 1000a)}, for a = 0.2 (dashed lines), 0.25 (solid lines) and 0.3 (dotted
lines). RMSE(.) denotes root mean squared error.

7

−1/2

of mn was chosen to make r̂n δn (θ0 ) approximately constant in the a = 0.25
case, as shown in Figure 1b. The same quantity grows with n for a = 0.2 and
shrinks with n for a = 0.3.
Figure 1c shows the root mean squared error in the Laplace estimator, and
as expected, the estimator is consistent for all three choices of a. The root
mean squared error of the Laplace estimator divided by that of the maximum
likelihood estimator, shown in Figure 1d, grows in the a = 0.2 case, stays
approximately constant if a = 0.25, and shrinks towards 1 if a = 0.3.
Figure 1e shows the actual coverage of likelihood ratio type confidence intervals for θ, of nominal 90% coverage. The upper three lines show the coverage
of the intervals constructed using the exact likelihood, which have very close to
nominal coverage for each a. The lower three lines show the coverage of the approximate likelihood intervals, which decreases with n for a = 0.2, and increases
towards the nominal 90% level for a = 0.3.
Figure 1f shows the total variation distance between π̃(θ|y) and π(θ|y), with
prior π(θ) ∝ 1/θ. The distance between the approximate posterior and the
exact posterior grows in the a = 0.2 case, stays approximately constant if a =
0.25, and shrinks towards zero if a = 0.3. The behaviour for a = 0.2 refutes
the conjecture of Rue et al. (2009) that the error in the approximate posterior
distribution should shrink to zero provided that mn grows with n.

3.2
3.2.1

Reduced dependence approximation to the likelihood
in an Ising model
An Ising model

We consider a simple Ising model for n = rc variables yi ∈ {−1, 1}, arranged on
an r × c lattice, with parameters θ = (α, β), so that
pr(Y = y; θ) = {Zr,c (θ)}−1 exp{αV0 (y) + βV1 (y)},
P
P
j indicates that i and j
where V0 (y) = i yi , and V1 (y) = i∼j yi yj . Here i ∼P
have an edge between them in the lattice, and Zr,c (θ) = y∈{−1,1}n exp{αV0 (y)+
βV1 (y)} is a normalizing constant. The likelihood function L(θ; y) = pr(Y =
y; θ) depends on Zr,c (θ), and it is the computation of this normalizing constant
which makes evaluation of the likelihood function difficult. By using variable
elimination (e.g. Jordan, 2004), Zr,c (θ) may be computed at cost O(rc2min{r,c} ),
which remains infeasibly expensive if both r and c are large.
3.2.2

Reduced dependence approximations

Many methods for approximating Zr,c (θ) have been proposed. We will study
properties of inference using the reduced dependence approximations introduced
by Friel et al. (2009), a family of approximations controlled by a positive integer
(k)
tuning parameter, which we call k. The approximation for fixed k is Z̃r,c (θ) =
{Zk,c (θ)}r−k+1 /{Zk−1,c (θ)}r−k .
We consider the case r = c = m, using a reduced dependence approximation
(k)
(k)
at level k to approximate to the likelihood, giving L̃m (β) = {Z̃m,m (θ)}−1 exp{αV0 (y)+
2 m
βV1 (y)}. The exact likelihood may be computed at cost O(m 2 ), and the reduced dependence approximation at level k at cost O(m2 + km2k ). The aim

8

is to understand how k = km should vary with m to give asymptotically valid
inference as m → ∞. The error in the log-likelihood,
(k)
(k)
ǫ(k)
m (β) = ℓ̃m (β) − ℓm (β) = log Zm,m (β) − log Z̃m,m (β),

does not depend on the data y, so we do not need to consider the statements in
probability: all of the errors are deterministic in this case.
3.2.3

A special case where the exact likelihood is available

If α = 0, and the lattice has periodic boundary conditions, so that the top row
of variables are joined to the bottom row, and the left row joined to the right,
Kaufman (1949) provides a relatively simple expression for Zr,c (0, β), so that it
is possible to compute the likelihood exactly, even for large lattices.
We restrict the parameter space
√ to α = 0, and assume β ∈ [0, 0.43]. This
guarantees that β < βc = log(1 + 2)/2 ≈ 0.44, where βc is a critical value at
which the behaviour of the Ising model suddenly changes, so that for β > βc
large areas of all plus ones or all minus ones are observed. If β0 = βc , the
maximum likelihood estimator may not have a normal limiting distribution, so
our results do not apply to this case.
The information provided by the data about β grows at rate rm = m2 . In
(k)
the supplementary material, we show that if k → ∞ as m → ∞, then δm (β) =
O{m2 k exp(−bβ k)} + o(1), where bβ = 2 cosh−1 {−1 + cosh(2β)2 / sinh(2β)}.
(k)
For any choice of k which grows with m, supβ∈[0,0.43] δm (β) = o(m2 ), so the
approximate likelihood estimator will be consistent.
In order to meet the conditions of Theorem 2, k = km should be chosen
(k)
so that δm {Bt (β0 )} = o(m). For any t < βc − β0 , for sufficiently large m,
(k)
(k)
δm {Bt (β0 )} = δm (β0 + t) = O{m2 k exp(−bβ0 +t k)} + o(1). Since bβ is a
continuous function of β, any k such that m2 k exp(−bβ0 k) = o(m) will meet
this condition for sufficiently small t. This may be achieved by taking km =
−1
cβ0 log m, for any cβ0 > b−1
β0 . Figure 2a shows how bβ varies with β.
(k)

For any k → ∞ as m → ∞, γm {Bt (β0 )} = o(m2 ), so the reduced dependence approximation to the likelihood with km = cβ0 log m will provide
asymptotically valid inference for any cβ0 > b−1
β0 . The cost of computing this
approximation is O(m log m mcβ0 log 2 ) < O(mcβ0 log 2+2 ), polynomial in the size
of the model, but increasing rapidly as β0 approaches the critical value.
3.2.4

Extension to the more general case

If α 6= 0, no simple expression for Zr,c (α, β) is known. Instead, we investigate
(k)
(K)
the behaviour of m−1 δm (α, β) numerically, by using ℓ̃m for some large K as a
proxy for the exact log-likelihood. Care is needed to choose K sufficiently large
to ensure that the results are not sensitive to the choice of K.
(k)
Taking K = 16, a contour plot of this approximation to log{m−1 δm (0.1, 0.3)}
against log m and k is shown in Figure 2b, for k = 2, . . . , 12, m = 50, . . . 300.
The same plots with K = 15 (dashed) and K = 14 (dot-dash) are overlaid,
and the differences are barely visible at this scale. Given this stability, it seems
(k)
reasonable to assume a contour plot of the exact log{m−1 δm (0.1, 0.3)} would
look very similar to this. To obtain asymptotically valid inference, km should

9

20

10

-4

-3

15

-2

k

b−1
β

8

10

6

-1
0

5
4

1
2

0
0.10

0.20

0.30

0.40

β

4.0

3

4.5

5.0
log(m)

5.5

(k)

−1
(a) b−1
δm (0.1, 0.3)}
β , for a restricted Ising model with (b) Contour plot of log{m
against log m and k.
α = 0.

Figure 2: Inference for an Ising model with a reduced dependence approximation
to the likelihood.
be chosen so that this rescaled error shrinks with m, which seems to occur if km
grows at rate c(0.1, 0.3) log m, for c(0.1, 0.3) larger than about 1.5. This pattern
of behaviour is very similar to the α = 0 case. In both cases, reduced dependence approximations with an appropriate choice of k will give asymptotically
valid inference at cost polynomial in the size of the model, in contrast with the
exponential cost of computing the likelihood exactly.

4

Discussion

The results obtained here can also be applied to other approximations to the
likelihood and to other types of model. The conditions on the approximate
likelihood, such as showing that δn∞ {Bt (θ0 )} = op (rn ), may be difficult to verify
in practice, as the exact likelihood is assumed unavailable. If the approximation
to the likelihood is a truncation of a series expansion for the exact likelihood,
as is the case for the Laplace approximation, the conditions can be checked by
examining the contribution from higher-order terms in the expansion. In other
cases, it may only be possible to investigate the size of the errors numerically,
by using used a more accurate and expensive approximation to the likelihood
as a proxy for the exact likelihood, as described in Section 3.2.4 for the Ising
model example.
Many approximation methods have a tuning parameter, k say, where increasing k allows computation of a more accurate likelihood approximation at
increased cost. The reduced dependence approximation at level k described in
Section 3.2.1 is one example of this. In order for the approximate likelihood
inference to be close to the exact likelihood inference, k = kn should be allowed
to vary with n. Given an understanding of how the error in the score function
varies with k and n, the results obtained here could be applied to determine
how kn should scale with n. This has the potential to allow the construction
of approximate likelihoods which match the inference with the exact likelihood
10

closely for all n, but which scale well to large data sizes.

Acknowledgements
I am grateful to David Firth, Nancy Reid, Cristiano Varin and two referees
for helpful comments which have greatly improved this paper. This work was
supported by the Engineering and Physical Sciences Research Council through
Intractable Likelihood: New Challenges from Modern Applications [grant number EP/K014463/1].

Appendix: proofs of results
Proof of Theorem 1
Proof. We apply Theorem 5.9 of van der Vaart (2000), with Ψn (θ) = rn−1 ũn (θ)
and Ψ(θ) = ū(θ). Then
supkΨn (θ) − Ψ(θ)k = supkūn (θ) + rn−1 ∇θ ǫn (θ) − ū(θ)k

θ∈Θ

θ∈Θ

≤ supkūn (θ) − ū(θ)k + rn−1 δn
θ∈Θ

= op (1)
The conditions of Theorem 5.9 of van der Vaart (2000) are met, so since Ψn (θ̃n ) =
0, we have that θ̃n → θ0 in probability, as n → ∞.

Proof of Lemma 1
Proof. Taking a Taylor expansion of ūn about θ̂n ,
ūn (θ) = ūn (θ̂n ) − (θ − θ̂n )T J¯n (θn∗ ) = −(θ − θ̂n )T J¯n (θn∗ ),
¯n (θ) = rn−1 ũn (θ). Then
for some θn∗ between θ and θ̂n . Define ũ
ũ¯n (θ) = ūn (θ) + rn−1 ∇θ ǫn (θ) = −(θ − θ̂n )T J¯n (θn∗ ) + rn−1 ∇θ ǫn (θ),
¯n (θ̃n ) = 0 solves
so any θ̃n solving ũ
θ̃n − θ̂n = J¯n−1 (θn∗ )rn−1 ∇θ ǫn (θ̃n ),
for some θn∗ between θ̃n and θ̂n .
But θn∗ is a consistent estimator of θ, because θ̃n is by Theorem 1, so
¯
Jn (θn∗ ) → I(θ0 ) in probability. So
θ̃n − θ̂n = Op (rn−1 δn (θ̃n )).
Write An for the event {θ̃n ∈ Bt (θ0 )}. Then pr(An ) → 1 as n → ∞ since θ̃n
is consistent, and conditional on An ,


θ̃n − θ̂n = Op rn−1 δn∞ {Bt (θ0 )} = op (rn−1 an ).
11

For any ǫ > 0,
pr(kθ̃n − θ̂n k ≥ ǫan rn−1 )

C
= pr(kθ̃n − θ̂n k ≥ ǫan rn−1 |An )pr(An ) + pr(kθ̃n − θ̂n k ≥ ǫan rn−1 |AC
n )pr(An )

≤ pr(kθ̃n − θ̂n k ≥ ǫan rn−1 |An ) + pr(AC
n)

→0

as n → ∞.

Proof of Theorem 2
1/2

Proof. Applying Lemma 1 with an = rn ,
θ̃n − θ̂n = op (rn1/2 rn−1 ) = op (rn−1/2 ).
So
rn1/2 (θ̃n − θ0 ) = rn1/2 (θ̂n − θ0 ) + op (1) → N (0, I(θ0 )−1 )
in distribution, as required.

Proof of Theorem 3
Proof. We have
(Λ̃n − Λn )/2 = {ℓ̃n (θ̃n ) − ℓ̃n (θ̃nR )} − {ℓn (θ̂n ) − ℓn (θ̂nR )}

= {ℓ̃n (θ̃n ) − ℓ̃n (θ̂n )} + {ℓ̃n (θ̂n ) − ℓn (θ̂n )} + {ℓ̃n (θ̂nR ) − ℓ̃n (θ̃nR )} + {ℓn (θ̂nR ) − ℓ̃n (θ̂nR )}.

Taking a Taylor expansion of the first term,
ℓ̃n (θ̂n ) − ℓ̃n (θ̃n ) = (θ̂n − θ̃n )T ũn (θ̃n ) + (θ̂n − θ̃n )T J˜n (θ̄n )(θ̂n − θ̃n )
= (θ̂n − θ̃n )T J˜n (θ̄n )(θ̂n − θ̃n )
for some θ̄n between θ̂n and θ̃n , so ℓ̃n (θ̂n ) − ℓ̃n (θ̃n ) = op (1), since kθ̂n − θ̃n k =
−1/2
op (rn ) and J˜n (θ̄n ) = Jn (θ̄n ) + op (rn ) = Op (rn ). Similarly, ℓ̃n (θ̂nR ) − ℓ̃n (θ̃nR ) =
op (1), so
(Λ̃n − Λn )/2 = ǫn (θ̂n ) − ǫn (θ̂nR ) + op (1)

= (θ̂n − θ̂nR )T ∇θ ǫn (θn∗ ) + op (1)
−1/2

for some θn∗ between θ̂n and θ̂nR . But θ̂n − θ̂nR = Op (rn ), and, for sufficiently
large n, k∇θ ǫn (θn∗ )k ≤ δn∞ {Bt (θ0 )}, so Λ̃n − Λn = op (1), as required.

Results needed to prove Theorem 4
In order to prove Theorem 4, it is helpful to first consider properties of inference
with the penalized loglikelihood ℓπn (θ) = ℓ(θ) + log π(θ) with log-prior penalty.
We write θ̂nπ for the corresponding penalized likelihood estimator, which is the
posterior mode. Similarly, write Aπ for the version of the quantity A computed
12

with ℓπn (.) in place of ℓn (.), and Ãπ for the approximate version. Since ǫπn (θ) =
ǫn (θ), all the error terms remain unchanged.
Under the regularity conditions assumed on the model, the penalized likelihood estimator θ̂nπ will be consistent, and for any bn = op (rn ), the posterior
probability that θ ∈ Bb−1/2 (θ0 ) will tend to one as n → ∞.
n
To prove Theorem 4, we will use the following lemma, which says that the
error in the penalized log-likelihood ratio may be approximated in terms of the
error in the score function.
Lemma 2. Suppose that δn∞ = op (rn ), and that there exists t > 0 such that
1/2
δn∞ {Bt (θ0 )} = op (rn ). Then θ̂nπ and θ̃nπ are consistent estimators of θ0 , and
−1/2
θ̂nπ − θ̃nπ = op (rn ).
If γn∞ {Bt (θ0 )} = op (bn ), for some bn = o(rn ), then
sup
π
−1/2 (θ̂n )
bn

θ∈B

 

 π π
{ℓ̃n (θ̃n )− ℓ̃πn (θ)}−{ℓπn (θ̂nπ )−ℓπn (θ)} − (θ̂nπ −θ)T ∇θ ǫn (θ) = op (1).
−1/2

Proof. That θ̂nπ − θ̃nπ = op (rn
Theorem 2. Write

) follows by a similar argument to the proof of

Cn (θ) = {ℓ̃πn (θ̃nπ ) − ℓ̃πn (θ)} − {ℓπn (θ̂nπ ) − ℓπn (θ)}

= {ℓ̃πn (θ̂nπ ) − ℓπn (θ̂nπ )} − {ℓ̃πn (θ) − ℓπn (θ)} − {ℓ̃πn (θ̂nπ ) − ℓ̃πn (θ̃nπ )}

= ǫn (θ̂nπ ) − ǫn (θ) − {ℓ̃πn (θ̂nπ ) − ℓ̃πn (θ̃nπ )}.
Then

ℓ̃n (θ̂nπ ) − ℓ̃n (θ̃nπ ) = (θ̂nπ − θ̃nπ )T J˜nπ (θ̄n )(θ̂nπ − θ̃nπ )

for some θ̄n between θ̂nπ and θ̃nπ , and
ǫn (θ̂nπ ) − ǫn (θ) = (θ̂nπ − θ)T ∇θ ǫn (θ) + (θ̂nπ − θ)T ∇Tθ ∇θ ǫn (θn∗ (θ))(θ̂nπ − θ)
for some θn∗ (θ) between θ̂nπ and θ. Write
Dn (θ) = Cn (θ) − (θ̂nπ − θ)T ∇θ ǫn (θ)

= (θ̂nπ − θ)T ∇Tθ ∇θ ǫn (θn∗ (θ))(θ̂nπ − θ) − (θ̂nπ − θ̃nπ )T J˜nπ (θ̄n )(θ̂nπ − θ̃nπ ).

Then supθ∈B

π
−1/2 (θ̂n )
bn

sup
π
−1/2 (θ̂n )
bn

θ∈B

≤

n

sup
π)
θ∈B −1/2 (θ̂n
bn

Dn (θ) may be expressed as

(θ̂nπ − θ)T ∇Tθ ∇θ ǫn (θn∗ (θ))(θ̂nπ − θ) − (θ̂nπ − θ̃nπ )T J˜nπ (θ̄n )(θ̂nπ − θ̃nπ )
n

(θ̂nπ − θ)T ∇Tθ ∇θ ǫn (θn∗ (θ))(θ̂nπ − θ)

o

o

+ (θ̂nπ − θ̃nπ )T J˜nπ (θ̄n )(θ̂nπ − θ̃nπ )

∞
≤ b−1
n δn {Bt (θ0 )} + op (1)

for n sufficiently large. But δn∞ {Bt (θ0 )} = op (bn ), so supθ∈B

π
−1/2 (θ̂n )
bn

op (1), as required.

13

Dn (θ) =

Proof of Theorem 4
Proof. The normalized exact and approximate posterior distributions
are π(θ|y) =
R
Ln (θ)π(θ)/Zn and π̃(θ|y) = L̃n (θ)π(θ)/Z̃n , where Zn = Ln (θ)π(θ)dθ and
R
Z̃n = L̃n (θ)π(θ)dθ.
First, we find a Laplace approximation
Ẑn = (2π)−p/2 |Jnπ (θ̂nπ )|−1/2 Lπn (θ̂nπ )
to Zn . Then log Zn − log Ẑn = op (1), because Zn is a p-dimensional integral,
where p remains fixed as n → ∞. Similarly, Z̃n may be approximated by using
a Laplace approximation. Then
log Zn − log Z̃n = log Ẑn − log Z̃ˆn + op (1)
= (log |Jnπ (θ̂nπ )| − log |J˜nπ (θ̃nπ )|)/2 + ℓπn (θ̂nπ ) − ℓ̃πn (θ̃nπ ) + op (1).
Since γn∞ {Bt (θ0 )} = op (rn ), both rn−1 Jnπ (θ̂nπ ) and rn−1 J˜nπ (θ̂nπ ) converge towards
I(θ0 ), so
log |Jnπ (θ̂nπ )| − log |J˜nπ (θ̃nπ )| = log |rn−1 Jnπ (θ̂nπ )| − log |rn−1 J˜nπ (θ̃nπ )|
= log |I(θ0 ) + op (1)| − log |I(θ0 ) + op (1)|
= op (1).

So
log Zn − log Z̃n = ℓπn (θ̂nπ ) − ℓ̃πn (θ̃nπ ) + op (1).
−1/2

Since δn∞ {Bt (θ0 )} = op (rn ) and γn∞ {Bt (θ0 )} = op (rn ), we may choose bn =
−1/2
o(rn ) such that δn∞ {Bt (θ0 )} = op (bn ) and γn∞ {Bt (θ0 )} = op (bn ). Writing
Sn = Bb−1/2 (θ̂nπ ),
n

sup
θ∈Sn

n

o
log π̃(θ|y) − log π(θ|y) − (θ̂nπ − θ)T ∇θ ǫn (θ)
o
n
= sup ℓ̃πn (θ) − log Z̃n − ℓπn (θ) + log Zn − (θ̂nπ − θ)T ∇θ ǫn (θ)
θ∈Sn
o
n
≤ sup ℓ̃πn (θ) − ℓ̃πn (θ̃nπ ) − ℓπn (θ) + ℓ̂πn (θ̂nπ ) − (θ̂nπ − θ)T ∇θ ǫn (θ) + op (1)
θ∈Sn

= op (1),

14

by Lemma 2. Then
2dT V {π̃(θ|y), π(θ|y)} =

Z

=

ZΘ

≤

Z

π̃(θ|y) − π(θ|y) dθ

Sn

Sn

π̃(θ|y) − π(θ|y) dθ +

Z

π̃(θ|y) − π(θ|y) dθ +

Z

C
Sn

π̃(θ|y) − π(θ|y) dθ
2 π(θ|y) dθ +

C
Sn

Z

2 π̃(θ|y) dθ

C
Sn

π̃(θ|y)
− 1 π(θ|y)dθ + op (1)
π(θ|y)
Sn
π̃(θ|y)
≤ sup
− 1 + op (1)
θ∈Sn π(θ|y)
=

Z

= sup exp{log π̃(θ|y) − log π(θ|y)} − 1 + op (1)
θ∈Sn

≤ exp



sup | log π̃(θ|y) − log π(θ|y)| − 1 + op (1).

θ∈Sn

But
sup | log π̃(θ|y) − log π(θ|y)| ≤ sup |(θ̂nπ − θ)T ∇θ ǫn (θ)| + op (1)

θ∈Sn

θ∈Sn

= bn−1/2 δn∞ {Bt (θ0 )} + op (1)
= op (1).
So dT V {π̃(θ|y), π(θ|y)} = op (1), as claimed.

References
Douglas Bates, Martin Mächler, Ben Bolker, and Steve Walker. Fitting linear
mixed-effects models using lme4. Journal of Statistical Software, 67(1):1–48,
2015. doi: 10.18637/jss.v067.i01.
Randal Douc, Éric Moulines, and Tobias Rydén. Asymptotic properties of the
maximum likelihood estimator in autoregressive models with Markov regime.
The Annals of Statistics, 32(5):2254–2304, 2004.
Nial Friel, AN Pettitt, Robert Reeves, and Ernst Wit. Bayesian inference in
hidden markov random fields for binary data defined on large lattices. Journal
of Computational and Graphical Statistics, 18(2):243–261, 2009.
V. P. Godambe. An optimum property of regular maximum likelihood estimation. Ann. Math. Statist., 31(4):1208–1211, 12 1960. doi: 10.1214/aoms/
1177705693. URL http://dx.doi.org/10.1214/aoms/1177705693.
Michael I. Jordan.
Graphical models.
Statist. Sci.,
140–155, 02 2004.
doi:
10.1214/088342304000000026.
http://dx.doi.org/10.1214/088342304000000026.

19(1):
URL

Bruria Kaufman. Crystal statistics. II. Partition function evaluated by spinor
analysis. Physical Review, 76(8):1232, 1949.
15

Bruce G. Lindsay. Composite likelihood methods. Contemporary Mathematics,
80:221–239, 1988.
Francesco Pauli, Walter Racugno, and Laura Ventura. Bayesian composite
marginal likelihoods. Statistica Sinica, pages 149–164, 2011.
José C. Pinheiro and Douglas M. Bates. Approximations to the log-likelihood
function in the nonlinear mixed-effects model.
Journal of Computational and Graphical Statistics, 4(1):12–35, 1995. ISSN 10618600. URL
http://www.jstor.org/stable/1390625.
Mathieu Ribatet, Daniel Cooley, and Anthony C Davison. Bayesian inference
from composite likelihoods, with an application to spatial extremes. Statistica
Sinica, pages 813–845, 2012.
Håvard Rue, Sara Martino, and Nicolas Chopin. Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 71(2):319–392, 2009.
Håkon Tjelmeland and Haakon Michael Austad.
Exact and approximate recursive calculations for binary markov random fields defined on graphs.
Journal of Computational and Graphical Statistics, 21(3):758–780, 2012. doi: 10.1080/10618600.2012.632236. URL
http://dx.doi.org/10.1080/10618600.2012.632236.
Aad W van der Vaart. Asymptotic Statistics. Cambridge Series in Statistical
and Probabilistic Mathematics. Cambridge University Press, 2000.
Cristiano Varin, Nancy Reid, and David Firth. An overview of composite likelihood methods. Statistica Sinica, 21(1):5–42, 2011.

16

Supplementary material for
On asymptotic validity of naive inference with an
approximate likelihood
Helen Ogden
University of Southampton, UK

1

Finding δ(i) (θ) in Example 3.1

Recall Yi is the number of successes out of m = mn trials on item i. We study
how δ(i) (θ) = k(d/dθ)ǫ(i) (θ)k varies with m.
Write


f (b; yi ) = −yi log logit−1 (b) + (m − yi ) log 1 − logit−1 (b)
and

g(b; θ, yi ) = f (b; yi ) − log φ(b; 0, θ),
so that
L(i) (θ) =

Z

∞

exp{−g(b; θ, yi )}db.

−∞

In the following, we drop the data yi from the notation for convenience. Write
b̂(θ) for the maximizer of g(., θ), and
ĝr (θ) =

∂rg
(b̂(θ); θ).
∂bk

By equation (4) of Shun and McCullagh (1995), the error in the Laplace
approximation to the log-likelihood ℓ(i) (θ) is
ǫ(i) (θ) =

∞
X
1 X
n2 (P )(−1)v ĝ|p1 | (θ) . . . ĝ|pv | (θ){ĝ2 (θ)}−l ,
2l!
l=1

(S1)

P ∈P2l

where P = p1 | . . . |pv is a partition of 2l indices into v blocks of size 3 or more,
and n2 (P ) is the number of partitions Q of 2l indices into l blocks of size 2, such
that Q is complementary to P .
Write hP (θ) = ĝ|p1 | (θ) . . . ĝ|pv | (θ){ĝ2 (θ)}−l . Then hP (θ) = Op (mv−l ), since
ĝr (θ) = Op (m) for each r.
Differentiating (S1) gives
∞
X
d
1 X
d
ǫ(i) (θ) =
n2 (P )(−1)v hP (θ),
dθ
2l!
dθ
l=1

P ∈P2l

1

(S2)

and
v
v h
i
Y
Y
X
d
′
−l
lĝ|pj | (θ)ĝ2′ (θ){ĝ2 (θ)}−(l+1) ,
ĝ|p
(θ)
ĝ
(θ){ĝ
(θ)}
−
hP (θ) =
2
|pj |
i|
dθ
j=1
i=1
j6=i

(S3)
where ĝr′ (θ) = (d/dθ)ĝr (θ).
For each r, we have
o
d n (r)
g (b̂(θ); θ)
dθ
∂g (r)
db̂(θ)
=
(b̂(θ); θ) +
ĝr+1 (θ).
∂θ
dθ

ĝr′ (θ) =

(S4)

We now study the size of each of the terms in (S4). We have
∂
∂g (r)
(b; θ) =
∂θ
∂θ



∂r
− r log φ(b; 0, θ) = Op (1),
∂b

(S5)

and
ĝr+1 (θ) = Op (m).

(S6)

For each θ, b̂(θ) satisfies g1 (b̂(θ); θ) = 0. Differentiating this with respect to θ,
db̂(θ)
∂g1
g2 (b̂(θ); θ) +
(b̂(θ); θ) = 0.
dθ
∂θ
But

∂g1
(b; θ) = −2bθ−3 ,
∂θ

so
db̂(θ)
= 2b̂(θ)θ−3 {ĝ2 (θ)}−1 = Op (m−1 ).
(S7)
dθ
Substituting (S5), (S6) and (S7) into (S4) gives that ĝr′ (θ) = Op (1) for each
r. From (S3) we then have
d
hP (θ) = Op (mv−l−1 ).
dθ
The highest order terms in (S2) come from partitions with (l, v) = (2, 1) or
(3, 2), and so δ(i) (θ) = Op (m−2 ).

2

(k)

Finding δm (β) in Example 3.2

Kaufman (1949) provides an exact expression for the normalizing constant for
an Ising model on an n × m lattice, with α = 0 and periodic boundary, as
Zn×m (0, β) = {2 sinh(2β)}

nm/2

Ān,m (β)/2,

where
(2)
(3)
(4)
Ān,m (β) = A(1)
n,m (β) + An,m (β) + An,m (β) + An,m (β),

2

and
A(1)
n,m (β) =
A(3)
n,m (β) =

n
Y

q=0
n
Y

2 cosh (m a2q+1,n (β)/2) , A(2)
n,m (β) =

2 sinh (m a2q+1,n (β)/2) ,

q=0

A(4)
n,m (β) =

2 cosh (m a2q,n (β)/2) ,

q=0

where

n
Y

n
Y

2 sinh (m a2q,n (β)/2)

q=0


al,n (β) = cosh−1 cosh(2β)2 / sinh(2β) − cos(πl/n)

for l ≥ 1, and a0,n (β) = a0 (β) = 2β + log {tanh(β)}.
(k)
Using the approximation Zm×m (β) to Zm×m (β), the error in the log-likelihood
is
ǫ(k)
m (β) = (m − k + 1) log Āk,m (β) − (m − k) log Āk−1,m (β) − log Ām,m (β).
Differentiating this with respect to β,
d 
d 
d (k)
d 
log Āk,m (β) −(m−k)
log Āk−1,m (β) −
log Ām,m (β) ,
ǫ (β) = (m−k+1)
dβ m
dβ
dβ
dβ
(S8)
and
4
o
X
d 
d n
(i)
log A(i)
log Ān,m (β) =
n,m (β) rn,m (β),
dβ
dβ
i=1
(i)

(i)

where rn,m (β) = An,m (β)/Ān,m (β).
We have
n
o
X
d n
log A(1)
(β)
=
m/2
a′2q+1,n (β) tanh(ma2q+1,n (β)/2)
n,m
dβ
q=0
= m/2

n
X

a′2q+1,n (β) + O(m exp{−a0 (β)m})

q=0

as m → ∞, since tanh(x) = 1 + O(exp{−2x}) as x → ∞, and a2q+1,n (β) ≥
a0 (β) > 0.
n
o
(i)
d
log An,m (β) , and
Similar expressions may be obtained for the other dβ
combining these gives
d 
(o)
(e)
log Ān,m (β) = mSn(o) (β)rn,m
(β)+mSn(e) (β)rn,m
(β)+O(m exp{−a0 (β)m})
dβ
Pn
Pn
(1)
(o)
(e)
(o)
′
′
=
where Sn =
q=0 a2q,n (β), rn,m (β) = rn,m (β) +
q=0 a2q+1,n (β), Sn
(2)

(e)

(3)

(4)

rn,m (β) and rn,m (β) = rn,m (β) + rn,m (β). Define

f (x; β) = dβ {−1 + cβ − cos(x)}−1/2 {1 + cβ − cos(x)}−1/2
where dβ = 4 cosh(2β) − 2 cosh(2β) coth(2β)2 and cβ = cosh(2β)2 / sinh(2β).
(o)
(e)
Then a′j,n (β) = f (jπ/n; β), and n−1 Sn (β) and n−1 Sn (β) are both trapezium
R
2π
(o)
1
−1 (o)
rule approximations to I(β) = 2π
Sn (β) −
0 f (x; β)dx. Write Rn (β) = n
(e)
−1 (e)
I(β) and Rn (β) = n Sn (β) − I(β) for the error in each of these approximations to the integral.
3

(e)

(o)

Lemma 3. For each β < βc , Rn (β) = max{|Rn (β)|, |Rn (β)|} = O(exp{−bβ n}),
where bβ = 2 cosh−1 {−1 + cosh(2β)2 / sinh(2β)}.
Proof. We apply the results of Trefethen and Weideman (2014) to show exponentially fast convergence of these trapezium rule approximations to I(β).
These results depend on properties of the integrand f (z, β), considered as a
function of complex-valued z. There are a branch points of f (z, β) at a distance
aβ = cosh−1 {−1 + cosh(2β)2 / sinh(2β)} from the real axis, and the function is
analytic for −aβ < Im z < aβ , so by Theorem 3.2 of Trefethen and Weideman
(o)
(2014), |Rn (β)| = O(exp{−2aβ n}) = O(exp{−bβ n}).
(e)
(o)
The same argument holds with Rn (β) in place of Rn (β), so Rn (β) =
O(exp{−bβ n}), as required.
We now prove the main result.
(k)

Lemma 4. If k → ∞ as m → ∞, δm (β) = O(m2 k exp{−bβ k}) + o(1).
Proof. We have
d 
log Ān,m (β) = mnI(β) + mntn,m (β) + O(m exp{−a0 (β)m})
dβ
(o)

(o)

(e)

(e)

where tn,m (β) = Rn (β)rn,m (β) + Rn (β)rn,m (β).
Substituting this into (S8), the contributions from the mnI(β) terms cancel,
and the combined remainder terms are always o(1), since m2 exp{−a0 (β)m} =
o(1). We are left with
d (k)
ǫ (β) = (m − k + 1)mtk,m (β) − (m − k)mtk−1,m (β) − mtm,m (β) + o(1).
dβ m
Then
d (k)
ǫ (β)
dβ m
≤ (m − k + 1)mk|tk,m (β)| + (m − k)m(k − 1)|tk−1,m (β)| + m2 |tm,m (β)| + o(1)

(k)
δm
(β) =

≤ (m − k + 1)mkRk (β) + (m − k)m(k − 1)Rk−1 (β) + m2 Rm (β) + o(1)
(o)

(o)

(e)

(e)

since |tn,m (β)| ≤ |Rn (β)|rn,m (β) + |Rn (β)|rn,m (β) ≤ Rn (β)
= O(m2 k exp{−bβ k}) + o(1)
by Lemma 3, as required.

References
Bruria Kaufman. Crystal statistics. II. Partition function evaluated by spinor
analysis. Physical Review, 76(8):1232, 1949.
Zhenming Shun and Peter McCullagh. Laplace approximation of high dimensional integrals. Journal of the Royal Statistical Society. Series B (Methodological), 57(4):749–760, 1995.
Lloyd N Trefethen and JAC Weideman. The exponentially convergent trapezoidal rule. SIAM Review, 56(3):385–458, 2014.
4

