# Bayesian leave-one-out cross-validation for large data

Måns Magnusson  Michael Riis Andersen  Johan Jonasson  Aki Vehtari

###### Abstract

Model inference, such as model comparison, model checking, and model
selection, is an important part of model development. Leave-one-out cross-
validation (LOO) is a general approach for assessing the generalizability of a
model, but unfortunately, LOO does not scale well to large datasets. We
propose a combination of using approximate inference techniques and
probability-proportional-to-size-sampling (PPS) for fast LOO model evaluation
for large datasets. We provide both theoretical and empirical results showing
good properties for large data.

Approximate inference, Leave-one-out cross-validation, Model inference,
Subsampling

##  1 Introduction

Model inference, such as model comparison, checking, and selection, is an
integral part of developing new models. From a Bayesian decision-theoretic
point of view (see Vehtari & Ojanen (2012) for an extended discussion) we want
to make a choice $a\in\mathcal{A}$, in our case a model $p_{M}$, that maximize
our _expected utility_ for a utility function $u(a,\cdot)$ as

$\bar{u}(a)=\int u(a,\tilde{y_{i}})p_{t}(\tilde{y}_{i})d\tilde{y}_{i}\,,$

where $p_{t}(\tilde{y}_{i})$ is the true probability distribution generating
observation $\tilde{y}_{i}$.

A common scenario is to study how well a model _generalizes_ to unseen data
(Box, 1976; Vehtari & Ojanen, 2012; Vehtari et al., 2017). A popular utility
function $u$ with good theoretical properties for probabilistic models is the
log score function (Bernardo, 1979; Robert, 1996; Vehtari & Ojanen, 2012). The
log score function give rise to using the _expected log predictive density_
(elpd) for model inference, defined as

$\overline{\text{elpd}}_{M}=\int\log
p_{M}(\tilde{y}_{i}|y)p_{t}(\tilde{y}_{i})d\tilde{y}_{i}\,,$

where $\log p_{M}(\tilde{y}_{i}|y)$ is the log predictive density for a new
observation for the model $M$.

Leave-one-out cross-validation (LOO-CV) is one approach to estimate the elpd
for a given model, and is the method of focus in this paper (Bernardo & Smith,
1994; Vehtari & Ojanen, 2012; Vehtari et al., 2017). Using LOO-CV we can treat
our observations as pseudo-Monte Carlo samples from $p_{t}(\tilde{y}_{i})$ and
estimate the $\overline{\text{elpd}}_{\text{loo}}$ as

$\displaystyle\overline{\text{elpd}}_{\text{loo}}$
$\displaystyle=\frac{1}{n}\sum^{n}_{i=1}\log p_{M}(y_{i}|y_{-i})$ (1)
$\displaystyle=\frac{1}{n}\sum^{n}_{i=1}\log\int
p_{M}(y_{i}|\theta)p_{M}(\theta|y_{-i})d\theta$
$\displaystyle=\frac{1}{n}\,\text{elpd}_{\text{loo}}\,,$

where $n$ is the number of observations (that may be very large),
$p_{M}(y_{i}|\theta)$ is the likelihood, and $p_{M}(\theta|y_{-i})$ is the
posterior for $\theta$ where we hold out observation $i$. This will henceforth
be called the LOO posterior and $p_{M}(\theta|y)$ will be referred to as the
full posterior. In this paper both $\overline{\text{elpd}}_{\text{loo}}$ and
$\text{elpd}_{\text{loo}}$ will be quantities of interest, depending on the
situation.

Bayesian LOO-CV has many appealing theoretical properties compared to other
common model evaluation techniques. The popular $k$-fold cross-validation is,
in general, a biased estimator of $\text{elpd}_{M}$, since each model is only
trained using a subset of the full data (Vehtari & Ojanen, 2012). The LOO-CV
is, just as the Watanabe-Akaike Information criteria (WAIC), a consistent
estimate of the true $\text{elpd}_{M}$ for regular and singular models
(Watanabe, 2010). A model is regular if the map taking the parameters to the
probability distribution is one-to-one and the Fisher information is positive-
definitive. If a model is not regular, then the model is singular (Watanabe,
2010). Since many models, such as neural networks, normal mixture models,
hidden Markov models, and topic models, are singular, we need consistent
methods to estimate the elpd for singular models (Watanabe, 2010). Although
the WAIC and LOO-CV have the same asymptotic properties, recent research has
shown that the LOO-CV is more robust than WAIC in the finite data domain
(Vehtari et al., 2017).

In addition to the theoretical properties, the LOO-CV also gives an intuitive
framework for evaluating models where the user easily can use different
utility functions of interest as well as easily taking hierarchical data
structures into account by using leave-one-group-out or leave-one-cluster out
cross-validation (see Merkle et al. (2018) for a discussion). Taken together,
LOO-CV has many very good properties, both empirical and theoretical. In this
paper, we will focus on LOO-CV as a way of evaluating models.

Modern probabilistic machine learning techniques need to scale to massive
data. In a data-rich regimes, we often want complex models, such as
hierarchical and non-linear models. Model comparison and model evaluation are
important for model development, but little focus has been put into finding
ways of scaling LOO-CV to larger data. The main problem is that a straight-
forward implementation means that $n$ models need to be estimated. Even if
this problem is solved, for example using importance sampling (see below), we
still have two problems.

First (1), many posterior approximation techniques, such as Markov Chain Monte
Carlo (MCMC), does not generally scale to large $n$ or is computationally very
costly. Second (2), computing $\text{elpd}_{\text{loo}}$ still needs to be
computed over $n$ observations. If it is costly to estimate individual
contributions (i.e. $\log p_{M}(y_{i}|y_{-i})$), computing the total
$\text{elpd}_{\text{loo}}$ may be very costly for very large models.

###  1.1 Pareto-smoothed importance sampling

If we would implement LOO-CV naively, inference needs to be repeated $n$ times
for each model. Gelfand (1996) propose the use of importance sampling to solve
this problem. The idea is to estimate $p_{M}(y_{i}|y_{-i})$ in Eq. (1) using
the importance sampling approximation

$\displaystyle\log\hat{p}(y_{i}|y_{-i})=\log\left(\frac{\frac{1}{S}\sum_{s=1}^{S}p_{M}(y_{i}|\theta_{s})r(\theta_{s})}{\frac{1}{S}\sum_{s=1}^{S}r(\theta_{s})}\right),$
(2)

where $\theta_{s}$ are $s\in 1,...,S$ draws from the full posterior
$p(\theta|y)$, and

$\displaystyle r(\theta_{s})$
$\displaystyle=\frac{p_{M}(\theta_{s}|y_{-i})}{p_{M}(\theta_{s}|y)}$
$\displaystyle\propto\frac{1}{p_{M}(y_{i}|\theta_{s})}\,,$

where the last step is a well-known result of Gelfand (1996). The ratios
$r(\theta_{s})$ can be unstable due to a long right tail, but this can be
resolved using Pareto-smoothed importance sampling (PSIS) (Vehtari et al.,
2015). Using PSIS we fit a generalized Pareto distribution to the largest
weights $r(\theta_{s})$ and replace the largest importance sample ratios with
order statistics from the estimated generalized Pareto distribution,
decreasing the variance by introducing a small bias. PSIS also has the benefit
that we can use the estimated shape parameter $\hat{k}$ from the generalized
Pareto distribution to determine the reliability of the estimate. For data-
points with $\hat{k}>0.7$ the estimates of $\log p(y_{i}|y_{-i})$ _can_ be
unreliable and hence $\hat{k}$ can be used as a diagnostic (Vehtari et al.,
2017).

However, PSIS-LOO has the same scaling problem as LOO-CV in general since it
requires (1) samples from the true posterior (e.g. using MCMC) and (2) the
estimation of the $\text{elpd}_{\text{loo}}$ contributions from all
observations (Gelfand, 1996; Vehtari et al., 2017). Both of these requirements
can be costly in a data-rich regime and are the main problems we address in
this paper.

###  1.2 Contributions and limitations

In this paper, we focus on the problems of LOO-CV for large datasets and our
contributions are three-fold. First, we extend the method of Gelfand (1996) to
posterior approximations by including a correction term to the importance
sampling weights. Second, we propose sampling individual
$\text{elpd}_{\text{loo}}$ components with probability-proportional-to-size
sampling (PPS) to estimate $\text{elpd}_{\text{loo}}$. Third, we show
theoretically that these contributions have very favorable asymptotic
properties as $n\rightarrow\infty$. We show that the proposed estimator for
$\text{elpd}_{\text{loo}}$ is consistent for any consistent posterior
approximation $q$ (such as Laplace approximations, mean-field, and full-rank
variational inference posterior approximations). We also show that the
variance due to subsampling will decrease as the number of observations $n$
grows. In the limit, and given the assumptions in Section 2.3, we only need
one subsampled observation, and one draw from the full posterior, to estimate
$\overline{\text{elpd}}_{\text{loo}}$ with zero variance. Taken together this
introduces a new, fast, and theoretically motivated approach to model
evaluation for large datasets.

The limitations of our approach are the same as using PSIS-LOO (Vehtari et
al., 2017) as well as the requirement that the approximate posterior needs to
be sufficiently close to the true posterior (see Yao et al., 2018, for a
discussion).

##  2 Bayesian leave-one-out cross-validation for large data sets

Leave-one-out cross-validation (LOO-CV) has very good theoretical and
practical properties. This makes it relevant to develop tools to scale LOO-CV.
We solve this problem using scalable posterior approximations, such as Laplace
and variational approximations and using probability-proportional-to-log-
predictive-density subsampling inspired by Hansen & Hurwitz (1943).

###  2.1 Estimating the elpd using posterior approximations

Laplace and variational posterior approximations are attractive for fast model
comparisons due to their computational scalability. Laplace approximation
approximate the posterior distribution with multivariate normal distribution
$q_{Lap}(\theta|y)$ with the mean being the mode of the posterior and the
covariance the inverse Hessian at the mode (Azevedo-Filho & Shachter, 1994).

In variational inference, we minimize the Kullback-Leibler (KL) divergence
between an approximate family $\mathcal{Q}$ of densities and the true
posterior $p(\theta|y)$ (Jordan et al., 1999; Blei et al., 2017). Hence we
find the approximation $q$ that is closest to the true posterior in a KL
divergence sense. Here we let $\mathcal{Q}$ be a family of multivariate normal
distributions with a diagonal covariance structure (mean-field) or a full
covariance structure (full-rank). Hence we will work with a mean-field
variational approximation $q_{MF}(\theta|y)$ and a full-rank variational
approximation $q_{FR}(\theta|y)$.

Although, all these posterior approximations,
$q_{Lap}(\theta|y),q_{MF}(\theta|y)$, and $q_{FR}(\theta|y)$, will, in
general, be different than the true posterior distribution, we can use them as
a proposal distribution in an importance sampling scheme. In this scheme we
use a posterior approximation $q_{M}(\theta|y)$ for a model $M$ as the
proposal distribution and $p_{M}(\theta|y_{-i})$, the LOO posterior, as our
target distribution. The expectation of interest is the same as in the
standard PSIS-LOO given by Eq. (2), but we also propose to correct for the
posterior approximation error. Hence we change $r(\theta)$ to

$\displaystyle r(\theta_{s})$
$\displaystyle=\frac{p_{M}(\theta_{s}|y_{-i})}{q_{M}(\theta_{s}|y)}$
$\displaystyle=\frac{p_{M}(\theta_{s}|y_{-i})}{p_{M}(\theta_{s}|y)}\frac{p_{M}(\theta_{s}|y)}{q_{M}(\theta_{s}|y)}$
(3)
$\displaystyle\propto\frac{1}{p_{M}(y_{i}|\theta_{s})}\frac{p_{M}(\theta_{s}|y)}{q_{M}(\theta_{s}|y)}\,.$

The two-factor approach is needed to take the posterior approximation into
account. The factorization in Eq. (2.1) shows that the importance correction
contains two parts, the correction from the full posterior to the LOO
posterior and the correction from the full approximate distribution to the
full posterior. Both components often have lighter tailed proposal
distribution than the corresponding target distribution which can increase the
variance of the importance sampling estimate (Geweke, 1989; Gelfand, 1996).

Pareto-smoothed importance sampling can be used to both stabilize the weights
in estimating the contributions to the $\text{elpd}_{\text{loo}}$ and in
evaluating variational inference approximations using $\hat{k}$ as a
diagnostic (Vehtari et al., 2015; Yao et al., 2018). Hence we use PSIS to
stabilize the weights with the additional benefit that we can use $\hat{k}$,
the shape parameter in the generalized Pareto distribution, to diagnose how
well the approximation is working (Vehtari et al., 2015).

