Unbiased estimation of log normalizing constants with applications to Bayesian cross-validation

Setting

In this article we propose a new estimator of $Z$, which combines unbiased Markov chain Monte Carlo (Jacob et al., 2017) with the path sampling identity (Gelman and Meng, 1998; see also Chapter 5 of Chen et al., 2000), also known as thermodynamic integration (Kirkwood, 1935; Neal, 2005; Calderhead and Girolami, 2009). The specificity of the proposed estimator is its unbiasedness for the logarithm of $Z$, i.e. the expectation of the proposed estimator is exactly $\log Z$. Existing estimators based on Markov chain Monte Carlo (Chen et al., 1997) are only asymptotically unbiased, while existing estimators based on annealed importance samplers (Neal, 2001) and sequential Monte Carlo samplers (Del Moral et al., 2006) are unbiased for $Z$ and not for $\log Z$.

Leveraging unbiasedness for $\log Z$, we consider a Bayesian cross-validation (CV) criterion based on the logarithmic scoring rule (Alqallaf and Gustafson, 2001; Bornn et al., 2010; Vehtari et al., 2017, e.g). In cross-validation, one randomly splits the available data into training and validation, then the posterior distribution given the training data is numerically approximated, and finally the predictive performance on the validation data is assessed e.g. with the logarithmic scoring rule (Parry et al., 2012). We propose an estimator that is directly unbiased for these Bayesian cross-validation objectives, which can be averaged over independent copies to obtain consistent estimators and asymptotically exact confidence intervals from the central limit theorem for i. i. d. variables.

The rest of the document is structured as follows. Section 2 introduces the proposed estimators, and their tuning parameters are discussed. Numerical experiments in simple examples can be found in Section 3. Section 4 discusses our findings and future directions. The code to reproduce the experiments of the article is available at https://github.com/pierrejacob/unbiasedpathsampling.

Proposed estimators

We propose an unbiased estimator of $\log Z$ in Section 2.1, and obtain an unbiased estimator of a Bayesian cross-validation criterion in Section 2.2. Our implementation relies on the unbiased MCMC estimators of Jacob et al., 2017, which are briefly reviewed in Section 2.3, while Section 2.4 discusses tuning choices.

The thermodynamic integration or path sampling identity relies on the following interchange between differentiation and integration (Kirkwood, 1935),

By introducing an arbitrary density $\lambda\mapsto q(\lambda)$, strictly positive on $(0,1)$, we obtain the path sampling identity:

This is useful if we can approximate integrals with respect to $\lambda$ by Monte Carlo or numerical integration, and if we can approximate the inside expectation $\E_{\lambda}[\nabla_{\lambda}U_{\lambda}(X)]$ by Markov chain Monte Carlo (MCMC, Robert and Casella, 2004), for instance.

Here we denote by $E(\lambda)=-\E_{\lambda}[\nabla_{\lambda}U_{\lambda}(X)]$ the inner expectation in (3), and we introduce $\hat{E}(\lambda)$, an unbiased estimator of $E(\lambda)$ that we can generate for any $\lambda\in$; we defer the construction of such estimators to Section 2.3. We can then define an estimator of $r_{01}=\log(Z_{1}/Z_{0})$ with the following procedure.

Draw $\lambda\sim q(\mathop{d\lambda})$, a distribution supported on $$.

Given $\lambda$, generate a variable $\hat{E}(\lambda)$ with expectation $E(\lambda)=-\E_{\lambda}[\nabla_{\lambda}U_{\lambda}(X)]$.

Return $\hat{r}_{01}=\hat{E}(\lambda)/q(\lambda)$.

The random variable $\hat{r}_{01}$ has expectation $r_{01}$ by the law of iterated expectations, and we refer to it as an unbiased path sampling estimator (UPS). Note that sequential Monte Carlo samplers and related methods (Del Moral et al., 2006) would provide unbiased estimators of $Z$ and not of $\log Z$. Thus these estimators will not be unbiased for $r_{01}$. We will now see that the lack of bias on the logarithmic scale can be exploited to propose new estimators of Bayesian cross-validation criteria.

2 Unbiased Bayesian cross-validation

A number of articles discuss the computational difficulties associated with Bayesian cross-validation, e.g. Alqallaf and Gustafson, 2001; Bhattacharya and Haslett, 2007; Bornn et al., 2010; Lamnisos et al., 2012; McVinish et al., 2013; Vehtari et al., 2017. We first define the object of interest, before presenting our estimator. Let $x$ denote an unknown parameter with prior density $\prob(x)$, and let $y_{1:n}=\{y_{1},\dotsc,y_{n}\}$ denote the data composed of $n$ units. The likelihood function is denoted by $x\mapsto\prob(y_{1:n}\mid x)$. Cross-validation consists in randomly splitting $y_{1:n}$ into $T$ and $V$, where $T$ stands for training and $V$ for validation. The sets $T,V$ form a partition of $y_{1:n}$, $T\cap V=\emptyset$ and $T\cup V=y_{1:n}$. Denote by $n_{T}$ and $n_{V}$ the numbers of elements in $T$ and $V$; for instance, if $n_{T}=n-1$, the procedure is termed “leave-one-out” cross-validation. Given a split of the data $y_{1:n}=(T,V)$, we introduce a measure of accuracy in predicting $V$ using the training data $T$. A typical choice is the logarithmic score $-\log\prob(V\mid T)$ (see Parry et al., 2012, for a discussion on the choice of scoring rule) where $\prob(V\mid T)=\int\prob(V\mid T,x)\prob(x\mid T)\mathop{dx}$ is the posterior predictive density given $T$ and evaluated on $V$. Note that $\prob(V\mid T,x)$ simplifies to $\prob(V\mid x)$ if the data are modeled as conditionally independent given $x$. The cross-validation objective, “CV” below, is defined as an average over all splits $(T,V)$ of size $(n_{T},n_{V})$,

Given a split $T,V$, we can estimate $\log\prob(V\mid T)$ using the path sampling identity and the unbiased estimators of the previous section. Indeed, that quantity is a log-ratio of the normalizing constants $\prob(T)$ and $\prob(T,V)$. By introducing the path

This motivates the following strategy: sample a split $T,V$ uniformly from $\mathcal{S}$, and then obtain an unbiased estimator of $-\log\prob(V\mid T)$ given $T,V$. The resulting estimator is directly unbiased for CV in (4), by the law of iterated expectations. We summarize the procedure below.

Sample index sets $T,V$ uniformly at random over $\mathcal{S}$, the set of partitions of $\{1,\dotsc,n\}$ into a set of size $n_{T}$ and a set of size $n_{V}=n-n_{T}$.

Note how the lack of bias on the logarithmic scale is important for the above procedure to produce an unbiased estimator of CV. We could also extend the above procedure to allow for non-uniform sampling of the partitions from $\mathcal{S}$.

3 Reminders on unbiased MCMC

4 Tuning choices

A number of choices have to be made for the proposed estimators to be operational. The first choice is that of a path of distributions. There are generic choices such as the geometric path, and choices motivated by algorithmic considerations on a case-by-case basis. We will discuss the choice of paths through examples, in Section 3.

Given a path of distributions $(\pi_{\lambda})$, algorithms approximating expectations $\E_{\lambda}$ with respect to $\pi_{\lambda}$ typically involve tuning parameters. We describe the tuning of unbiased MCMC in Section 2.4.1. Then we discuss choices of distribution $q(\mathop{d\lambda})$ in Section 2.4.2.

The unbiased MCMC estimators described in Section Section 2.3 require the specification of a Markov kernel $P$, a coupled kernel $\bar{P}$, and an initial distribution for the chains. Specifying these objects is typically difficult, but not specific to the setting of normalizing constant estimation. Therefore we defer to the large literature on MCMC algorithms (Robert and Casella, 2004; Brooks et al., 2011), as well as the relevant discussions in Jacob et al., 2017 in the context of unbiased MCMC. Ultimately we will care about the expected cost and the variance of the proposed unbiased estimators, in order to maximize the efficiency of the proposed estimators, as discussed in the next section.

4.2 Tuning of the distribution q⁡(d​λ)q(\mathop{d\lambda})

If $C_{\lambda}$ is constant over $\lambda$, then the solution of the above minimization problem is given by $\lambda\mapsto q^{\star}(\lambda)\propto\sqrt{m_{2}(\lambda)}$. In Gelman and Meng, 1998 that solution is given, and then the verification that this is indeed a solution is done via Cauchy-Schwarz. Here we provide an informal derivation of the solution, in the case where $C_{\lambda}$ is constant over $\lambda$. We write the function to minimize as $\int m_{2}(\lambda)/q(\lambda)\mathop{d\lambda}$, and introduce the Lagrangian