###  2.2 Probability-proportional-to-size subsampling and Hansen-Hurwitz
estimation

Using PSIS we can estimate each $\log\hat{p}(y_{i}|y_{-i})$ term and sum them
to estimate $\text{elpd}_{\text{loo}}$. Estimating every
$\log\hat{p}(y_{i}|y_{-i})$ can be costly, especially as $n$ grows. In some
situations using PSIS-LOO, estimating $\text{elpd}_{\text{loo}}$ can take even
longer than computing the full posterior once, due to the computational burden
of computing $\log\hat{p}(y_{i}|y_{-i})$, estimating $\hat{k}$ and using the
generalized Pareto distribution to stabalize the weights for each individual
observation. To handle this problem we suggest using a sample of the
$\text{elpd}_{\text{loo}}$ components to estimate $\text{elpd}_{\text{loo}}$.

Estimating totals, such as $\text{elpd}_{\text{loo}}$, has a long tradition in
sampling theory (see Cochran, 1977). If we have auxiliary variables that are a
good approximation of our variable of interest, we can use a probability-
proportional-to-size (PPS) sampling scheme to reduce the sampling variance in
the estimate of $\text{elpd}_{\text{loo}}$ using the _unbiased_ Hansen-Hurwitz
(HH) estimator (Hansen & Hurwitz, 1943). When evaluating models, we can often
easily compute $\log p_{M}(y_{i}|y)$, the full posterior log predictive
density, for all observations. We then sample $m<n$ observations proportional
to $\tilde{\pi}_{i}\propto\pi_{i}=-\log p_{M}(y_{i}|y)=-\log\int
p_{M}(y_{i}|\theta)p_{M}(\theta|y)d\theta$. We here assume that all $\log
p_{M}(y_{i}|y)<0$, but this assumption is only for convenience.

In the case of regular models and large $n$, we can also approximate $\log
p_{M}(y_{i}|y)\approx\log p_{M}(y_{i}|\hat{\theta})$ where $\hat{\theta}$ can
be a Laplace posterior mean estimate $\hat{\theta}_{q}$ or a VI mean estimate
${\mathbb{E}}_{\theta\sim q}[\theta]$. In the case of VI and Laplace
approximations, this further speeds up the computation of the
$\tilde{\pi}_{i}$s since we do not need to integrate over $\theta$ for all $n$
observations. Using a sampling with probability-proportional-to-size scheme,
the estimator for $\text{elpd}_{\text{loo}}$ can be formulated as

$\displaystyle\widehat{\overline{\text{elpd}}}_{\text{loo},q}=\frac{1}{n}\frac{1}{m}\sum^{m}_{i}\frac{1}{\tilde{\pi}_{i}}\log\hat{p}(y_{i}|y_{-i})\,,$
(4)

where $\tilde{\pi}_{i}$ is the probability of subsampling observation $i$,
$\log\hat{p}(y_{i}|y_{-i})$ is the (self-normalized) importance sampling
estimate of $\log p(y_{i}|y_{-i})$ given by Eq. (2) and (2.1), and $m$ is the
subsample size. The variance estimator can be expressed as (see Cochran, 1977,
Theorem 9A.2.).

$\displaystyle v(\widehat{\overline{\text{elpd}}}_{\text{loo},q})=$
$\displaystyle\frac{1}{n^{2}m(m-1)}\sum^{m}_{i=1}\left(\frac{\log\hat{p}(y_{i}|y_{-i})}{\tilde{\pi}_{i}}-n\widehat{\overline{\text{elpd}}}_{\text{loo}}\right)^{2}\,.$
(5)

The benefits of the HH estimator are many. First, if the probabilities are
proportional to the variable of interest ($\log\hat{p}(y_{i}|y_{-i})$ here),
the variance in Eq. (2.2) will go to zero, a property of use in the asymptotic
analysis in Section 2.3. Second, the estimator of
$\overline{\text{elpd}}_{\text{loo}}$ is not limited to posterior
approximation methods, but can also be used with MCMC (but without the
importance sampling correction factor). Third, PPS sampling has the benefit
that we can use Walker-Alias multinomial sampling (Walker, 1977). By building
an Alias table in $O(n)$ time we can then sample a new observation in $O(1)$
time. This means that can continue to sample observations until we have
sufficient precision in $\overline{\text{elpd}}_{\text{loo}}$ for our model
comparison purposes, independent of the number of observations $n$. Fourth,
the estimator is unbiased for all $\tilde{\pi}_{i}$. So by using $\log
p(y_{i}|\hat{\theta})$ instead of $\log p(y_{i}|y)$ we would expect a small
increase in variance since we would expect that for finite $n$, $\log
p(y_{i}|y)$ would be a better approximation of $\log p(y_{i}|y_{-i})$ than
$\log p(y_{i}|\hat{\theta})$, but at a greater computational cost.

To compare models, we are often also interested in the variance of
$\overline{\text{elpd}}_{\text{loo}}$, or for the dataset, henceforth called
$\sigma^{2}_{\text{loo}}$. To estimate $\sigma^{2}_{\text{loo}}$ we can use
the same observations as sampled previously, as

$\displaystyle\hat{\sigma}^{2}_{\text{loo}}$
$\displaystyle=\frac{1}{nm}\sum^{m}_{i}\frac{\hat{p}^{2}_{i}}{\tilde{\pi}_{i}}+$
(6)
$\displaystyle\frac{1}{n^{2}m(m-1)}\sum^{m}_{i}\left(\frac{\hat{p}_{i}}{\tilde{\pi}_{i}}-\frac{1}{m}\sum^{m}_{i}\frac{\hat{p}_{i}}{\tilde{\pi}_{i}}\right)^{2}-$
$\displaystyle\left(\frac{1}{nm}\sum^{m}_{i}\frac{\hat{p}_{i}}{\tilde{\pi}_{i}}\right)^{2}$

where $\hat{p}_{i}=\log\hat{p}(y_{i}|y_{-i})$. For a proof of unbiasedness of
the $\hat{\sigma}^{2}_{\text{loo}}$ estimator for $\sigma^{2}_{\text{loo}}$ in
Eq. (6), see the supplementary material. Also, note that here
$\sigma^{2}_{\text{loo}}=\frac{1}{n}\sum_{i}^{n}(\hat{p}^{2}_{i}-(\frac{1}{n}\sum_{i}^{n}\hat{p}^{2}_{i})^{2})$,
which in itself is not an unbiased estimate for the true
$\sigma^{2}_{\text{loo}}$ (Bengio & Grandvalet, 2004).

Although the variance estimator is unbiased, it is not as efficient as the
estimator of $\text{elpd}_{\text{loo}}$. This is part due to the fact that
$\tilde{\pi}_{i}$ is not proportional to $\hat{p}^{2}_{i}$ in the first line
in Eq. (6). This can be solved by sampling in two steps both proportional to
$\hat{p}_{i}$ and $\hat{p}^{2}_{i}$.

###  2.3 Asymptotic properties

For larger data sets the asymptotic properties of the method are crucial and
we derive asymptotic properties for the methods as follows. We consider a
generic Bayesian model; a sample $(y_{1},y_{2},\ldots,y_{n})$,
$y_{i}\in\mathcal{Y}\subseteq{\mathbb{R}}$, is drawn from a true density
$p_{t}=p(\cdot|\theta_{0})$ for some true parameter $\theta_{0}$. The
parameter $\theta_{0}$ is assumed to be drawn from a prior $p(\theta)$ on the
parameter space $\Theta$, which we assume to be an open and bounded subset of
${\mathbb{R}}^{d}$. A number of conditions are used. They are as follows.

  * (i)

the likelihood $p(y|\theta)$ satisfies that there is a function
$C:\mathcal{Y}\rightarrow{\mathbb{R}}_{+}$, such that ${\mathbb{E}}_{y\sim
p_{t}}[C(y)^{2}]<\infty$ and such that for all $\theta_{1}$ and $\theta_{2}$,
$|p(y|\theta_{1})-p(y|\theta_{2})|\leq
C(y)p(y|\theta_{2})\|\theta_{1}-\theta_{2}\|$.

  * (ii)

$p(y|\theta)>0$ for all $(y,\theta)\in\mathcal{Y}\times\Theta$,

  * (iii)

There is a constant $M<\infty$ such that $p(y|\theta)<M$ for all $(y,\theta)$,

  * (iv)

all assumptions needed in the Bernstein-von Mises (BvM) Theorem (Walker,
1969),

  * (v)

for all $\theta$, $\int_{\mathcal{Y}}(-\log p(y|\theta))p(y|\theta)dy<\infty$.

Of these assumptions, (i) and (iv) are the most restrictive. The assumption
that the parameter space is bounded is not very restrictive in practice since
we can approximate any proper prior arbitrarily well with a truncated
approximation.

######  Proposition 1.

Let the subsampling size $m$ and the number of posterior draws $S$ be fixed at
arbitrary integer numbers, let the sample size $n$ grow, assume that (i)-(v)
hold and let $q=q_{n}(\cdot|y)$ be any consistent approximate posterior. Write
$\hat{\theta}_{q}=\arg\max\\{q(\theta):\theta\in\Theta\\}$ and assume further
that $\hat{\theta}_{q}$ is a consistent estimator of $\theta_{0}$. Then

$|\widehat{\overline{\text{elpd}}}_{\text{loo}}(m,q)-\overline{\text{elpd}}_{loo}|\rightarrow
0$

in probability as $n\rightarrow\infty$ for any of the following choices of
$\pi_{i}$, $i=1,\ldots,n$.

  * (a)

$\pi_{i}=-\log p(y_{i}|y)$,

  * (b)

$\pi_{i}=-{\mathbb{E}}_{y}[\log p(y_{i}|y)]$,

  * (c)

$\pi_{i}=-{\mathbb{E}}_{\theta\sim q}[\log p(y_{i}|\theta)]$,

  * (d)

$\pi_{i}=-\log p(y_{i}|{\mathbb{E}}_{\theta\sim q}[\theta])$,

  * (e)

$\pi_{i}=-\log p(y_{i}|\hat{\theta}_{q})$.

###### Proof.

See the supplementary material. ∎

This proposition has three main points. First, the estimator of the
$\widehat{\overline{\text{elpd}}}_{\text{loo}}$ is consistent for any
consistent posterior approximation. In the limit, the mean-field variational
approximation will also estimate the true
$\overline{\text{elpd}}_{\text{loo}}$. Second, the estimator is also
consistent irrespective of the sub-sampling size $m$ and the number of draws,
$S$, from the posterior. This is a very good scaling characteristic. Third,
the estimator is consistent also if we approximate $\tilde{\pi}_{i}$ with
$\log p(y_{i}|\hat{\theta}_{q})$. This means that for larger data we can plug
in point estimates to quickly compute $\tilde{\pi}_{i}$ and still have the
consistency property.

The main limitations with Proposition 2 are that it is based on the
consistency of the posterior approximations and the proposition does only hold
for regular models for which $q$ are consistent. This is mainly due to the
fact that Laplace and VI are not, in general, consistent for singular models.

###  2.4 Computational complexity

In the large $n$ domain it is also of interest to study the computational
complexity of our approach. Assuming that the additional cost of computing
$p(y_{i}|y_{-i})$ compared to the point log predictive density (lpd) at
$\hat{\theta}$, $\log p(y_{i}|\hat{\theta})$, is $O(S)$, where $S$ is the
number of samples from the full posterior. Then the cost of computing the full
$\text{elpd}_{\text{loo}}$ is

$O(nS)\,.$

If we instead use our proposed method we would have the complexity

$O(n+mS)\,,$

where $m$ is the subsampling size. Using the proposed approach, we get an
_unbiased_ estimate of $\text{elpd}_{\text{loo}}$ together with the variance
$v(\widehat{\text{elpd}}_{\text{loo}})$ of that estimate, giving us
information on the precision of the method for a given $m$.

Finally, we could, for large $n$ just use the same lpd as an approximation
with complexity

$O(n)\,.$

This estimate is though _biased_ for all finite $n$, and we have no diagnostic
indicating how good or bad the approximation is.

This shows the large-scale characteristic of our proposed approach. By adding
a small cost ($mS$), we will have a good estimate of the true
$\text{elpd}_{\text{loo}}$ at the same cost as computing just the lpd. If
using the point lpd is a good approximation we would need less $m$. On the
other hand, if the point lpd would be a bad approximation, we would need a
larger $m$. The variance estimator in Eq. (2.2) would in these situations
serve as an indicator, with a higher variance estimate.

###  2.5 Method summary

We have presented a method for estimating the elpd efficiently using posterior
approximation and PPS subsampling. One of the attractive properties of the
method is that we can diagnose if the method is working. Using PSIS-LOO we can
diagnose the estimation of each individual $\log\hat{p}(y_{i}|y_{-i})$ as well
as the overall posterior approximation using the $\hat{k}$ diagnostic. Then
the variance of the HH estimator in Eq. (2.2) captures the effect of the
subsampling in the finite $n$ case. Our approach for large-scale LOO can be
summarized in the following steps.

  1. 1.

Estimate the models of interest using any consistent posterior approximation
technique.

  2. 2.

Compute the $\hat{k}$ diagnostic for the posterior to asses the general
overall posterior approximation. See (Yao et al., 2018) for an example for
variational inference.

  3. 3.

Compute $\tilde{\pi}_{i}\propto-\log p(y_{i}|y)$ for all $n$ observations. For
regular models this can be approximated with $\tilde{\pi}_{i}\propto-\log
p(y_{i}|\hat{\theta})$ for large data.

  4. 4.

Sample $m$ observations using PPS sampling and compute
$\log\hat{p}(y_{i}|y_{-i})$ using Eq. (2) and (2.1) for the sampled
observations. Use $\hat{k}$ to diagnose the estimation of each individual
$\log\hat{p}(y_{i}|y_{-i})$.

  5. 5.

Estimate $\widehat{\text{elpd}}_{\text{loo}}$,
$v(\widehat{\text{elpd}}_{\text{loo}})$, and $\hat{\sigma}^{2}_{\text{loo}}$
using Eq. (4), (2.2), and (6) to compare model predictive performance.

  6. 6.

Repeat step 3 and 4 until sufficient precision is reached.

The downside is that the $\hat{k}$ diagnostic can be too conservative for our
purpose. In the case of a correlated posterior and mean-field variational
inference, $\hat{k}$ may indicate a poor approximation even though the
estimation of $\text{elpd}_{\text{loo}}$ is still consistent and may work
well. In this situation, we would get a result indicating that all
$\log\hat{p}(y_{i}|y_{-i})$ are problematic, even though the estimation
actually work well, something we will see in the experiments.

##  3 Experiments

To study the characteristic of the proposed approach we study multiple models
and datasets. We use simulated datasets used to fit a Bayesian linear
regression model with $D$ variables and $N$ observations. The data is
generated such that so we get either a correlated (c) or an independent (i)
posterior for the regression parameters by construction. This will enable us
to study the effect of the mean-field assumptions in variational posterior
approximations. In addition, we use data from the radon example of Lin et al.
(1999) to show performance on a larger dataset with multiple models.

All posterior computations uses Stan 2.18 (Carpenter et al., 2017; Stan
Development Team, 2018) and all models used can be found in the supplementary
material. The methods has been implemented using the loo R package (Vehtari et
al., 2018) framework for Stan and is available as supplementary material. We
use mean-field and full-rank Automatic Differentiation Variational Inference
(ADVI) (Kucukelbir et al., 2017) and Laplace approximations as implemented in
Stan. ADVI automatically handles constrained variables and uses stochastic
variational inference.

###  3.1 Estimating $\text{elpd}_{\text{loo}}$ using posterior approximations

Data |  | ADVI(FR) | ADVI(MF) | Laplace | MCMC
---|---|---|---|---|---
LR(c) 100D | $\text{elpd}_{\text{loo}}$ | -14249 | -14267 | -14247 | -14247
|  $\hat{k}>0.7$ (%) | 100 | 100 | 0 | 0
LR(c) 10D | $\text{elpd}_{\text{loo}}$ | -14271 | -14271 | -14272 | -14272
|  $\hat{k}>0.7$ (%) | 0 | 100 | 0 | 0
LR(i) 100D | $\text{elpd}_{\text{loo}}$ | -14193 | -14239 | -14238 | -14239
|  $\hat{k}>0.7$ (%) | 100 | 0 | 0 | 0
LR(i) 10D | $\text{elpd}_{\text{loo}}$ | -14202 | -14202 | -14202 | -14203
|  $\hat{k}>0.7$ (%) | 0 | 0 | 0 | 0
Table 1: Estimation of $\text{elpd}_{\text{loo}}$ using posterior
approximations. For all models and posterior approximations,
$\sigma_{\text{LOO}}\approx 70$. No subsampling is used and MCMC is gold
standard.

Table 1 contains the estimatied values of $\text{elpd}_{\text{loo}}$ for
different posterior approximations. We used we used 100 000 iterations for
ADVI and 1000 warmup iterations and 2000 samples from 2 chains for the MCMC.
From the table, we can see that using PSIS-LOO and posterior approximations to
estimate $\text{elpd}_{\text{loo}}$ works well or diagnostic correctly
indicates the failure. As we would expect, the mean-field approximation for
the correlated posterior does not approximate the true posterior very well
when the posterior has correlated parameters and the $\hat{k}$ values are too
high for all observations. In spite of the high $\hat{k}$ values, the estimate
of the $\text{elpd}_{\text{loo}}$ is not very far from the (gold-standard)
MCMC estimate, showing the consistency result in Prop. 2 for mean-field
approximations - even when the true posterior covariance structure is not in
the variational family.

The second result is that the full-rank VI approximation has a poor fit for a
large number of parameters ($D=100$). This comes from that the full rank ADVI
needs to approximate the full posterior covariance structure (with 5 000
parameters) based on stochastic gradients. The increased perturbation in the
estimate of the covariance matrix has the effect of increasing the overall
$\hat{k}$, especially for larger dimensions indicating a less good
approximation of the posterior.

###  3.2 Subsampling using PPS sampling

Data | $m$ | Method | $\widehat{\text{elpd}}_{\text{loo}}$ | SE($\widehat{\text{elpd}}_{\text{loo}}$) | $\hat{\sigma}_{\text{loo}}$
---|---|---|---|---|---
LR(c) | - | True | -14247 | 0 | 71
100D | 10 | PPS(1) | -14236 | 11.8 | 23
|  | PPS(2) | -14244 | 13.6 | 67
|  | SRS | -14234 | 2197.8 | 70
| 100 | PPS(1) | -14248 | 4.3 | 67
|  | PPS(2) | -14249 | 4.5 | 73
|  | SRS | -13823 | 598.8 | 60
| 1000 | PPS(1) | -14245 | 1.5 | 69
|  | PPS(2) | -14248 | 1.5 | 71
|  | SRS | -14068 | 212.9 | 67
LR(c) | - | True | -14272 | 0 | 71
10D | 10 | PPS(1) | -14272 | 3.2 | 88
|  | PPS(2) | -14269 | 3.6 | 81
|  | SRS | -18096 | 3310.2 | 105
| 100 | PPS(1) | -14272 | 1.2 | 87
|  | PPS(2) | -14272 | 1.1 | 78
|  | SRS | -13921 | 669.9 | 67
| 1000 | PPS(1) | -14272 | 0.4 | 75
|  | PPS(2) | -14272 | 0.4 | 69
|  | SRS | -14266 | 223.3 | 71
Table 2: Effect of subsampling proportional to log predictive density. The
result are based on MCMC draws and $\hat{\theta}$ is the posterior mean for
the parameters. PPS(1) is subsampling proportional to $-\log(p(y_{i}|y))$,
PPS(2) is subsampling proportional to $-\log(p(y_{i}|\hat{\theta}))$, and SRS
is simple random sampling.

Table 2 shows empirical results on the effect of using subsampling
proportional to the predictive density compared to simple random sampling. The
results are much in line with what we would expect from the theory presented
in Section 2.3. We can see that the proposed method, sampling proportional to
$-\log(p(y_{i}|y))$ (PPS(1)) and sampling proportional to
$-\log(p(y_{i}|\hat{\theta}))$ (PPS(2)) outperforms simple random sampling
with orders of magnitude. Using just a sample size of $m=10$ observations
using our proposed method is much more precise than using $m=1000$
observations with simple random sampling, although we can see that the
estimate $\hat{\sigma}_{\text{loo}}$ is less reliable for such small sample
sizes. This results can be explained by the sampling probabilities used in the
subsampling procedure. Figure 1 show the distribution of sampling
probabilities where we can see that the probabilities are highly skewed,
indicating the reason for the inefficiency of the SRS compared to the proposed
approach. Table 2 also show that sampling with
$\tilde{\pi}\propto-\log(p(y_{i}|\hat{\theta}))$ does not cost us very much in
precision when estimating $\text{elpd}_{\text{loo}}$. In many situations with
larger data we would expect that using a point estimate, $\hat{\theta}$, of
the parameters in computing the likelihood values would be much faster than
computing $-\log(p(y_{i}|y))$.

Figure 1: Sampling probabilities ($\tilde{\pi}$) for the LR(cor) 100D data and $\pi_{i}=-\log p(y_{i}|y)$. The results are very similar for the other LR data and for $\pi_{i}=-\log p(y_{i}|\hat{\theta})$. $m$ | $n$ | PPS | SRS | $\sigma_{\text{loo}}$
---|---|---|---|---
10 | 100 | 3.6 | 16 | 8
| 1000 | 2.7 | 361 | 22
| 10000 | 3.8 | 9351 | 72
| 100000 | 17.5 | 19721 | 225
100 | 100 | 1.2 | 10 | 8
| 1000 | 1.3 | 84 | 22
| 10000 | 1.4 | 1144 | 72
| 100000 | 6.2 | 5894 | 225
Table 3: Standard errors, SE($\widehat{\text{elpd}}_{\text{loo}}$), for PPS
and SRS subsampling in relationship with $\sigma_{\text{loo}}$. The result are
based on MCMC draws and $\hat{\theta}$ is the posterior mean for the
parameters. PPS is subsampling proportional to $-\log(p(y_{i}|\hat{\theta}))$,
and SRS is simple random sampling.

Table 3 shows empirical results on the scaling characteristics of the proposed
method. We can see that for the PPS estimator, as the size of the data, $n$,
increases, the variance of the estimator is more or less constant. Using a SRS
sampling scheme, on the other hand, clearly show that to estimate the total
$\text{elpd}_{\text{loo}}$, we would need to increase the sample size, $m$, as
the number of observations, $n$, increases.

###  3.3 Hierarchical models for radon measurements

As an example of how the proposed method can be used, we exemplify with the
dataset of (Lin et al., 1999), used as an example of hierarchical modeling by
Gelman & Hill (2006).111We base our example and data on the Stan case study by
Chris Fonnesbeck at https://mc-stan.org/users/documentation/case-
studies/radon.html The data make up a total of 12 573 home radon measurements
in a total of 386 counties with a different number of observations per county.
This example is enlightening for a number of reasons. First, it is large
enough to actually take some computing time to analyze, but is still small
enough so we can use MCMC to compute the full data $\text{elpd}_{\text{loo}}$
as gold standard. The models used here are also both regular and singular,
showing the usability in a broader class of models. Finally, this example also
shows how we can mix different approximation techniques for different models
when doing model comparisons.

We compare seven different linear models of predicting the log radon levels in
individual houses based on floor measurements and county uranium levels. The
seven models are a pooled simple linear model (model 1), a non-pooled model
with one intercept estimated per county (model 2), a partially pooled model
with a hierarchical mean parameter per county (model 3), a variable intercept
model per county (model 4), a variable slope model per county (model 5), a
variable intercept and slope model (model 6), and finally a model with a
county level features and county level intercepts using the log uranium level
in the county. We use vague priors based on the Stan prior choice
recommendations222See https://github.com/stan-dev/stan/wiki/Prior-Choice-
Recommendations with $N(0,10)$ priors on regression coefficients and
intercepts and half-$N(0,1)$ for variance parameters. We ran all models using
Laplace, ADVI(FR), ADVI(MF) and MCMC. We ran the ADVI approximations for 100
000 iterations and for MCMC we use the Stan standard dynamic HMC algorithm
using 2 chains and 500 warmup iterations per chain. For all models we used
1000 posterior samples to compute the $\text{elpd}_{\text{loo}}$.