We would like to differentiate with respect to $q$ and set the derivative to zero. Introduce the directional derivative $\frac{\mathop{d}}{\mathop{d\varepsilon}}(q(\lambda)+\varepsilon v(\lambda))$ where $v(\lambda)$ is a function. Replacing $q(\lambda)$ by $q(\lambda)+\varepsilon v(\lambda)$ and differentiating with respect to $\varepsilon$ in the Lagrangian yields

Setting $\varepsilon$ to zero yields $\int\{-m_{2}(\lambda)/q(\lambda)^{2}+\xi\}v(\lambda)\mathop{d\lambda}$, and trying to set that expression to zero simultaneously for every choice of $v$, we obtain $-m_{2}(\lambda)/(q(\lambda)^{2})+\xi=0$, i.e. $q(\lambda)\propto\sqrt{m_{2}(\lambda)}$. This gives the candidate solution.

4.3 Proposed tuning procedure

We now combine the above sections into practical guidelines for the proposed estimators.

For $\lambda$ in a grid of $L+1$ values $0=\lambda^{}\leq\ldots\leq\lambda^{[L]}=1$, construct and tune an unbiased MCMC (initial distribution, Markov kernel $P$, and coupled kernel $\bar{P}$) targeting $\pi_{\lambda}$, and draw independent samples of the associated meeting times $\tau_{\lambda}$.

Use these estimates to define a distribution $q(\mathop{d\lambda})$, such that $q(\lambda)$ is approximately proportional to $\sqrt{m_{2}(\lambda)}$ for all $\lambda$ in $$.

We describe a concrete way of performing step 5, for completeness. Given a grid of values $0=\lambda^{}\leq\ldots\leq\lambda^{[L]}=1$ and associated estimates $(\hat{m}_{2}(\lambda^{[l]}))^{1/2}$ of $(m_{2}(\lambda^{[l]}))^{1/2}$ for $l\in\{0,\dotsc,L\}$ obtained in step 4, we can define a distribution $q(\mathop{d\lambda})$ that is piecewise uniform on the intervals $[\lambda^{[l]},\lambda^{[l+1]}]$, and such that

Sampling from such a distribution can be done in order $\mathcal{O}(\log L)$ operations, by first selecting an interval $[\lambda^{[l]},\lambda^{[l+1]}]$ with probability $\int_{\lambda^{[l]}}^{\lambda^{[l+1]}}q(\mathop{d\lambda})$, and then sampling uniformly from that interval.

After the preliminary phase described in the five steps above, the generation of estimators $\hat{r}_{01}$ can proceed as follows. First, $\lambda$ is drawn from $q(\mathop{d\lambda})$ obtained in step 5 above. We then find the nearest value $\lambda^{[l]}$ in the grid, with index $l\in\{0,\dotsc,L\}$. We can look up tuning parameters corresponding to $\lambda^{[l]}$ for the unbiased MCMC estimators, stored during step 2 above, and the values of $k_{\lambda^{[l]}}$ and $m_{\lambda^{[l]}}$ stored during step 3 above. Using these tuning values we can generate an unbiased estimator $\hat{E}(\lambda)$ of $E(\lambda)=-\E_{\lambda}[\nabla_{\lambda}U_{\lambda}(X)]$. The estimator $\hat{r}_{01}=\hat{E}(\lambda)/q(\lambda)$ is finally returned.

Numerical experiments

The numerical experiments are structured as follows. Section 3.1 contains toy examples of unbiased path sampling estimators. Section 3.2 considers logistic regressions with different choices of paths and of unbiased MCMC estimators, and an example taken from Epifani et al., 2008; Vehtari et al., 2017. Section 3.3 considers linear regressions with examples taken from Alqallaf and Gustafson, 2001; Peruggia, 1997; Vehtari et al., 2017. Throughout the experiments, 95% confidence intervals for an estimand $\mu$ are obtained as $\bar{\mu}\pm 1.96s/\sqrt{M}$, where $\bar{\mu}$ is the mean of $M$ independent unbiased estimators of $\mu$ and $s$ is their sample standard deviation. These confidence intervals are justified asymptotically as $M\to\infty$ by the central limit theorem for i. i. d. random variables, provided that the variance of the unbiased estimators is finite. On parallel machines and under budget constraints, valid confidence intervals can be constructed following Glynn and Heidelberger, 1991; see also related remarks in Jacob et al., 2017.