$M$ | Laplace | ADVI(FR) | ADVI(MF)
---|---|---|---
1 | 0.24 | 0.23 | 0.34
2 | 0.93 | 5.11 | 0.45
3 | 1.44 | 4.19 | 0.28
4 | 2.05 | 6.99 | 0.71
5 | - | 5.62 | 1.04
6 | - | 10.99 | 2.98
7 | 1.70 | 7.39 | 0.89
Table 4: Posterior $\hat{k}$ values for the different radon models and the
model approximations. Laplace was not possible for model 5 and 6.

Table 4 show the $\hat{k}$ values for different approximations for the
different models. We can see that for the simplest model (1) we get a good
posterior approximation with just Laplace approximation, but as the models
become more complex (and singular) we need better approximation techniques
such as ADVI. We can see that ADVI (mean-field) in general perform well and
can be used for inference in many models, even though more complex models
(such as model 4-7) is not approximated sufficiently well using the mean-field
approximation. ADVI full-rank, again, have problems due to the larger number
of parameters in the more complex models.

$M$ | Method | $\widehat{\text{elpd}}_{\text{loo}}$ | $SE$ | $\text{elpd}_{\text{loo}}$ | $\text{elpd}_{\text{kcv}}$
---|---|---|---|---|---
1 | Laplace | -18560 | 0.3 | -18560 | -18561
| ADVI(FR) | -18562 | 1.0 | -18559 | -18560
| ADVI(MF) | -18564 | 2.1 | -18559 | -18558
| MCMC | -18559 | 0.3 | -18560 | -18559
2 | Laplace | -17058 | 31.9 | -17049 | -17142
| ADVI(MF) | -17068 | 50.1 | -17059 | -17105
| MCMC | -17064 | 30.7 | -17067 | -17110
3 | Laplace | -17035 | 33.0 | -17017 | -17117
| ADVI(MF) | -17097 | 20.7 | -17090 | -17090
| MCMC | -17068 | 17.4 | -17085 | -17102
4 | Laplace | -17003 | 66.6 | -16866 | -17057
| ADVI(MF) | -16990 | 19.8 | -17013 | -17034
| MCMC | -17043 | 19.5 | -17022 | -17049
5 | ADVI(MF) | -18223 | 37.6 | -18225 | -18285
| MCMC | -18295 | 52.2 | -18253 | -18308
6 | ADVI(MF) | -16656 | 90.8 | -16603 | -16869
| MCMC | -16835 | 63.2 | -16798 | -16864
7 | Laplace | -17096 | 45.8 | -17063 | -17140
| ADVI(MF) | -16996 | 25.4 | -16957 | -17035
| MCMC | -17126 | 27.3 | -17136 | -17050
Table 5: The estimated $\widehat{\text{elpd}}_{\text{loo}}$ using a subsample
of size $m=500$ and its standard error (SE). The full
$\text{elpd}_{\text{loo}}$ based on all observations is also included as well
as $\text{elpd}_{\text{kcv}}$, an estimation of the elpd using 10-fold cross-
validation. The $\sigma_{\text{loo}}\approx 90$ for all approximations and
models. Less than 1% of the observation have problematic $\hat{k}$ using MCMC,
making it a good gold standard.

Based on these approximate posteriors we can analyze the
$\text{elpd}_{\text{loo}}$ for the different models. Table 5 shows the
$\text{elpd}_{\text{loo}}$ and an estimate,
$\widehat{\text{elpd}}_{\text{loo}}$, based on a subsample of size 500. As a
comparison we also compute $\text{elpd}_{\text{kcv}}$, computing an estimate
of elpd using a 10-fold cross-validation scheme (without bias correction). For
the simple baseline model 1, we can use Laplace approximation and a subsample
to estimate the $\text{elpd}_{\text{loo}}$ in roughly 2 seconds with a
sufficient precision for most purposes. Using MCMC and computing the full elpd
take roughly 35 seconds for this medium-sized dataset.

Table 5 also shows that ADVI (mean-field) work well both for regular and
singular models. Using ADVI(MF) for the singular models 3 and 4, where the
$\hat{k}$ values indicating a good posterior approximation, the approach works
really well. We can also see that the $\hat{k}$ diagnostic works well as an
indicator. The Laplace and ADVI(MF) approximations with high $\hat{k}$ values
can be quite off, see model 7 for an example.

The results of Table 5 also give us an idea of how the subsampling can be
used. By comparing the SE of our estimates with $\sigma_{\text{loo}}$, that is
roughly 90 for all models, we see how far a subsample with 500 observations
takes us. For most models, our SE is small enough to help us decide between
models, while for the more complex models. The precision needed depends on the
specific use case and if we need better precision we can simply add more
subsamples to get the precision needed.

If we study Table 5 we see that using ADVI(MF), Laplace and a subsample of
size 500 we can get quite far comparing these models. We could quickly rule
out model 1 and 5, but where we would need to use MCMC for model 5, due to the
high $\hat{k}$ for the ADVI approximations. Model 4 and 6 are the most
promising but we need to estimate the models using MCMC due to the high
$\hat{k}$ values for the ADVI approximations. Although, based on just the
subsample, we can see that model 6, the variable intercept and slope model,
seem to be the most promising model for this data. Comparing the fully
computed $\text{elpd}_{\text{loo}}$ for the different models we could compute
the difference in $\text{elpd}_{\text{loo}}$ between model 6 and 4 to 220 with
a standard error of 26, clearly indicating that model 6 is the one to prefer
in this situation. Using 10-fold cross-validation (see
$\text{elpd}_{\text{kcv}}$), we arrive at a similar result, but at the cost of
re-estimating the model 10 times.

##  4 Conclusions

In this work we solve the two major hurdles for using leave-one-out cross-
validation for large data, namely using posterior approximations to estimate
the $\text{elpd}_{\text{loo}}$ for individual observations and efficient
subsampling. We prove the consistency in $n$ and also show that for regular
models we have consistency also for common posterior approximations such as
Laplace and ADVI, even for mean-field ADVI in situations with correlated
posteriors, making the results promising for large-scale model evaluations.
Finally, our proposed method also comes with diagnostics to assess if the
quality of the subsampling and posterior approximations. We can use the
$\hat{k}$ diagnostic to asses the posterior approximations and
$v(\widehat{\text{elpd}}_{\text{loo}})$, the variance of the HH estimator, to
give us a good measure of the uncertainty due to subsampling.

## References

  * Azevedo-Filho & Shachter (1994) Azevedo-Filho, Adriano and Shachter, Ross D.  Laplace’s method approximations for probabilistic inference in belief networks with continuous variables.  In _Uncertainty Proceedings 1994_ , pp. 28–36. Elsevier, 1994.
  * Bengio & Grandvalet (2004) Bengio, Yoshua and Grandvalet, Yves.  No unbiased estimator of the variance of k-fold cross-validation.  _Journal of machine learning research_ , 5(Sep):1089–1105, 2004.
  * Bernardo (1979) Bernardo, José M.  Expected information as expected utility.  _the Annals of Statistics_ , pp. 686–690, 1979.
  * Bernardo & Smith (1994) Bernardo, José M and Smith, Adrian FM.  _Bayesian theory_.  IOP Publishing, 1994.
  * Blei et al. (2017) Blei, David M, Kucukelbir, Alp, and McAuliffe, Jon D.  Variational inference: A review for statisticians.  _Journal of the American Statistical Association_ , 112(518):859–877, 2017.
  * Box (1976) Box, George EP.  Science and statistics.  _Journal of the American Statistical Association_ , 71(356):791–799, 1976.
  * Carpenter et al. (2017) Carpenter, Bob, Gelman, Andrew, Hoffman, Matthew D, Lee, Daniel, Goodrich, Ben, Betancourt, Michael, Brubaker, Marcus, Guo, Jiqiang, Li, Peter, and Riddell, Allen.  Stan: A probabilistic programming language.  _Journal of statistical software_ , 76(1), 2017.
  * Cochran (1977) Cochran, William G.  _Sampling Techniques, 3rd Edition._ John Wiley, 1977.
  * Gelfand (1996) Gelfand, Alan E.  Model determination using sampling-based methods.  _Markov chain Monte Carlo in practice_ , pp. 145–161, 1996.
  * Gelman & Hill (2006) Gelman, Andrew and Hill, Jennifer.  _Data analysis using regression and multilevel/hierarchical models_.  Cambridge university press, 2006.
  * Geweke (1989) Geweke, John.  Bayesian inference in econometric models using Monte Carlo integration.  _Econometrica: Journal of the Econometric Society_ , pp. 1317–1339, 1989.
  * Hansen & Hurwitz (1943) Hansen, Morris H. and Hurwitz, William N.  On the theory of sampling from finite populations.  _The Annals of Mathematical Statistics_ , 14(4):333–362, 12 1943.
  * Jordan et al. (1999) Jordan, Michael I, Ghahramani, Zoubin, Jaakkola, Tommi S, and Saul, Lawrence K.  An introduction to variational methods for graphical models.  _Machine learning_ , 37(2):183–233, 1999.
  * Kucukelbir et al. (2017) Kucukelbir, Alp, Tran, Dustin, Ranganath, Rajesh, Gelman, Andrew, and Blei, David M.  Automatic differentiation variational inference.  _The Journal of Machine Learning Research_ , 18(1):430–474, 2017.
  * Lin et al. (1999) Lin, Chia-yu, Gelman, Andrew, Price, Phillip N, and Krantz, David H.  Analysis of local decisions using hierarchical modeling, applied to home radon measurement and remediation.  _Statistical Science_ , pp. 305–328, 1999.
  * Merkle et al. (2018) Merkle, EC, Furr, D, and Rabe-Hesketh, S.  Bayesian model assessment: Use of conditional vs marginal likelihoods.  _arXiv preprint arXiv:1802.04452_ , 2018.
  * Robert (1996) Robert, Christian P.  Intrinsic losses.  _Theory and decision_ , 40(2):191–214, 1996.
  * Stan Development Team (2018) Stan Development Team.  The Stan Core Library, 2018.  URL http://mc-stan.org/.  Version 2.18.0.
  * Vehtari & Ojanen (2012) Vehtari, Aki and Ojanen, Janne.  A survey of Bayesian predictive methods for model assessment, selection and comparison.  _Statistics Surveys_ , 6:142–228, 2012.
  * Vehtari et al. (2015) Vehtari, Aki, Gelman, Andrew, and Gabry, Jonah.  Pareto smoothed importance sampling.  _arXiv preprint arXiv:1507.02646_ , 2015.
  * Vehtari et al. (2017) Vehtari, Aki, Gelman, Andrew, and Gabry, Jonah.  Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC.  _Statistics and Computing_ , 27(5):1413–1432, 2017.
  * Vehtari et al. (2018) Vehtari, Aki, Gelman, Andrew, Gabry, Jonah, Yao, Yuling, Piironen, Juho, and Goodrich, Ben.  loo: Efficient leave-one-out cross-validation and WAIC for Bayesian models.  _R package version 2.0.0_ , 2018.
  * Walker (1977) Walker, Alastair J.  An efficient method for generating discrete random variables with general distributions.  _ACM Transactions on Mathematical Software (TOMS)_ , 3(3):253–256, 1977.
  * Walker (1969) Walker, Andrew M.  On the asymptotic behaviour of posterior distributions.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , pp. 80–88, 1969.
  * Wang & Blei (2018) Wang, Yixin and Blei, David M.  Frequentist consistency of variational Bayes.  _Journal of the American Statistical Association_ , (just-accepted):1–85, 2018.
  * Watanabe (2010) Watanabe, Sumio.  Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory.  _Journal of Machine Learning Research_ , 11(Dec):3571–3594, 2010.
  * Yao et al. (2018) Yao, Yuling, Vehtari, Aki, Simpson, Daniel, and Gelman, Andrew.  Yes, but did it work?: Evaluating variational inference.  In _Proceedings of the 35th International Conference on Machine Learning_ , pp. 5581–5590, 2018.

## Appendix

## Proof of Proposition 1

A generic Bayesian model is considered; a sample $(y_{1},y_{2},\ldots,y_{n})$,
$y_{i}\in\mathcal{Y}\subseteq{\mathbb{R}}$, is drawn from a true density
$p_{t}=p(\cdot|\theta_{0})$ for some true parameter $\theta_{0}$. The
parameter $\theta_{0}$ is assumed to be drawn from a prior $p(\theta)$ on the
parameter space $\Theta$, which we assume to be an open and bounded subset of
${\mathbb{R}}^{d}$.

A number of conditions are used. They are as follows.

  * (i)