We start with a grid of values of $\lambda$: $\lambda^{[l]}=l/L$ for $l\in\{0,\dotsc,L\}$, with $L=10$. For each $\lambda^{[l]}$, we run coupled MH chains until they meet, 100 times independently. We obtain a distribution of meeting times $\tau$ for each $\lambda$, represented on 1(a). The overlaid full line represents the $99\%$ quantiles, which we denote by $k_{0},\dotsc,k_{L}$. We also compute the average meeting times for each $\lambda^{[l]}$, which we denote $\bar{\tau}_{0},\dotsc,\bar{\tau}_{L}$. We then define

Given values of $k_{l}$ and $m_{l}$, for each $\lambda^{[l]}$ in the grid of $L+1$ values defined above, we approximate the first and second moments of $\pi_{\lambda}$ with $100$ independent estimators. We use these moments to redefine the initial distribution of the Markov chains, which we set to a Normal distribution adapted to $\pi_{\lambda}$, and to tune the proposal standard deviation, which we set to be the estimated standard deviation of $\pi_{\lambda}$. At this point we could sample meeting times again and choose new values for $k_{\lambda}$ and $m_{\lambda}$, but we omit this here. Next, we estimate $\sqrt{m_{2}(\lambda)}$ for each $\lambda^{[l]}$ in the grid, and define $q(\mathop{d\lambda})$ accordingly, following step 5 in Section 2.4. The estimates of $\sqrt{m_{2}(\lambda)}$ are shown in 1(b). This completes the tuning phase, and we can now generate unbiased estimators $\hat{r}_{01}=\hat{E}(\lambda)/q(\lambda)$ of $r_{01}=\log(Z_{1}/Z_{0})$. We show these estimates against $\lambda$ in 1(c). These are generated $5,000$ times independently. Concretely, they yield the confidence interval $[-0.11,0.12]$ for the estimand $r_{01}=0$ at level $95\%$.

1.2 Double-well example

We perform similar experiments on a path of two-dimensional distributions linking the potential $U_{0}:x\mapsto(x_{1}+2)^{2}+(x_{2}^{2}/2)$, corresponding to a Normal distribution centered at $(-2,0)$ and with diagonal variances $(1/2,1)$, to the potential $U_{1}:x\mapsto(1/10)(((x_{1}-1)^{2}-x_{2}^{2})^{2}+10(x_{1}^{2}-5)^{2}+(x_{1}+x_{2})^{4}+(x_{1}-x_{2})^{4})$. The latter is a double-well potential, with modes around $(-2,0)$ and $(2,0)$. By numerical integration we find $\log(Z_{1}/Z_{0})$ to be approximately $-6.9$. We introduce the geometric path $U_{\lambda}(x)=(1-\lambda)U_{0}(x)+\lambda U_{1}(x)$. For each $\lambda$, we start chains from a Normal centered at $(-2,-2)$ and with covariance matrix $\mathbf{I}_{2}$, the identity matrix of size $2\times 2$. We consider random walk MH schemes with Normal proposal, with covariance $2\mathbf{I}_{2}$; the coupled version relies on maximal couplings of the proposals, as in the previous section.

We draw $1,000$ meeting times independently, for $\lambda^{[l]}=l/L$ with $l\in\{0,\dotsc,L\}$ and $L=10$. The distributions are shown in violin plots in 2(a). We observe much larger meeting times for $\lambda$ close to one, which corresponds to the MH chains struggling to explore both modes of the double-well potential. We thus conservatively set $k_{l}$ to be twice the $99\%$ quantiles of the meeting times, instead of the quantiles themselves.

We follow the same heuristics as in Section 3.1.1 for the choice of $m_{l}$. Without modifying the initial distribution nor the proposal distribution of the MH chains, we estimate $\sqrt{m_{2}(\lambda)}$ for each $\lambda^{[l]}$ in the grid, based on 100 independent copies, and define $q(\mathop{d\lambda})$ following again step 5 in Section 2.4. The estimates of $\sqrt{m_{2}(\lambda)}$ are shown in 2(b). Finally we generate unbiased estimators $\hat{r}_{01}=\hat{E}(\lambda)/q(\lambda)$ of $r_{01}=\log(Z_{1}/Z_{0})$, and represent these estimates against $\lambda$ in 2(c). These are generated $1,000$ times independently and result in the $95\%$ confidence interval $[-7.55,-6.37]$ for the estimand $r_{01}\approx-6.9$.

2 Logistic regression

With basic manipulations this is equivalent to the following simpler form

The Pólya-Gamma Gibbs (PGG) sampler (Polson et al., 2013; Choi and Hobert, 2013) is a Gibbs sampler that targets $\prob(\beta\mid Y)$ through the introduction of auxiliary variables $W$. First, we recall that the Pólya-Gamma distribution with parameters $(1,c)$, denoted by $\text{PG}(1,c)$, has a density $x\mapsto\pg(x;c)$ defined for all $c\geq 0$, $x>0$ as

Introduce $n$ auxiliary variables $W=(W_{1},\dotsc,W_{n})$, independent of each other given $\beta$, such that $W_{i}$ follows PG$(1,\left\lvert d_{i}^{\intercal}\beta\right\rvert)$ for all $1\leq i\leq n$. An extended target distribution is defined as $\prob(\beta,\omega\mid Y)\propto\prob(\beta\mid Y)g(\omega\mid\beta)$, where $\omega$ denotes a realization of $W$, and $g(\omega\mid\beta)=\prod_{i=1}^{n}\pg(\omega_{i};\left\lvert d_{i}^{\intercal}\beta\right\rvert)$. The appeal of this extension is that we can write the target as

and therefore the conditional of $\beta$ given $\omega,Y$ simplifies to

given $(\beta^{(t)},W^{(t)})$, draw $\beta^{(t+1)}\sim\normal\left(\mu(W^{(t)}),\Sigma(W^{(t)})\right)$,

draw $W_{i}^{(t+1)}\sim\text{PG}(1,\left\lvert d_{i}^{\intercal}\beta^{(t+1)}\right\rvert)$, independently for all $i\in\left\{1,\dotsc,n\right\}$.

In the experiments below, we initialize the chains from the prior distribution $\mathcal{N}(b,B)$.

We first remark that the above reasoning holds when replacing the covariates $D$ by $\lambda D$ for any $\lambda\in$. This corresponds to the likelihood $\beta\mapsto\prod_{i=1}^{n}\exp(\lambda d_{i}^{\intercal}\beta y_{i})/(1+\exp(\lambda d_{i}^{\intercal}\beta))$, for $\lambda\in$. In the case $\lambda=0$, the likelihood is equal to $2^{-n}$ for all $\beta$, while with $\lambda=1$, we retrieve the original likelihood. For all $\lambda$, we can introduce Pólya-Gamma variables $W_{i}$ following PG$(1,\left\lvert\lambda d_{i}^{\intercal}\beta\right\rvert)$ for all $1\leq i\leq n$, and obtain a corresponding PGG sampler.

This enables normalizing constant estimators for the logistic regression model with little tuning, since the PGG sampler itself has no tuning parameters. Here, for all $\lambda,\beta$, we define

We consider a synthetic data set with $n=1000$ rows and $p=7$ columns. The covariates are generated from a standard Normal distribution and the outcome is generated from the model with $\beta^{\star}=(0,0.1,0.2,0.3,0.4,0.5,0.6)$. The prior mean $b$ is set to zero and the covariance $B$ to a diagonal matrix with entries equal to $10$. We start by gridding the interval $$, and for each value$\lambda^{[l]}=l/L$with$L=10$, we set$k$as the$99\%$quantile of the meeting times for the coupled PGG sampler, based on 1000 independent runs. We set$m$as in the previous sections, to make the average cost approximately constant over$\lambda$. Next we estimate the second moments of$\hat{E}(\lambda)$on the grid of values of$\lambda$, we design a proposal$q(\mathop{d\lambda})$ following step 5 in Section 2.4, and obtain the estimates of 3(a).

We obtain the $1,000$ independent estimators $\hat{E}(\lambda)/q(\lambda)$ shown in 3(b), leading to a $95\%$ confidence interval of $$on$r_{01}=\log(Z_{1}/Z_{0})$. The actual value is found to be close to$70$using importance sampling based on a Laplace approximation to the posterior, accurate in the present example (see below, and also related discussions in Bardenet et al., 2017). We can see from 3(b) that the estimates take very large values for$\lambda$close to zero. This suggests that, instead of choosing an equispaced grid of values of$\lambda$on$$when designing$q(\mathop{d\lambda})$, we could aim at a higher resolution towards the left end of the interval$$.