the likelihood $p(y|\theta)$ satisfies that there is a function
$C:\mathcal{Y}\rightarrow{\mathbb{R}}_{+}$, such that ${\mathbb{E}}_{y\sim
p_{t}}[C(y)^{2}]<\infty$ and such that for all $\theta_{1}$ and $\theta_{2}$,
$|p(y|\theta_{1})-p(y|\theta_{2})|\leq
C(y)p(y|\theta_{2})\|\theta_{1}-\theta_{2}\|$.

  * (ii)

$p(y|\theta)>0$ for all $(y,\theta)\in\mathcal{Y}\times\Theta$,

  * (iii)

There is a constant $M<\infty$ such that $p(y|\theta)<M$ for all $(y,\theta)$,

  * (iv)

all assumptions needed in the Bernstein-von Mises (BvM) Theorem (Walker,
1969),

  * (v)

for all $\theta$, $\int_{\mathcal{Y}}(-\log p(y|\theta))p(y|\theta)dy<\infty$.

Remarks.

  * •

There are alternatives or relaxations to (i) that also work. One is to assume
that there is an $\alpha>0$ and $C$ with ${\mathbb{E}}_{y}[C(y)^{2}]<\infty$
such that $|p(y|\theta_{1})-p(y|\theta_{2})|\leq
C(y)p(y|\theta_{2})\|\theta_{1}-\theta_{2}\|^{\alpha}$. There are many
examples when (i) holds, e.g. when $y$ is normal, Laplace distributed or
Cauchy distributed with $\theta$ as a one-dimensional location parameter.

  * •

The assumption that $\Theta$ is bounded will be used solely to draw the
conclusion that ${\mathbb{E}}_{y,\theta}\|\theta-\theta_{0}\|\rightarrow 0$ as
$n\rightarrow\infty$, where $y$ is the sample and $\theta$ is either
distributed according to the true posterior (which is consistent by BvM) or
according to a consistent approximate posterior. The conclusion is valid by
the definition of consistency and the fact that the boundedness of $\Theta$
makes $\|\theta-\theta_{0}\|$ a bounded function of $\theta$. If it can be
shown by other means for special cases that
${\mathbb{E}}_{y,\theta}\|\theta-\theta_{0}\|\rightarrow 0$ despite $\Theta$
being unbounded, then our results also hold.

  * •

We can (and will) without loss of generality assume that $M=1/2$ is sufficient
in (iii), for if not then simply transform data and consider $z_{i}=2My_{i}$
instead of $y_{i}$.

The main quantity of interest is the mean expected log pointwise predictive
density, which we want to use for model evaluation and comparison.

######  Definition 1 ($\overline{\text{elpd}}$).

The _mean expected log pointwise predictive density_ for a model $p$ is
defined as

$\overline{\text{elpd}}=\int p_{t}(x)\log p(x)\,dx$

where $p_{t}(x)=p(x|\theta_{0})$ is the _true_ density at a new unseen
observation $x$ and $\log p(x)$ is the log predictive density for observation
$x$.

We estimate $\overline{\text{elpd}}$ using leave-one-out cross-validation
(loo).

######  Definition 2 (Leave-one-out cross-validation).

The loo estimator $\overline{\text{elpd}}_{loo}$ is given by

$\overline{\text{elpd}}_{loo}=\frac{1}{n}\sum^{n}_{i=1}\log p(y_{i}|y_{-i}),$
(7)

where $p(y_{i}|y_{-i})=\int p(y_{i}|\theta)p(\theta|y_{-i})d\theta$.

To estimate $\overline{\text{elpd}}_{loo}$ in turn, we use importance sampling
and the Hansen-Hurwitz estimator. Definitions follow.

######  Definition 3.

The Hansen-Hurwitz estimator is given by

$\widehat{\overline{\text{elpd}}}_{\text{loo}}(m,q)=\frac{1}{m}\frac{1}{n}\sum_{j=1}^{m}\frac{1}{\tilde{\pi}_{j}}\log\hat{p}(y_{j}|y_{-j})$

where $\tilde{\pi}_{i}$ is the probability of subsampling observation $i$,
$\log\hat{p}(y_{i}|y_{-i})$ is the (self-normalized) importance sampling
estimate of $\log p(y_{i}|y_{-i})$ defined as

$\log\hat{p}(y_{i}|y_{-i})=\log\left(\frac{\frac{1}{S}\sum_{s=1}^{S}p(y_{i}|\theta_{s})r(\theta_{s})}{\frac{1}{S}\sum_{s=1}^{S}r(\theta_{s})}\right)\,,$

where

$\displaystyle r(\theta_{s})=$
$\displaystyle\frac{p(\theta_{s}|y_{-i})}{p(\theta_{s}|y)}\frac{p(\theta_{s}|y)}{q(\theta_{s}|y)}$
$\displaystyle\propto$
$\displaystyle\frac{1}{p(y_{i}|\theta_{s})}\frac{p(\theta_{s}|y)}{q(\theta_{s}|y)}$

and where $q(\theta|y)$ is an approximation of the posterior distribution,
$\theta_{s}$ is a sample from the approximate posterior distribution
$q(\theta|y)$ and S is the total posterior sample size.

######  Proposition 2.

Let the subsampling size $m$ and the number of posterior draws $S$ be fixed at
arbitrary integer numbers, let the sample size $n$ grow, assume that (i)-(vi)
hold and let $q=q_{n}(\cdot|y)$ be any consistent approximate posterior. Write
$\hat{\theta}_{q}=\arg\max\\{q(\theta):\theta\in\Theta\\}$ and assume further
that $\hat{\theta}_{q}$ is a consistent estimator of $\theta_{0}$. Then

$|\widehat{\overline{\text{elpd}}}_{\text{loo}}(m,q)-\overline{\text{elpd}}_{loo}|\rightarrow
0$

in probability as $n\rightarrow\infty$ for any of the following choices of
$\pi_{i}$, $i=1,\ldots,n$.

  * (a)

$\pi_{i}=-\log p(y_{i}|y)$,

  * (b)

$\pi_{i}=-{\mathbb{E}}_{y}[\log p(y_{i}|y)]$,

  * (c)

$\pi_{i}=-{\mathbb{E}}_{\theta\sim q}[\log p(y_{i}|\theta)]$,

  * (d)

$\pi_{i}=-\log p(y_{i}|{\mathbb{E}}_{\theta\sim q}[\theta])$,

  * (e)

$\pi_{i}=-\log p(y_{i}|\hat{\theta}_{q})$.

Remark. By the variational BvM Theorems of Wang and Blei, (Wang & Blei, 2018),
$q$ can be taken to be either $q_{Lap}$, $q_{MF}$ or $q_{FR}$, i.e. the
approximate posteriors of the Laplace, mean-field or full-rank variational
families respectively in Proposition 2, provided that one adopts the mild
conditions in their paper.

The proof of Proposition 2 will be focused on proving (a) and then (b)-(e)
will follow easily. We begin with the following key lemma.

######  Lemma 3.

With all quantities as defined above,

${\mathbb{E}}_{y\sim p_{t}}|\pi_{i}-\log p(y_{i}|\theta_{0})|\rightarrow 0,$
(8)

with any of the definitions (a)-(e) of $\pi_{i}$ of Proposition 2.
Furthermore,

${\mathbb{E}}_{y\sim p_{t}}|\log p(y_{i}|y_{-i})-\log
p(y_{i}|\theta_{0})|\rightarrow 0,$ (9)

and

${\mathbb{E}}_{y\sim p_{t}}|\log\hat{p}(y_{i}|y)-\log
p(y_{i}|\theta_{0})|\rightarrow 0.$ (10)

as $n\rightarrow\infty$.

###### Proof.

To avoid burdening the notation unnecessarily, we write throughout the proof
${\mathbb{E}}_{y}$ for ${\mathbb{E}}_{y\sim p_{t}}$. For now, we also write
${\mathbb{E}}_{\theta}$ as shorthand for ${\mathbb{E}}_{\theta\sim
p(\cdot|y_{-i})}$. Recall that $x_{+}=\max(x,0)=ReLU(x)$.

Hence

$\displaystyle{\mathbb{E}}_{y}\left[\left(\log\frac{p(y_{i}|y_{-i})}{p(y_{i}|\theta_{0})}\right)_{\operatorname*{\scalerel*{+}{\sum}}}\right]$
$\displaystyle={\mathbb{E}}_{y}\left[\left(\log\frac{{\mathbb{E}}_{\theta}[p(y_{i}|\theta)]}{p(y_{i}|\theta_{0})}\right)_{\operatorname*{\scalerel*{+}{\sum}}}\right]$
$\displaystyle\leq{\mathbb{E}}_{y}\left[\log\left(1+\frac{{\mathbb{E}}_{\theta}\left[C(y_{i})p(y_{i}|\theta_{0})\|\theta-\theta_{0}\|\right]}{p(y_{i}|\theta_{0})}\right)\right]$
$\displaystyle\leq{\mathbb{E}}_{y,\theta}[C(y_{i})\|\theta-\theta_{0}\|]$
$\displaystyle\leq\left({\mathbb{E}}_{y_{i}}[C(y_{i})^{2}]{\mathbb{E}}_{y,\theta}\left[\|\theta-\theta_{0}\|^{2}\right]\right)^{1/2}$
$\displaystyle\rightarrow 0\mbox{ as }n\rightarrow\infty.$

Here the first inequality follows from condition (i) and the second inequality
from the fact that $\log(1+x)<x$ for $x\geq 0$. The third inequality is
Schwarz inequality. The limit conclusion follows from the consistency of the
posterior $p(\cdot|y_{-i})$ and the definition of weak convergence, since
$\|\theta-\theta_{0}\|^{2}$ is a continuous bounded function of $\theta$
(recall that $\Theta$ is bounded) and that the first factor is finite by
condition (i).

For the reverse inequality,

$\displaystyle{\mathbb{E}}_{y}\left[\left(\log\frac{p(y_{i}|\theta_{0})}{p(y_{i}|y_{-i})}\right)_{\operatorname*{\scalerel*{+}{\sum}}}\right]$
$\displaystyle={\mathbb{E}}_{y}\left[\left(\log{\mathbb{E}}_{\theta}\left[\frac{p(y_{i}|\theta_{0})]}{p(y_{i}|\theta)}\right]\right)_{\operatorname*{\scalerel*{+}{\sum}}}\right]$
$\displaystyle\leq{\mathbb{E}}_{y}\left[\log\left(1+{\mathbb{E}}_{\theta}\left[\frac{C(y_{i})p(y_{i}|\theta)\|\theta-\theta_{0}\|}{p(y_{i}|\theta)}\right]\right)\right]$
$\displaystyle\leq\left({\mathbb{E}}_{y_{i}}[C(y_{i})^{2}]{\mathbb{E}}_{y,\theta}\left[\|\theta-\theta_{0}\|^{2}\right]\right)^{1/2}$
$\displaystyle\rightarrow 0\mbox{ as }n\rightarrow\infty.$

This proves (9) and an identical argument proves (8) for $\pi_{i}=p(y_{i}|y)$.

For $\pi_{i}=-{\mathbb{E}}_{y}[\log p(y_{i}|y)]$, note first that

$\displaystyle{\mathbb{E}}_{y}\left|{\mathbb{E}}_{y}[\log
p(y_{i}|y)]-{\mathbb{E}}_{y}[\log p(y_{i}|y_{-i})]\right|$
$\displaystyle=\left|{\mathbb{E}}_{y}[\log p(y_{i}|y)-\log
p(y_{i}|y_{-i})]\right|$ $\displaystyle\leq{\mathbb{E}}_{y}\left|\log
p(y_{i}|y)-\log p(y_{i}|y_{-i})]\right|$

which goes to $0$ by (9) and (a). Hence we can replace
$\pi_{i}=-{\mathbb{E}}[\log p(y_{i}|y)]$ with $\pi_{i}=-{\mathbb{E}}[\log
p(y_{i}|y_{-i})]$ when proving (b). To that end, observe that

$\displaystyle\left({\mathbb{E}}_{y}[\log p(y_{i}|y_{-i})]-\log
p(y_{i}|\theta_{0})\right)_{+}$
$\displaystyle=\left({\mathbb{E}}_{y_{i}}\left[{\mathbb{E}}_{y_{-i}}\left[\log\frac{p(y_{i}|y_{-i})}{p(y_{i}|\theta_{0})}\right]\right]\right)_{\operatorname*{\scalerel*{+}{\sum}}}$
$\displaystyle\leq{\mathbb{E}}_{y}\left[\left(\log\frac{p(y_{i}|y_{-i})}{p(y_{i}|\theta_{0})}\right)_{\operatorname*{\scalerel*{+}{\sum}}}\right].$

where the inequality is Jensen’s inequality used twice on the convex function
$x\rightarrow x_{+}$. Now everything is identical to the proof of (9) and the
reverse inequality is analogous.

The other choices of $\pi_{i}$ follow along very similar lines. For
$\pi_{i}=-\log p(y_{i}|\hat{\theta}_{q})$, we have on mimicking the above that

$\displaystyle{\mathbb{E}}_{y}\left[\left(\log\frac{p(y_{i}|\hat{\theta}_{q})}{p(y_{i}|\theta_{0})}\right)_{\operatorname*{\scalerel*{+}{\sum}}}\right]$
$\displaystyle\leq\left({\mathbb{E}}_{y_{i}}[C(y_{i})^{2}]{\mathbb{E}}_{y}\left[\|\hat{\theta}_{q}-\theta_{0}\|^{2}\right]\right)^{1/2}$

and ${\mathbb{E}}_{y}[\|\hat{\theta}_{q}-\theta_{0}\|^{2}]\rightarrow 0$ as
$n\rightarrow\infty$ by the assumed consistency of $\hat{\theta}_{q}$. The
reverse inequality is analogous and (8) for
$\pi_{i}=p(y_{i}|\hat{\theta}_{q})$ is established.

For the case $\pi_{i}=-\log p(y_{i}|{\mathbb{E}}_{\theta\sim q}\theta)$, the
analogous analysis gives

$\displaystyle{\mathbb{E}}_{y}\left[\left(\log\frac{p(y_{i}|{\mathbb{E}}_{\theta\sim
q}\theta)}{p(y_{i}|\theta_{0})}\right)_{\operatorname*{\scalerel*{+}{\sum}}}\right]$
$\displaystyle\leq{\mathbb{E}}_{y_{i}}[C(y_{i})^{2}]{\mathbb{E}}_{y}[\|{\mathbb{E}}_{\theta\sim
q}\theta-\theta_{0}\|^{2}].$

Since $x\rightarrow\|x-\theta_{0}\|^{2}$ is convex, the second factor on the
right hand side is bounded by ${\mathbb{E}}_{y,\theta\sim
q}[\|\theta-\theta_{0}\|^{2}]$ which goes to 0 by the consistency of $q$ and
the boundedness of $\Theta$. The reverse inequality is again analogous.

Finally for $\pi_{i}=-{\mathbb{E}}_{\theta\sim q}[\log p(y_{i}|\theta)]$,

$\displaystyle{\mathbb{E}}_{y}\left[\left({\mathbb{E}}_{\theta\sim q}[\log
p(y_{i}|\theta)]-\log
p(y_{i}|\theta_{0})\right)_{\operatorname*{\scalerel*{+}{\sum}}}\right]$
$\displaystyle={\mathbb{E}}_{y}\left[\left({\mathbb{E}}_{\theta\sim
q}\left[\log\frac{p(y_{i}|\theta)}{p(y_{i}|\theta_{0})}\right]\right)_{\operatorname*{\scalerel*{+}{\sum}}}\right]$
$\displaystyle\leq{\mathbb{E}}_{y,\theta\sim
q}\left[\left(\log\frac{p(y_{i}|\theta)}{p(y_{i}|\theta_{0})}\right)_{\operatorname*{\scalerel*{+}{\sum}}}\right]$
$\displaystyle\leq\left({\mathbb{E}}_{y_{i}}[C(y_{i})^{2}]{\mathbb{E}}_{y,\theta\sim
q}[\|\theta-\theta_{0}\|^{2}]\right)^{1/2}\rightarrow 0$

as $n\rightarrow\infty$ by the consistency of $q$. Here the first inequality
is Jensen’s inequality applied to $x\rightarrow x_{+}$ and the second
inequality follows along the same lines as before.

For (10), write
$r^{\prime}(\theta_{s})=r(\theta_{s})/\sum_{j=1}^{S}r(\theta_{j})$ for the
random weights given to the individual $\theta_{s}$:s in the expression for
$\hat{p}(y_{i}|y_{-i})$. Then we have, with
$\theta=(\theta_{1},\ldots,\theta_{S})$ chosen according to $q$,

$\displaystyle{\mathbb{E}}_{y}\left[\left(\log\frac{\hat{p}(y_{i}|y_{-i})}{p(y_{i}|\theta_{0})}\right)_{\operatorname*{\scalerel*{+}{\sum}}}\right]$
$\displaystyle={\mathbb{E}}_{y,\theta}\left[\left(\log\frac{\sum_{s=1}^{S}r^{\prime}(\theta_{s})p(y_{i}|\theta_{s})}{p(y_{i}|\theta_{0})}\right)_{\operatorname*{\scalerel*{+}{\sum}}}\right]$
$\displaystyle\leq{\mathbb{E}}_{y,\theta}\left[\log\left(1+\frac{\sum_{s=1}^{S}r^{\prime}(\theta_{s})|p(y_{i}|\theta_{s})-p(y_{i}|\theta_{0})|}{p(y_{i}|\theta_{0})}\right)\right]$
$\displaystyle\leq{\mathbb{E}}_{y,\theta}\left[\log\left(1+C(y_{i})\sum_{s=1}^{S}r^{\prime}(\theta_{s})\|\theta_{s}-\theta_{0}\|\right)\right]$
$\displaystyle\leq{\mathbb{E}}_{y,\theta}\left[\log\left(1+C(y_{i})\sum_{s=1}^{S}\|\theta_{s}-\theta_{0}\|\right)\right]$
$\displaystyle\leq{\mathbb{E}}_{y,\theta}\left[C(y_{i})\sum_{s=1}^{S}\|\theta_{s}-\theta_{0}\|\right]$
$\displaystyle\leq\left({\mathbb{E}}_{y_{i}}[C(y_{i})^{2}]{\mathbb{E}}_{y,\theta}\left[\left(\sum_{s=1}^{S}\|\theta_{s}-\theta_{0}\|\right)^{2}\right]\right)^{1/2},$

where the second inequality is condition (i) and the limit conclusion follows
from the consistency of $q$. For the reverse inequality to go through
analogously, observe that

$\displaystyle\frac{\left|p(y_{i}|\theta_{0})-\sum_{s}r^{\prime}(\theta_{s})p(y_{i}|\theta_{s})\right|}{\sum_{s}r^{\prime}(\theta_{s})p(y_{i}|\theta_{s})}$
$\displaystyle\leq\frac{\sum_{s}r^{\prime}(\theta_{s})|p(y_{i}|\theta_{s})-p(y_{i}|\theta_{0})|}{\sum_{s}r^{\prime}(\theta_{s})p(y_{i}|\theta_{s})}$
$\displaystyle\leq\frac{\sum_{s}r^{\prime}(\theta_{s})p(y_{i}|\theta_{s})\|\theta_{s}-\theta_{0}\|}{\sum_{s}r^{\prime}(\theta_{s})p(y_{i}|\theta_{s})}$
$\displaystyle\leq\max_{s}\|\theta_{s}-\theta_{0}\|$
$\displaystyle\leq\sum_{s}\|\theta_{s}-\theta_{0}\|.$

Equipped with this observation, mimic the above.

∎

For convenience we will write
$\hat{e}:=\hat{e}_{m,q}=\widehat{\overline{\text{elpd}}}_{loo}$, which for our
purposes is more usefully expressed as

$\hat{e}=\frac{1}{n}\frac{1}{m}\sum_{i=1}^{n}\sum_{j=1}^{m}I_{ij}\frac{1}{\bar{\pi}_{i}}\log\hat{p}(y_{i}|y_{-i}),$

where $I_{ij}$ is the indicator that sample point $y_{i}$ is chosen in draw
$j$ for the subsample used in $\hat{e}$. Write also

$e=\frac{1}{n}\frac{1}{m}\sum_{i=1}^{n}\sum_{j=1}^{m}I_{ij}\frac{1}{\bar{\pi}_{i}}\log
p(y_{i}|y_{-i}).$

In other words, $e$ is the HH estimator with $\hat{p}$ replaced with $p$.

######  Lemma 4.

With the notation as just defined and $\pi_{i}=-\log p(y_{i}|y)$,

${\mathbb{E}}|\hat{e}-e|\rightarrow 0$

as $n\rightarrow\infty$.

###### Proof.

We have, with expectations with respect to all sources of randomness involved
in $\hat{e}$ and $e$

$\displaystyle{\mathbb{E}}|\hat{e}-e|$ $\displaystyle\leq\,$
$\displaystyle\frac{1}{m}\frac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{m}{\mathbb{E}}\left[{\mathbb{E}}\left[I_{ij}\frac{1}{\bar{\pi_{i}}}|\log\hat{p}(y_{i}|y_{-i})-\log
p(y_{i}|y_{-i})|\Big{|}y\right]\right]$ $\displaystyle=\,$
$\displaystyle{\mathbb{E}}\left[\frac{1}{n}\frac{1}{m}\sum_{i=1}^{n}\sum_{j=1}^{m}|\log\hat{p}(y_{i}|y_{-i})-\log
p(y_{i}|y_{-i})|\right]$ $\displaystyle=\,$
$\displaystyle{\mathbb{E}}|\log\hat{p}(y_{i}|y_{-i})-\log p(y_{i}|y_{-i})|.$

The result now follows from (9), (10) and the triangle inequality. ∎

Proof of Proposition 2. As stated before, we start with a focus on (a), which
means that for now we have $\pi_{i}=-\log p(y_{i}|y)$ By Lemma 4, it suffices
to prove that $|e-\overline{\text{elpd}}_{loo}|\rightarrow 0$ in probability
with $\pi_{i}$ chosen according to any of (a)-(e). The variance of a HH
estimator is well known and some easy manipulation then tells us that the
conditional variance of $e$ given $y$ is given by

$V(e)={\mathbb{V}}{\rm
ar}(e|y)=\frac{1}{n^{2}}\frac{1}{m}(S_{\pi}S_{2}-S_{p}^{2}),$

where $S_{p}=\sum_{i=1}^{n}p_{i}$, $S_{\pi}=\sum_{i=1}^{n}\pi_{i}$ and
$S_{2}=\sum_{i=1}^{n}(p_{i}^{2}/\pi_{i})$. We claim that for any $\delta>0$,
for $n$ sufficiently large, ${\mathbb{P}}_{y}(V(e)<\delta)>1-\delta$. To this
end, observe first that

$\displaystyle{\mathbb{E}}_{y}[-\log p(y_{i}|y)]$ $\displaystyle\leq\,$
$\displaystyle{\mathbb{E}}_{y}[-\log
p(y_{i}|\theta_{0})]+{\mathbb{E}}_{y}|\log p(y_{i}|y)-\log
p(y_{i}|\theta_{0})|$ $\displaystyle\leq\,$
$\displaystyle{\mathbb{E}}_{y}[-\log p(y_{i}|\theta_{0})]+\delta<\infty$

for sufficiently large $n$, since the first term is finite by condition (v).
Let $A=A_{n}={\mathbb{E}}_{y}[-\log p(y_{i}|y)]$.

Now,

${\mathbb{E}}_{y}\left[\frac{1}{n}|S_{p}-S_{\pi}|\right]={\mathbb{E}}_{y}\left[\frac{1}{n}\left|\sum_{i=1}^{n}\pi_{i}-\sum_{i=1}^{n}p_{i}\right|\right]\rightarrow
0$

as $n\rightarrow\infty$ by (8) and (9). Hence for arbitrary $\alpha>0$,
${\mathbb{P}}_{y}(|S_{p}-S_{\pi}|<\alpha^{2}n)>1-\alpha$ for $n$ large enough.
Also

$\frac{p_{i}^{2}}{\pi_{i}}\leq\frac{(\pi_{i}+|p_{i}-\pi_{i}|)^{2}}{\pi_{i}}<\pi_{i}+4|\pi_{i}-p_{i}|$

(the last inequality using condition (iii): $\pi_{i}\geq-\log(1/2)>1/2$), so
$n^{-1}{\mathbb{E}}_{y}|S_{\pi}-S_{2}|\rightarrow 0$ and so
${\mathbb{P}}_{y}(|S_{p}-S_{2}|<\alpha^{2}n)>1-\alpha$ for sufficiently large
$n$. Hence with probability exceeding $1-2\alpha$, $y$ will be such that for
sufficiently large $n$,