Therefore we consider a grid of values of $\lambda$ equispaced on the logarithmic scale: $\lambda^{[l]}=\exp(-L+l)$ for $l=0,\dotsc,L$ with $L=10$. Going through the exact same tuning steps, we obtain the $1,000$ estimators of 3(c), leading to the narrower confidence interval $[64.7,74.7]$ at level $95\%$ (with a width of 10 instead of 36 for the previous one). This illustrates the potential gains obtained by carefully choosing the distribution $q(\mathop{d\lambda})$.

We conclude this section by noting that more dramatic gains can be obtained by changing the path of distributions. In the context of logistic regression with $n\gg p$, the Laplace approximation of the posterior, defined as $\normal(\hat{\beta}_{\text{MLE}},\hat{V})$ where $\hat{\beta}_{\text{MLE}}$ is the maximum likelihood estimator and $\hat{V}$ is the inverse of minus the Hessian of the log-likelihood evaluated at $\hat{\beta}_{\text{MLE}}$, seems to be very accurate. We thus introduce a geometric path $(\pi_{\lambda})$ between the Laplace approximation and the posterior distribution. We use a random walk MH algorithm to target $\pi_{\lambda}$ for all $\lambda\in$, with proposal covariance matrix equal to $\hat{V}/p$ where the dimension $p$ is equal to $7$. To couple the MH algorithms, we use strategy that combines reflection and maximal couplings, as described in Jacob et al., 2017. The initial distribution of the chains is chosen to be the Laplace approximation. For $\lambda=0$, we obtain $k=127$ as the $99\%$ quantile of the meeting times, and we set $m=5k$. We use these values of $k$ and $m$ for all $\lambda$, and we choose $q(\mathop{d\lambda})$ to be uniform on $$. With$100$independent estimators we obtain a confidence interval of$[70.24,70.26]$at$95\%$for$\log Z_{1}+n\log 2$. This is orders of magnitude narrower than the previous intervals, for a smaller computational cost. The choice of paths can thus play a critical role in the efficiency of the proposed estimators, and approximations of the posterior distribution can be used to construct such paths.

2.2 Cross-validation

We now consider the approximation of CV in (4). We consider a leave-one-out criterion, with $n_{T}=n-1$ and $n_{V}=1$. We thus construct paths between the posterior given the training data $T$, with normalizing constant $\prob(T)$, and the posterior given all the data $(T,V)$, with normalizing constant $\prob(T,V)$.

Our first path follows the reasoning of the previous section: we can multiply the covariates in the validation set by $\lambda\in$ to preserve the original structure of the likelihood and thus to enable a similar PGG sampler. The unnormalized densities are then

Again we see that this is essentially a linear function of $\beta$ and thus its moments under $\pi_{\lambda}$ are finite for all $\lambda$.

To tune the procedure, we obtain meeting times for the coupled PGG sampler based on the full data set, and choose $k$ as a $99\%$ quantile (here equal to $8$), and $m=5k=40$. Recall that the PGG sampler itself has no tuning parameters. Then, drawing a validation set at random $1,000$ time independently, generating $\lambda$ uniformly on $$and obtaining the associated estimator$\hat{E}(\lambda)$, we obtain unbiased estimators of CV in (4). We plot a histogram of these estimators in 4(a). A$95\%$confidence interval for the CV objective is obtained as$[-0.62,-0.56]$.

Alternatively, we introduce a geometric path between the posterior given $T$ and given $T,V$, which corresponds to the unnormalized densities

For this path, we use random walk MH as in the previous section, with initial distribution and proposal covariance tuned using a Laplace approximation of the posterior distribution. We obtain a $99\%$ quantile of meetings at $k=135$ and set $m=5k$. Over $1,000$ independent experiments we obtain unbiased estimators of the CV objective shown in 4(b). The associated $95\%$ confidence interval for CV is $[-0.60,-0.56]$. Thus, this second approach appears to be marginally more efficient than the first one; the cost comparison is made slightly difficult by the fact that PGG and MH have different costs per iteration.

2.3 Leukemia survival data