$\displaystyle V(e)$
$\displaystyle\leq\frac{1}{n^{2}}\frac{1}{m}\left((S_{p}+\alpha^{2}n)^{2}-S_{p}^{2}\right)$
$\displaystyle=\frac{1}{n^{2}}\frac{1}{m}(2\alpha^{2}nS_{p}+\alpha^{4}n^{2}).$

We had ${\mathbb{E}}_{y}[S_{p}]=An$ and Markov’s inequality thus entails that
${\mathbb{P}}_{y}(S_{p}<An/\alpha)>1-\alpha$. Adding this piece of information
to the above, we get that with probability larger than $1-3\alpha$, $y$ will
for sufficiently large $n$ be such that

$V(e)\leq(2\alpha+\alpha^{4})n^{2}<3\alpha.$

For such $y$, Chebyshev’s inequality gives

${\mathbb{P}}(|e-{\mathbb{E}}[e|y]|>\alpha^{1/2}|y)<3\alpha^{1/2}.$

The HH estimator is unbiased, so
${\mathbb{E}}[e|y]=\overline{\text{elpd}}_{loo}$. We get for arbitrary
$\epsilon>0$ on taking $\alpha$ sufficiently small and $n$ correspondingly
large, taking all randomness into account

${\mathbb{P}}(|e-\overline{\text{elpd}}_{loo}|>\epsilon)<1-\epsilon$

which entails that $|e-\overline{\text{elpd}}_{loo}|\rightarrow 0$ in
probability. As observed above, this proves (a).

For the remaining parts, write $e_{p}$ when taking $\pi_{i}$ in $e$ according
to statement (p) in the proposition. By (8),
${\mathbb{E}}|e_{p}-e_{a}|\rightarrow 0$ for $p=b,c,d,e$ and we are done.

$\Box$

## Unbiasness of using the Hansen-Hurwitz estimator

###  4.1 On the Hansen-Hurwitz estimator

Let $\mathcal{Y}=\left\\{y_{1},y_{2},\ldots,y_{N}\right\\}$ be a set of non-
negative observations, $y_{i}>0$ and let
$\pi=\left\\{\pi_{1},\pi_{2},\ldots,\pi_{N}\right\\}$ be a probability vector
s.t. $\sum\pi_{j}=1$. Furthermore, let $a_{k}\in\left\\{1,2,\ldots,N\right\\}$
be i.i.d. samples from a multinomial distribution with probabilities $\pi$,
i.e. $a_{k}\stackrel{{\scriptstyle
iid}}{{\sim}}\text{Multinomial}\left(\pi\right)$.

We want to estimate the total

$\displaystyle\tau=\sum_{n=1}^{N}y_{i}$ (11)

using the Hansen-Hurwitz estimator given by

$\displaystyle\hat{\tau}=\frac{1}{M}\sum_{m=1}^{M}\frac{x_{m}}{p_{m}},$ (12)

where $x_{m}\equiv y_{a_{m}}$, $p_{m}\equiv\pi_{a_{m}}$, and
$a_{m}\sim\text{Multinomial}\left(\pi\right)$.

We can decompose $x_{m}$ and $p_{m}$ as follows

$\displaystyle x_{m}\equiv y_{a_{m}}$
$\displaystyle=\sum_{j=1}^{N}\mathbb{I}\left[a_{m}=j\right]y_{j}$ (13)
$\displaystyle p_{m}\equiv p_{a_{m}}$
$\displaystyle=\sum_{j=1}^{N}\mathbb{I}\left[a_{m}=j\right]\pi_{j}$ (14)

###  4.2 The Hansen-Hurwitz estimator is unbiased

First, we will show that the HH estimator, $\hat{\tau}$, is unbiased. We have,

$\displaystyle\mathbb{E}\left[\hat{\tau}\right]$
$\displaystyle=\mathbb{E}\left[\frac{1}{M}\sum_{m=1}^{M}\frac{x_{m}}{p_{m}}\right]=\frac{1}{M}\sum_{m=1}^{M}\mathbb{E}\left[\frac{x_{m}}{p_{m}}\right]$
(15)

Using the definitions in eq. (13) and (14) yields

$\displaystyle\mathbb{E}\left[\hat{\tau}\right]$
$\displaystyle=\frac{1}{M}\sum_{m=1}^{M}\mathbb{E}\left[\frac{\sum_{j=1}^{N}\mathbb{I}\left[a_{m}=j\right]y_{j}}{\sum_{j=1}^{N}\mathbb{I}\left[a_{m}=j\right]\pi_{j}}\right]$
$\displaystyle=\frac{1}{M}\sum_{m=1}^{M}\mathbb{E}\left[\sum_{j=1}^{N}\frac{y_{j}}{\pi_{j}}\mathbb{I}\left[a_{m}=j\right]\right]$
$\displaystyle=\frac{1}{M}\sum_{m=1}^{M}\sum_{j=1}^{N}\frac{y_{j}}{\pi_{j}}\mathbb{E}\left[\mathbb{I}\left[a_{m}=j\right]\right]$
$\displaystyle=\frac{1}{M}\sum_{m=1}^{M}\sum_{j=1}^{N}\frac{y_{j}}{\pi_{j}}\pi_{j}$
(16)

since
$\pi_{j}=\mathbb{P}\left[a_{m}=j\right]=\mathbb{E}\left[\mathbb{I}\left[a_{m}=j\right]\right]$.

Now it follows that

$\displaystyle\mathbb{E}\left[\hat{\tau}\right]$
$\displaystyle=\frac{1}{M}\sum_{m=1}^{M}\sum_{j=1}^{N}y_{j}=\sum_{j=1}^{N}y_{j}=\tau.$
(17)

###  4.3 An unbiased estimator of $\sigma^{2}_{\text{loo}}$

We also want to estimate the variance of the population $\mathcal{Y}$, i.e.

$\displaystyle\sigma_{y}^{2}=\frac{1}{N}\sum_{n=1}^{N}\left(y_{n}-\bar{y}\right)^{2},$
(18)

where $\bar{y}=\frac{1}{N}\sum y_{n}$.

First, we decompose the above as follows

$\displaystyle\sigma_{y}^{2}=\frac{1}{N}\sum_{n=1}^{N}y^{2}_{n}-\bar{y}^{2}.$
(19)

We will consider estimators for the two terms,
$\frac{1}{N}\sum_{n=1}^{N}y^{2}_{n}$ (1) and $\bar{y}^{2}$ (2), separately.
First, we will show that the following is an unbiased estimate of the first
term,

$\displaystyle T_{1}=\frac{1}{NM}\sum_{m=1}^{M}\frac{x_{m}^{2}}{p_{m}}.$ (20)

We have

$\displaystyle\mathbb{E}\left[T_{1}\right]=\mathbb{E}\left[\frac{1}{NM}\sum_{m=1}^{M}\frac{x_{m}^{2}}{p_{m}}\right]=\frac{1}{NM}\sum_{m=1}^{M}\mathbb{E}\left[\frac{x_{m}^{2}}{p_{m}}\right]$
(21)

Again, we use the representations in eq. (13) and (14) to get

$\displaystyle\mathbb{E}\left[\frac{1}{NM}\sum_{m=1}^{M}\frac{x_{m}^{2}}{p_{m}}\right]$
$\displaystyle=\frac{1}{NM}\sum_{m=1}^{M}\mathbb{E}\left[\frac{\sum_{j=1}^{N}\mathbb{I}\left[a_{m}=j\right]y_{j}^{2}}{\sum_{j=1}^{N}\mathbb{I}\left[a_{m}=j\right]\pi_{j}}\right]$
$\displaystyle=\frac{1}{NM}\sum_{m=1}^{M}\mathbb{E}\left[\sum_{j=1}^{N}\mathbb{I}\left[a_{m}=j\right]\frac{y_{j}^{2}}{\pi_{j}}\right]$
$\displaystyle=\frac{1}{NM}\sum_{m=1}^{M}\sum_{j=1}^{N}\frac{y_{j}^{2}}{\pi_{j}}\mathbb{E}\left[\mathbb{I}\left[a_{m}=j\right]\right]$
$\displaystyle=\frac{1}{NM}\sum_{m=1}^{M}\sum_{j=1}^{N}\frac{y_{j}^{2}}{\pi_{j}}\pi_{j}$
$\displaystyle=\frac{1}{N}\sum_{j=1}^{N}y_{j}^{2}.$ (22)

This completes the proof of for the first term.

For the second term, we use the estimator $T_{2}$ given by

$\displaystyle T_{2}=$
$\displaystyle\frac{1}{M(M-1)}\sum_{m=1}^{M}\left[\frac{x_{m}}{Np_{m}}-\frac{1}{N}\sum_{k=1}^{M}\frac{x_{k}}{Mp_{k}}\right]^{2}$
$\displaystyle-\left[\frac{1}{N}\sum_{k=1}^{M}\frac{x_{k}}{Mp_{k}}\right]^{2}.$
(23)

We have

$\displaystyle\frac{1}{M(M-1)}\sum_{m=1}^{M}\left[\frac{x_{m}}{Np_{m}}-\sum_{k=1}^{M}\frac{x_{k}}{NMp_{k}}\right]^{2}-\left[\sum_{k=1}^{M}\frac{x_{k}}{NMp_{k}}\right]^{2}$
$\displaystyle=\frac{1}{N^{2}M(M-1)}\sum_{m=1}^{M}\frac{x^{2}_{m}}{p^{2}_{m}}-\frac{1}{N^{2}M(M-1)}\left[\sum_{k=1}^{M}\frac{x_{k}}{p_{k}}\right]^{2}$
(24)

We consider now the expectation of the first term in the equation above

$\displaystyle\mathbb{E}\left[\sum_{m=1}^{M}\frac{x^{2}_{m}}{p^{2}_{m}}\right]$
$\displaystyle=\sum_{m=1}^{M}\mathbb{E}\left[\frac{x^{2}_{m}}{p^{2}_{m}}\right]$
$\displaystyle=\sum_{m=1}^{M}\mathbb{E}\left[\frac{\sum_{j=1}^{N}\mathbb{I}\left[a_{m}=j\right]y_{j}^{2}}{\sum_{j=1}^{N}\mathbb{I}\left[a_{m}=j\right]\pi_{j}^{2}}\right]$
$\displaystyle=\sum_{m=1}^{M}\mathbb{E}\left[\sum_{j=1}^{N}\mathbb{I}\left[a_{m}=j\right]\frac{y_{j}^{2}}{\pi_{j}^{2}}\right]$
$\displaystyle=\sum_{m=1}^{M}\sum_{j=1}^{N}\mathbb{E}\left[\mathbb{I}\left[a_{m}=j\right]\right]\frac{y_{j}^{2}}{\pi_{j}^{2}}$
$\displaystyle=M\sum_{j=1}^{N}\frac{y_{j}^{2}}{\pi_{j}}$ (25)

and the second term

$\displaystyle\mathbb{E}\left[\left[\sum_{k=1}^{M}\frac{x_{k}}{p_{k}}\right]^{2}\right]$
$\displaystyle=\mathbb{E}\left[\sum_{k=1}^{M}\sum_{j=1}^{M}\frac{x_{k}}{p_{k}}\frac{x_{j}}{p_{j}}\right]$
$\displaystyle=\sum_{k=1}^{M}\sum_{j=1}^{M}\mathbb{E}\left[\frac{x_{k}}{p_{k}}\frac{x_{j}}{p_{j}}\right]$
$\displaystyle=\sum_{j\neq
k}^{M}\mathbb{E}\left[\frac{x_{k}}{p_{k}}\frac{x_{j}}{p_{j}}\right]+\sum_{k=1}^{M}\mathbb{E}\left[\frac{x_{k}^{2}}{p_{k}^{2}}\right]$
$\displaystyle=\sum_{j\neq
k}^{M}\mathbb{E}\left[\frac{x_{k}}{p_{k}}\right]\mathbb{E}\left[\frac{x_{j}}{p_{j}}\right]+\sum_{k=1}^{M}\sum_{j=1}^{N}\frac{y_{j}^{2}}{\pi_{j}}$
$\displaystyle=\sum_{j\neq
k}^{M}\mathbb{E}\left[\frac{x_{k}}{p_{k}}\right]\mathbb{E}\left[\frac{x_{j}}{p_{j}}\right]+M\sum_{j=1}^{N}\frac{y_{j}^{2}}{\pi_{j}}$
$\displaystyle=M(M-1)\tau^{2}+M\sum_{j=1}^{N}\frac{y_{j}^{2}}{\pi_{j}}.$ (26)

Substituting back, we get