We follow Vehtari et al., 2017 and consider the leukemia data presented in Feigl and Zelen, 1965 and used as illustration in Epifani et al., 2008. We use the data formatted as in the package BGPhazard, see Garcıa-Bueno and Nieto-Barajas, 2016. The outcome is taken to be one if the survival time (column time of leukemiaFZ) is larger or equal to $50$ weeks, zero otherwise, and the two covariates are the columns wbc and AG, corresponding to counts of white blood cells and the outcome of a test related to white blood cell characteristics. There are 31 patients in the sample, so $n=31$, and we consider leave-one-out cross-validation, i.e. $n_{T}=n-1$ and $n_{V}=1$.

We introduce a path of distributions amenable to PGG sampling, as in the previous sections. Sampling uniformly the index of the observation to be left out, then sampling $\lambda$ uniformly in $$, and finally running coupled PGG chains targeting$\pi_{\lambda}$, we record the meeting times. We do so 1,000 times independently, and show the results as a function of the index of the observation left out in 5(a).

Based on this plot we select $k=100$, conservatively, and $m=5k=500$ for all runs. We then generate $10,000$ unbiased estimators of CV. We plot the estimators against the index of the left-out observation in 5(b), and we note that the values are very different for one particular index, here equal to 17. In 5(c) we plot a histogram of the estimates of the CV objective, putting all the indices together. From these estimates we obtain a $95\%$ confidence interval $[-0.72,-0.66]$ for the leave-one-out CV objective. Thus we see that the proposed estimators can have a larger variance for certain splits of the data compared to others. Investigating further the behavior of the estimators for certain splits, one might be able to reduce the variance, for instance by tuning the proposal distribution $q(\mathop{d\lambda})$, or by changing the path. Our estimators of CV might also be considered satisfactory as they stand. In any case, they do not suffer from infinite variance issues typically associated with importance sampling, when using a proposal distribution that has lighter tails than the target distribution.

3 Linear regressions

We next consider linear regressions, which have been used to illustrate Bayesian cross-validation e.g. in Alqallaf and Gustafson, 2001; Peruggia, 1997; Vehtari et al., 2017.

The first example is taken from Alqallaf and Gustafson, 2001. The data comprise of $n=62$ observations, each corresponding to an animal (arctic fox, owl monkey, etc). For each animal, the data set contains the body weight and the brain weight. The covariate $d_{i}$ of animal $i$ is a vector, with first entry equal to $1$ and second entry equal to the logarithm of body weight, while the outcome $y_{i}$ is the logarithm of brain weight. As before we write $Y$ for the vector of outcomes and $D$ for the matrix of covariates, on which we condition throughout. The model is given by

To get the conditional distribution of $\beta$ given $\sigma^{2}$ under the posterior distribution, note that

where $\hat{\beta}=(D^{\intercal}D)^{-1}D^{\intercal}Y$. Thus, the conditional distribution is Normal with mean $\hat{\beta}$ and covariance matrix $\sigma^{2}(D^{\intercal}D)^{-1}$. The distribution of $\sigma^{2}$ given $\beta$ is inverse Gamma, where recall that

Then $\sigma^{2}$ given $\beta$ is inverse Gamma with $a=n/2$ and $b=\left\lVert Y-D\beta\right\rVert^{2}/2$. Coupling this algorithm can be done by maximal coupling of each of the conditional update of a Gibbs sampler.

In Alqallaf and Gustafson, 2001, the predictive performance in this example is measured by the mean squared error, defined conditional on a split as

We implement this procedure and draw $1,000$ independent coupled chains. We observe meeting times between $1$ and $5$. Thus we set $k=10$, $m=25$, and draw $1,000$ independent unbiased estimators of CV. We obtain a $95\%$ confidence interval of $[32.79,33.03]$, and standard error of $0.06$. By comparison, Alqallaf and Gustafson, 2001 use 200 splits, and run 125 iterations of MCMC for each split, discarding the first 100. The total number of Gibbs iterations performed is approximately the same, and Alqallaf and Gustafson, 2001 obtain standard errors that are similar. An advantage of our method is in its simplicity: if we want more precise results, we simply generate more independent estimators.

We now consider the criterion $-\log\prob(Y_{V}\mid Y_{T})$, instead of the point-prediction mean squared error as above. The sequence of distributions defined in (5) is still amenable to a Gibbs sampling strategy and we need to work out the conditional distributions. The joint posterior density is