$\displaystyle\frac{1}{M(M-1)}\sum_{m=1}^{M}\left[\frac{x_{m}}{Np_{m}}-\frac{1}{N}\sum_{k=1}^{M}\frac{x_{k}}{Mp_{k}}\right]^{2}-\left[\frac{1}{N}\sum_{k=1}^{M}\frac{x_{k}}{Mp_{k}}\right]^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{N^{2}M(M-1)}M\sum_{j=1}^{N}\frac{y_{j}^{2}}{\pi_{j}}-$
$\displaystyle\frac{1}{N^{2}M(M-1)}\left[M(M-1)\tau^{2}+M\sum_{j=1}^{N}\frac{y_{j}^{2}}{\pi_{j}}\right]$
$\displaystyle=$
$\displaystyle\frac{1}{N^{2}(M-1)}\sum_{j=1}^{N}\frac{y_{j}^{2}}{\pi_{j}}-$
$\displaystyle\frac{1}{N^{2}(M-1)}\left[(M-1)\tau^{2}+\sum_{j=1}^{N}\frac{y_{j}^{2}}{\pi_{j}}\right]$
$\displaystyle=$ $\displaystyle-\frac{1}{N^{2}(M-1)}(M-1)\tau^{2}$
$\displaystyle=$ $\displaystyle-\frac{\tau^{2}}{N^{2}}$ $\displaystyle=$
$\displaystyle-\bar{y}^{2}.$ (27)

Combining the two estimators $T_{1}$ and $T_{2}$ we have:

$\displaystyle\mathbb{E}(T_{1}+T_{2})$
$\displaystyle=\frac{1}{N}\sum^{N}_{j=1}y_{j}^{2}-\bar{y}^{2}$
$\displaystyle=\sigma^{2}_{y}$

Hence, we have shown that the estimator of $\sigma^{2}_{y}$ is unbiased using
the sum of the estimators $T_{1}$ in Eq. 14 and $T_{2}$ in Eq. 18.

## Hierarchical models for the radon dataset

We compare seven different models of predicting the radon levels in individual
houses (indexed by $i$) by county (indexed by $j$). First we fit a pooled
model (model 1)

$\displaystyle y_{ij}$ $\displaystyle=\alpha+x_{ij}\beta+\epsilon_{ij}$
$\displaystyle\epsilon_{ij}$ $\displaystyle\sim N(0,\sigma_{y})$
$\displaystyle\alpha,\beta$ $\displaystyle\sim N(0,10)$
$\displaystyle\sigma_{y}$ $\displaystyle\sim N^{+}(0,1)\,,$

where $y_{ij}$ is the log radon level in house $i$ in county $j$, $x_{ij}$ is
the floor measurement and $\epsilon_{ij}$ is $N^{+}(0,1)$ is a truncated
Normal distribution at the positive real line. We compare this to a non-pooled
model (model 2),

$\displaystyle y_{ij}$ $\displaystyle=\alpha_{j}+x_{ij}\beta+\epsilon_{ij}$
$\displaystyle\epsilon_{ij}$ $\displaystyle\sim N(0,\sigma_{y})$
$\displaystyle\alpha_{j},\beta$ $\displaystyle\sim N(0,10)$
$\displaystyle\sigma_{y}$ $\displaystyle\sim N^{+}(0,1)\,,$

a partially pooled model (model 3),

$\displaystyle y_{ij}$ $\displaystyle=\alpha_{j}+\epsilon_{ij}$
$\displaystyle\epsilon_{ij}$ $\displaystyle\sim N(0,\sigma_{y})$
$\displaystyle\alpha_{j}$ $\displaystyle\sim N(\mu_{\alpha},\sigma_{\alpha})$
$\displaystyle\mu_{\alpha}$ $\displaystyle\sim N(0,10)$
$\displaystyle\sigma_{y},\sigma_{\alpha}$ $\displaystyle\sim N^{+}(0,1)\,,$

a variable intercept model (model 4),

$\displaystyle y_{ij}$ $\displaystyle=\alpha_{j}+x_{ij}\beta+\epsilon_{ij}$
$\displaystyle\epsilon_{ij}$ $\displaystyle\sim N(0,\sigma_{y})$
$\displaystyle\alpha_{j}$ $\displaystyle\sim N(\mu_{\alpha},\sigma_{\alpha})$
$\displaystyle\mu_{\alpha},\beta$ $\displaystyle\sim N(0,10)$
$\displaystyle\sigma_{y},\sigma_{\alpha}$ $\displaystyle\sim N^{+}(0,1)\,,$

a variable slope model (model 5),

$\displaystyle y_{ij}$ $\displaystyle=\alpha+x_{ij}\beta_{j}+\epsilon_{ij}$
$\displaystyle\epsilon_{ij}$ $\displaystyle\sim N(0,\sigma_{y})$
$\displaystyle\beta_{j}$ $\displaystyle\sim N(\mu_{\beta},\sigma_{\beta})$
$\displaystyle\mu_{\beta},\alpha$ $\displaystyle\sim N(0,10)$
$\displaystyle\sigma_{y},\sigma_{\beta}$ $\displaystyle\sim N^{+}(0,1)\,,$

a variable intercept and slope model (model 6),

$\displaystyle y_{ij}$
$\displaystyle=\alpha_{j}+x_{ij}\beta_{j}+\epsilon_{ij}$
$\displaystyle\alpha_{j}$ $\displaystyle\sim N(\mu_{\alpha},\sigma_{\alpha})$
$\displaystyle\beta_{j}$ $\displaystyle\sim N(\mu_{\beta},\sigma_{\beta})$
$\displaystyle\mu_{\alpha},\mu_{\beta}$ $\displaystyle\sim N(0,10)$
$\displaystyle\sigma_{y},\sigma_{\alpha},\sigma_{\beta}$ $\displaystyle\sim
N^{+}(0,1)\,,$

and finally a model with county level covariates and county level intercepts

$\displaystyle y_{ij}$
$\displaystyle=\alpha_{j}+x_{ij}\beta_{1}+u_{j}\beta_{2}+\epsilon_{ij}$
$\displaystyle\alpha_{j}$ $\displaystyle\sim N(\mu_{\alpha},\sigma_{\alpha})$
$\displaystyle\beta,\mu_{\alpha}$ $\displaystyle\sim N(0,10)$
$\displaystyle\sigma_{y},\sigma_{\alpha}$ $\displaystyle\sim N^{+}(0,1)\,,$

where $u_{j}$ is the log uranium level in the county. The Stan code used can
be found below.

## Stan models

###  4.4 Linear regression model

⬇

data {

int <lower=0> N;

int <lower=0> D;

matrix [N,D] x ;

vector [N] y;

}

parameters {

vector [D] b;

real <lower=0> sigma;

}

model {

target += normal_lpdf(y | x * b, sigma);

target += normal_lpdf(b | 0, 1);

}

generated quantities{

real log_joint_density_unconstrained;

vector[N] log_lik;

// Compute the log likelihoods for loo

for (n in 1:N) {

log_lik[n] =

normal_lpdf(y[n] | x[n,] * b, sigma);

}

}

###  4.5 Radon pooled model (1)

⬇

data {

int<lower=0> N;

vector[N] x;

vector[N] y;

int<lower=0,upper=1> holdout[N];

}

parameters {

vector[2] beta;

real<lower=0> sigma_y;

}

model {

vector[N] mu;

// priors

sigma_y ~ normal(0,1);

beta ~ normal(0,10);

// likelihood

mu = beta[1] + beta[2] * x;

for(n in 1:N){

if(holdout[n] == 0){

target +=

normal_lpdf(y[n]|mu[n],sigma_y);

}

}

}

###  4.6 Radon pooled model (2)

⬇

data {

int<lower=0> N;

int<lower=0> J;

int<lower=1,upper=J> county[N];

vector[N] x;

vector[N] y;

int<lower=0,upper=1> holdout[N];

}

parameters {

vector[J] a;

real beta;

real<lower=0> sigma_y;

}

model {

vector[N] mu;

// Prior

sigma_y ~ normal(0,1);

a ~ normal(0,10);

// Likelihood

for(n in 1:N){

mu[n] = beta*x[n] + a[county[n]];

if(holdout[n] == 0){

target +=

normal_lpdf(y[n]|mu[n],sigma_y);

}

}

}

###  4.7 Radon partially pooled model (3)

⬇

data {

int<lower=0> N;

int<lower=0> J;

int<lower=1,upper=J> county[N];

vector[N] y;

int<lower=0,upper=1> holdout[N];

}

parameters {

vector[J] a;

real mu_a;

real<lower=0> sigma_a;

real<lower=0> sigma_y;

}

model {

vector[N] mu;

// priors

sigma_y ~ normal(0,1);

sigma_a ~ normal(0,1);

mu_a ~ normal(0,10);

// likelihood

a ~ normal (mu_a, sigma_a);

for(n in 1:N){

mu[n] = a[county[n]];

if(holdout[n] == 0){

target +=

normal_lpdf(y[n]|mu[n],sigma_y);

}

}

}

###  4.8 Variable intercept model (4)

⬇

data {

int<lower=0> J;

int<lower=0> N;

int<lower=1,upper=J> county[N];

vector[N] x;

vector[N] y;

int<lower=0,upper=1> holdout[N];

}

parameters {

vector[J] a;

real beta;

real mu_a;

real<lower=0> sigma_a;

real<lower=0> sigma_y;

}

model {

vector[N] mu;

// Prior

sigma_y ~ normal(0,1);

sigma_a ~ normal(0,1);

mu_a ~ normal(0,10);

beta ~ normal(0,10);

a ~ normal (mu_a, sigma_a);

for(n in 1:N){

mu[n] = a[county[n]] + x[n]*beta;

if(holdout[n] == 0){

target +=

normal_lpdf(y[n]|mu[n],sigma_y);

}

}

}

###  4.9 Variable slope model (5)

⬇

data {

int<lower=0> J;

int<lower=0> N;

int<lower=1,upper=J> county[N];

vector[N] x;

vector[N] y;

int<lower=0,upper=1> holdout[N];

}

parameters {

real a;

vector[J] beta;

real mu_beta;

real<lower=0> sigma_beta;

real<lower=0> sigma_y;

}

model {

vector[N] mu;

// Prior

a ~ normal(0,10);

sigma_y ~ normal(0,1);

sigma_beta ~ normal(0,1);

mu_beta ~ normal(0,10);

beta ~ normal(mu_beta,sigma_beta);

for(n in 1:N){

mu[n] = a + x[n] * beta[county[n]];

if(holdout[n] == 0){

target +=

normal_lpdf(y[n]|mu[n],sigma_y);

}

}

}

###  4.10 Variable intercept and slope model (6)

⬇

data {

int<lower=0> N;

int<lower=0> J;

vector[N] y;

vector[N] x;

int county[N];

int<lower=0,upper=1> holdout[N];

}

parameters {

real<lower=0> sigma_y;

real<lower=0> sigma_a;

real<lower=0> sigma_beta;

vector[J] a;

vector[J] beta;

real mu_a;

real mu_beta;

}

model {

vector[N] mu;

// Prior

sigma_y ~ normal(0,1);

sigma_beta ~ normal(0,1);

sigma_a ~ normal(0,1);

mu_a ~ normal(0,10);

mu_beta ~ normal(0,10);

a ~ normal(mu_a, sigma_a);

beta ~ normal(mu_beta, sigma_beta);

for(n in 1:N){

mu[n] = a[county[n]] + x[n]*beta[county[n]];

if(holdout[n] == 0){

target +=

normal_lpdf(y[n]|mu[n],sigma_y);

}

}

}

###  4.11 Hierarchical intercept model (7)

⬇

data {

int<lower=0> J;

int<lower=0> N;

int<lower=1,upper=J> county[N];

vector[N] u;

vector[N] x;

vector[N] y;

int<lower=0,upper=1> holdout[N];

}

parameters {

vector[J] a;

vector[2] beta;

real mu_a;

real<lower=0> sigma_a;

real<lower=0> sigma_y;

}

transformed parameters {

}

model {

vector[N] mu;

vector[N] m;

sigma_a ~ normal(0, 1);

sigma_y ~ normal(0, 1);

mu_a ~ normal(0, 10);

beta ~ normal(0, 10);

a ~ normal(mu_a, sigma_a);

for(n in 1:N){

m[n] = a[county[n]] + u[n] * beta[1];

mu[n] = m[n] + x[n] * beta[2];

if(holdout[n] == 0){

target += normal_lpdf(y[n] | mu[n], sigma_y);

}

}

}