so that $\beta$ given the rest is $\normal(\mu_{\lambda},\Lambda_{\lambda}^{-1})$. On the other hand $\sigma^{2}$ given the rest is inverse Gamma $(a,b)$ with

3.2 Stack loss data

We consider the stack loss data example, which was considered in Peruggia, 1997; Vehtari et al., 2017. In the former article, it is shown that importance sampling from the posterior given all the data to the posterior leaving one data point out can lead to infinite variance estimators. Here we use the stackloss data set of (R Core Team, 2015), with the outcome set to be the column stack.loss, and the covariates Air.Flow, Water.Temp, Acid.Conc., and a column of ones. The data are shown in 6(a). We consider leave-one-out cross-validation, with $n_{T}=n-1=20$ here. For simplicity we use the same model as in the previous section, with a flat prior on $\beta$ given $\sigma^{2}$, instead of the proper prior given in Peruggia, 1997.

Using the coupled Gibbs sampler described in the previous section, we find meeting times to be less than $10$ with large probability, thus we set $k=10$ and $m=25$. We obtain the CV estimators shown in 6(b), based on $10,000$ independent replicates, plotted against the index of the left-out observation. As in Section 3.2.3, we can see that the variance of the CV estimators varies across the different ways of partitioning the data into training and validation sets. These CV estimators yield the $95\%$ confidence interval $[2.78,2.82]$.

Discussion

Further work will be needed to compare the proposed estimators with state-of-the-art methods such as sequential Monte Carlo samplers for normalizing constant estimation (Lee and Whiteley, 2015; Zhou et al., 2016; Andrieu et al., 2016, e.g.), with alternative approaches such as the ones described in Chen et al., 1997; Johnson, 1999; Neal, 2005; Salomone et al., 2018 and references therein, and with the different existing approaches for Bayesian cross-validation (Alqallaf and Gustafson, 2001; Bornn et al., 2010; Vehtari et al., 2017, e.g.).

Our estimators combine the path sampling identity with unbiased estimators of intractable integrals. As such, they are expected to break if either path sampling or the unbiased estimators break. Path sampling can give poor results if the path of distributions is ill-chosen, thus the design of these paths remains crucial. We have seen in Section 3.2 that different paths can give orders of magnitude differences in efficiencies. We have also seen that the paths can benefit from approximations of the posterior distribution, such as Laplace approximations. Mixtures of distributions fitted on MCMC samples or variational approximations could also be considered. Conditional on a path, the choice of distribution $q(\mathop{d\lambda})$ is also important and can be guided by preliminary runs. Unbiased MCMC estimators themselves break either if the underlying MCMC algorithms mix poorly, or if the coupling strategy is ineffective; we defer to Jacob et al., 2017 for related discussions, and to Heng and Jacob, 2018 for the case of Hamiltonian Monte Carlo algorithms.

We note that the path sampling identity (3) is an instance of a nested Monte Carlo (MC) problem, as defined and discussed in Rainforth et al., 2016. The target of nested MC is an expectation $I$ of the form

where the functions $f$, $\phi$ and the joint distribution of $(x,\lambda)$ are problem-dependent choices. In the case of path sampling, we obtain $I=r_{01}$ by choosing:

In this case $f(\lambda,E\left\lparen\lambda\right\rparen)$ is linear in its second argument, thus, given $\lambda$, unbiased estimators of $E\left\lparen\lambda\right\rparen$ directly translate into unbiased estimators of $f(\lambda,E\left\lparen\lambda\right\rparen)$. We remark that unbiased estimators could also be obtained for functions $f$ that are nonlinear in the second argument. For instance we can get an unbiased estimator of $\{E(\lambda)\}^{k}$, by sampling $k$ independent estimators of $\hat{E}(\lambda)$ and taking their product. More generally we can obtain unbiased estimators of $f(\lambda,E\left\lparen\lambda\right\rparen)$ given $\lambda$ for functions $f$ that are polynomials in the second argument.

Finally it is possible to adapt the proposed approach to estimate the Bayesian cross-validation objective associated with some other scoring rules, such as the one proposed in Hyvärinen, 2005, and considered in the setting of model comparison in e.g. Dawid and Musio, 2015; Shao et al., 2018.

The authors are grateful to Jeremy Heng and Stephane Shao for helpful discussions.

References