Simple, Scalable and Accurate Posterior Interval Estimation

Introduction

We propose a posterior interval estimation algorithm for uncertainty quantification in massive data settings in which usual Bayesian sampling algorithms are too slow. Bayesian models quantify uncertainty via the joint posterior distribution of the model parameters and predictive distributions of new observations. As joint posteriors and predictives are difficult to visualize and use in practice, the focus is almost always on posterior summaries of one-dimensional functionals. For example, it is typical to report 95% posterior credible intervals for a variety of one-dimensional functionals of interest. In practice, by far the most common approach to estimate credible intervals relies on running a Markov chain Monte Carlo algorithm to obtain samples from the joint posterior, based on which estimating intervals for different one-dimensional functionals is trivial. Traditional Markov chain Monte Carlo algorithms are too slow to be practically useful in massive data applications. However, given their rich history and broad use, it would be appealing to be able to incorporate a simple fix-up, which would allow trivial modifications of existing code, solve the computational bottleneck, and enable provably accurate estimation of posterior quantiles for any one-dimensional functional of interest.

Current classes of analytic approximations, such as Gaussian/Laplace, variational Bayes , and expectation propagation , clearly do not provide a generally useful alternative to sampling methods in terms of accurate estimation of posterior credible intervals. Hence, in comparing with the literature, we focus on scalable sampling algorithms. There has been a recent interest in scaling up Bayesian sampling in general and Markov chain Monte Carlo algorithms in particular, with many different threads considered. Three of the most successful include (i) approximating expensive Markov chain Monte Carlo transition kernels with easier to sample surrogates; (ii) running Markov chain Monte Carlo on a single machine but with different subsets of the data used as sampling proceeds ; and (iii) running Markov chain Monte Carlo in parallel for different data subsets and then combining . Motivated by our goal of defining a very simple and theoretically supported algorithm, we focus on embarassingly parallel Markov chain Monte Carlo following strategy (iii).

The key question in embarassingly parallel Markov chain Monte Carlo is how to combine samples from the different subset posteriors. If each subset posterior were approximately Gaussian, then weighted averaging is well justified, motivating the consensus Monte Carlo algorithm . Outside of this restrictive setting, one can instead rely on the product equation representation to combine using kernel smoothing or multi-scale histograms . Such approaches have theory support in terms of accuracy as the number of samples increases, but rely heavily on the accuracy of density estimators for the subset posteriors, suffering badly when subset posteriors have even slightly non-overlapping supports. Moreover, the product equation representation obtained by splitting the prior is not invariant to model reparameterization. An alternative approach is to use data subsamples to define noisy approximations to the full data posterior, and then take an appropriate notion of geometric center, such as geometric median or mean of these approximations. These later approaches are invariant to model reparameterization, but they require a somewhat conceptually and computationally complex combining algorithm.

In this article, we propose a new scalable algorithm for posterior interval estimation. Our algorithm first runs Markov chain Monte Carlo or any alternative posterior sampling algorithm in parallel for each subset posterior, with the subset posteriors proportional to the prior multiplied by the subset likelihood raised to the full data sample size divided by the subset sample size. To obtain an accurate estimate of a posterior quantile for any one-dimensional functional of interest, we simply calculate the quantile estimates in parallel for each subset posterior and then average these estimates. Hence, our combining step is completely trivial conceptually and computationally. We also provide theory justifying the performance of the quantile estimates. We emphasize that we are not proposing a new Markov chain Monte Carlo algorithm, but we are instead developing a simple approach to scale up existing algorithms to datasets with large numbers of observations.

Our approach is related to the frequentist Bag of Little Bootstraps and provides a Bayesian interpretation. Bag of Little Bootstraps divides massive data into small subsets and obtains bootstrap confidence intervals for a one-dimensional parameter on every subset from weighted bootstrap samples. Then the confidence interval of the one-dimensional parameter based on the whole data is constructed by averaging lower and upper bounds of the bootstrap confidence intervals across all subsets. Similarly, our algorithm averages quantiles from all subset posteriors. Our theory leads to new insights into Bag of Little Bootstraps, showing that its confidence intervals correspond to the confidence intervals of the Wasserstein barycenter of bootstrap distributions across all subsets.

Preliminaries

Our algorithm is related to the concept of Wasserstein barycenter of subset posteriors , which depends on the notions of Wasserstein distance and Wasserstein barycenter. Suppose $\Theta\in\mathcal{R}^{d}$ and $\|\theta_{1}-\theta_{2}\|$ is the Euclidean distance between any $\theta_{1},\theta_{2}\in\Theta$. For any two measures $\nu_{1},\nu_{2}$ on $\Theta$, their Wasserstein-2 distance is defined as

which can be viewed as the geometric center of the $N$ measures $\nu_{1},\ldots,\nu_{N}$.

2 Wasserstein Posterior and Posterior Interval Estimation

and we denote their corresponding distribution functions as $\Pi_{n}(\theta\mid X)$ and $\Pi_{m}(\theta\mid X_{j})$, respectively. In the definition of subset posterior density $\pi_{m}(\theta\mid X_{j})$, we have raised the subset likelihood function to the $K$th power. As a stochastic approximation to the overall posterior $\pi_{n}\left(\theta\mid X\right)$, this modification rescales the variance of each subset posterior given $X_{j}$ to be roughly of the same order as the variance of the overall posterior $\Pi_{n}(\theta\mid X)$, as in and . Based on (2.2), runs Markov chain Monte Carlo algorithms on the $K$ subsets in parallel, producing draws from each $\Pi_{m}(\theta\mid X_{j})$, $j=1,\ldots,K$. Empirical estimates of $\Pi_{m}(\theta\mid X_{j})$ for all $K$ subsets are obtained from the Markov chain Monte Carlo draws, their Wasserstein barycenter is estimated via a linear program, and used as an approximation of the overall posterior $\Pi_{n}(\theta\mid X)$.

Suppose we are interested in a scalar parameter $\xi=h(\theta)\in\Xi$ with $h:\Theta\mapsto\Xi\subseteq\mathcal{R}$. We denote the overall posterior for $\xi$ by $\Pi_{n}(\xi\mid X)$ and the $j$th subset posterior for $\xi$ by $\Pi_{m}(\xi\mid X_{j})$. For theory development, we mainly focus on the linear functional $\xi=h(\theta)=a^{\top}\theta+b$ for some fixed $a\in\mathcal{R}^{d}$ and $b\in\mathcal{R}$, which includes the individual components in $\theta$ as special cases. We can define the $W_{2}$ distance and the set of measures $\mathcal{P}_{2}(\Xi)$ on the univariate space $\Xi$. If $\Pi_{m}\left(\xi\mid X_{j}\right)\in\mathcal{P}_{2}(\Xi)$ for all $j=1,\ldots,K$, then the one-dimensional Wasserstein posterior $\overline{\Pi}_{n}(\xi\mid X)$ is defined as the Wasserstein barycenter of $\Pi_{m}\left(\xi\mid X_{j}\right)$ as in (1):

In the one-dimensional case, the Wasserstein posterior has an explicit relation with the $K$ subset posteriors. Let $F^{-1}(u)=\inf\{x:F(x)\geq u\}$ be the quantile function of a generic univariate distribution function $F(x)$. Let $F_{1}$ and $F_{2}$ be two univariate distributions in $\mathcal{P}_{2}(\Xi)$, with quantile functions $F_{1}^{-1}(u)$ and $F_{2}^{-1}(u)$, for any $u\in(0,1)$, respectively. Then the $W_{2}$ distance between $F_{1}$ and $F_{2}$ has an explicit expression by Lemma 8.2 of :

Therefore, $\overline{\Pi}_{n}(\xi\mid X)$ in (3) is explicitly related to the subset posteriors $\Pi_{m}\left(\xi\mid X_{j}\right)$ by

where $\Pi_{m}^{-1}\left(u\mid X_{j}\right)$ and $\overline{\Pi}_{n}^{-1}(u\mid X)$ are the quantile functions of $\Pi_{m}(\xi\mid X_{j})$ and $\overline{\Pi}_{n}(\xi\mid X)$, respectively. This expression for the one-dimensional $W_{2}$ barycenter has been derived in from an optimal transport perspective. The relation indicates that for a scalar functional $\xi$, the average of subset posterior quantiles produces another quantile function that corresponds exactly to the one-dimensional Wasserstein posterior. Therefore, in our algorithm, to evaluate the Wasserstein posterior of $\xi$, we simply take the empirical quantiles based on posterior draws from each $\Pi_{m}(\xi\mid X_{j})$ and then average them over $j=1,\ldots,K$. Our algorithm is summarized in Algorithm 1.

Main Results

In this section, we develop theory supporting our approach. Under mild regularity conditions, we show that the one-dimensional Wasserstein posterior $\overline{\Pi}_{n}(\xi\mid X)$ is an accurate approximation to the overall posterior $\Pi_{n}(\xi\mid X)$. Essentially, as the subset sample size $m$ increases, the $W_{2}$ distance between them diminishes at a faster than parametric rate $o_{p}(m^{-1/2})$. Their biases, variances and quantiles are only different in high orders of $m$. This rate can be improved to $o_{p}(n^{-1/2})$ when the maximum likelihood estimator of $\xi$ is unbiased. Our results are improved relative to previous papers relying on combining subset posteriors, such as and , with more detailed description of the limiting behavior of the estimated posterior and weaker restrictions on the growth rate of the number of subsets $K$.

Our theory relies on the parametric Bernstein-von Mises theorem. The consensus Monte Carlo algorithm in also leverages approximate normality in their asymptotic justification and can be viewed as a different way of averaging subset posteriors. They used weighted averages of subset posterior samples as an approximate sample from the true posterior, where the weights were taken as the inverse covariance matrices based on each subset posterior samples. Their weighting strategy relies more heavily on the normality assumption than our strategy of averaging quantiles. In contrast to the heuristic arguments in , we provide formal justification for using normal approximations on a large number of subsets, and quantify the asymptotic orders of the induced approximation errors.

We make the following assumptions on the data generating process, the prior and the posterior.

$\theta_{0}$ is an interior point of $\Theta\in\mathcal{R}^{d}$, where $d$ is a fixed positive integer and does not depend on $n$. $P_{\theta}=P_{\theta_{0}}$ almost everywhere if and only if $\theta=\theta_{0}$. $X$ contains independent and identically distributed observations generated from $P_{\theta_{0}}$.

The support of $p(x\mid\theta)$ is the same for all $\theta\in\Theta$.

$\log p(x\mid\theta)$ is three times differentiable with respect to $\theta$ in a neighborhood $B_{\delta_{0}}(\theta_{0})\equiv\{\theta\in\Theta:\|\theta-\theta_{0}\|\leq\delta_{0}\}$ of $\theta_{0}$, for some constant $\delta_{0}>0$. $E_{P_{\theta_{0}}}\left\{p^{\prime}(X\mid\theta_{0})/p(X\mid\theta_{0})\right\}=0$. Furthermore, there exists an envelope function $M(x)$ such that $\sup_{\theta\in B_{\delta_{0}}(\theta_{0})}\left|\partial\log p(x\mid\theta)/\partial\theta_{l_{1}}\right|\leq M(x)$, $\sup_{\theta\in B_{\delta_{0}}(\theta_{0})}\left|\partial^{2}\log p(x\mid\theta)/\partial\theta_{l_{1}}\partial\theta_{l_{2}}\right|\leq M(x)$, $\sup_{\theta\in B_{\delta_{0}}(\theta_{0})}\left|\partial^{3}\log p(x\mid\theta)/\partial\theta_{l_{1}}\partial\theta_{l_{2}}\partial\theta_{l_{3}}\right|\leq M(x)$ for all $l_{1},l_{2},l_{3}=1,\ldots,d$, for all values of $x$, and $E_{P_{\theta_{0}}}M(X)^{4}<\infty$.

Let $\psi(X_{1})=E_{\Pi_{m}(\theta\mid X_{1})}Km\|\theta-\hat{\theta}_{1}\|^{2}$, where $E_{\Pi_{m}(\theta\mid X_{1})}$ is the expectation with respect to $\theta$ under the posterior $\Pi_{m}(\theta\mid X_{1})$. Then there exists an integer $m_{0}\geq 1$, such that $\left\{\psi(X_{1}):m\geq m_{0},K\geq 1\right\}$ is uniformly integrable under $P_{\theta_{0}}$. In other words, $\lim_{C\to+\infty}\sup_{m\geq m_{0},K\geq 1}E_{P_{\theta_{0}}}\psi(X_{1})I\{\psi(X_{1})\geq C\}=0$, where $I(\cdot)$ is the indicator function.

Assumptions 1-5 are standard and mild regularity conditions on the model $P(x\mid\theta)$, which are similar to the assumptions of Theorem 8.2 in Chapter 6 of and Theorem 4.2 in for showing the asymptotic normality of posteriors. Assumption 6 requires the prior to have a finite second moment, such that with high probability all the posterior distributions are in the $\mathcal{P}_{2}(\Theta)$ space and the $W_{2}$ distance is well defined. In models with heavy tailed priors, such as our example in Section D.1, one can replace Assumption 6 by assuming that the posterior distribution conditional on a fixed number of initial observations has finite second moment; see Example 8.5 in Chapter 6 of and our Proposition 3 in the Appendix. The uniform integrability of subset posteriors in Assumption 7 is an extra mild technical assumption that helps us to generalize the usual Bernstein-von Mises result on the subsets from the convergence in probability to the convergence in $L_{1}$ distance. We verify Assumption 7 for normal linear models and some exponential family distributions in the Appendix. A stronger condition that can replace Assumption 7 is $\sup_{m\geq m_{0},K\geq 1}E_{X_{1}}E_{\Pi_{m}(\theta\mid X_{1})}Km\|\theta-\hat{\theta}_{1}\|^{2}<+\infty$. The following theorems hold for the one-dimensional Wasserstein posterior defined in (3).

Suppose Assumptions 1–7 hold and $\xi=a^{\top}\theta+b$ for some fixed $a\in\mathcal{R}^{d}$ and $b\in\mathcal{R}$. Let $I_{\xi}(\theta_{0})=\left\{a^{\top}I^{-1}(\theta_{0})a\right\}^{-1}$. Let $\overline{\xi}=a^{\top}\overline{\theta}+b$, $\hat{\xi}=a^{\top}\hat{\theta}+b$. Let $\Phi(\cdot;\mu,\Sigma)$ be the normal distribution with mean $\mu$ and variance $\Sigma$. (i) As $m\to\infty$,

where the convergence is in $P_{\theta_{0}}$-probability. (ii) If $\hat{\theta}_{1}$ is an unbiased estimator for $\theta$, so $E_{P_{\theta_{0}}}\hat{\theta}_{1}=\theta_{0}$, then as $m\to\infty$,

Theorem 1 shows that both the one-dimensional Wasserstein posterior of $\xi$ from combining $K$ subset posteriors and the overall posterior of $\xi$ based on the full dataset are asymptotically close in the $W_{2}$ distance to their respective limiting normal distributions, with slightly different means and the same variance. Such convergence in the $W_{2}$ distance implies weak convergence and convergence of the second moment. Furthermore, the $W_{2}$ distance between the Wasserstein and full posteriors converges to zero in probability with rates $m^{1/2}$ and $n^{1/2}$, depending on the behavior of the maximum likelihood estimator $\hat{\theta}_{1}$.

Previous asymptotic justifications for embarrassingly parallel Markov chain Monte Carlo approaches focus on consistency or convergence rates , while the above theorem is stronger in providing a limiting distribution. In addition, our conditions are much weaker in only requiring the subset sample size $m$ to increase, while imposing no restrictions on the growth rates of $m$ and $K$. Hence, the number of subsets $K$ can grow polynomially in $n$, mimicking the case in which many computers are available but computational resources per computer are limited. For example, the theorem allows $K=O(n^{c})$, $m=O(n^{1-c})$ for any $c\in(0,1)$. Under this setup, the one-dimensional Wasserstein posterior, the overall posterior and their normal limits will all converge to $\theta_{0}$ at the same rate of $O_{p}(n^{-1/2})$, and their mutual difference is of order $o_{p}(m^{-1/2})$.

When the maximum likelihood estimator $\hat{\theta}_{1}$ is unbiased, Part (ii) of the theorem provides a sharper convergence rate of $O_{p}(n^{-1/2})$ compared to the $O_{p}(m^{-1/2})$ rate in Part (i), still with no explicit restrictions on the growth rates of $m$ and $K$. When $K$ increases very fast, for example $K\approx n^{1/2}$ and $m\approx n^{1/2}$, the $O_{p}(n^{-1/2})$ rate in Part (ii) is much faster than the $O_{p}(n^{-1/4})$ rate from Part (i). Moreover, $O_{p}(m^{-1/2})$ is suboptimal since it is the parametric rate based on only the subset data with size $m$, while $O_{p}(n^{-1/2})$ is the optimal parametric rate based on the full data with size $n$. The reason for the improvement in Part (ii) lies in the high order difference between the two means $\overline{\xi}$ and $\hat{\xi}$ of the limiting normal distributions of the one-dimensional Wasserstein posterior and the overall posterior. When the unbiasedness assumption does not hold and $K$ increases with $n$, the difference between the averaged maximum likelihood estimator $\overline{\xi}$ and the overall maximum likelihood estimator $\hat{\xi}$ is typically of order $o_{p}(m^{-1/2})$, which does not scale in the number of subsets $K$. However, when all subset maximum likelihood estimators are unbiased, this difference is reduced by a factor of $K^{1/2}$ due to the averaging effect over $K$ subset posteriors and decreases faster as $o_{p}(n^{-1/2})$. Hence, in models having unbiased maximum likelihood estimators, the one-dimensional Wasserstein posterior achieves high order accuracy in approximating the overall posterior with a difference $o_{p}(n^{-1/2})$.

Independently, has considered a nonparametric generalized linear model and shown a related Bernstein-von Mises theorem. Besides the difference between the form of models, we emphasize that our result in Theorem 1 does not rely on the strong requirement of a uniform normal approximation for all subset posteriors, as used in Shang and Cheng’s paper. Instead, to show Theorem 1, it is only necessary for the normal approximation to work well on average among all subset posteriors. As a result, we have no explicit constraint on the growth rate on the number of subsets $K$, while their paper needs to control $K$ explicitly depending on the posterior convergence rate.

Suppose Assumptions 1–7 hold. Let $\xi_{0}=a^{\top}\theta_{0}+b$ and $\hat{\xi}$ be the same as defined in Theorem 1. For a generic distribution $F$ on $\Xi$, let $\bias(F)=E_{F}(\xi)-\xi_{0}$ and $\var(F)$ be the variance of $F$. Let $u_{1}$ and $u_{2}$ be two arbitrary fixed numbers such that $0<u_{1}<u_{2}<1$. Then the following relations hold:

where $O_{p}$ and $o_{p}$ are in $P_{\theta_{0}}$-probability. Furthermore, if $\hat{\theta}_{1}$ is an unbiased estimator of $\theta_{0}$, then

Theorem 2 provides the order for the differences for the bias, the variance and the quantiles between the one-dimensional Wasserstein posterior and the overall posterior. Essentially the one-dimensional Wasserstein posterior has an asymptotic bias $\overline{\xi}-\hat{\xi}$ from the overall posterior, which is generally of order $o_{p}(m^{-1/2})$ and has higher order $o_{p}(n^{-1/2})$ when the subset maximum likelihood estimators are unbiased. The variances of the one-dimensional Wasserstein posterior and the overall posterior agree in the leading order. Similar to the biases, the difference between their quantiles has order $o_{p}(m^{-1/2})$ in the general case, and improves to a higher order $o_{p}(n^{-1/2})$ when the subset maximum likelihood estimators are unbiased. In our algorithm, when we take $K$ different subset posterior credible intervals and average them, the averages of the lower and upper quantiles are asymptotically close to the quantiles from the overall posterior in the leading order. Therefore, Theorem 2 also validates our algorithm in the sense of posterior uncertainty quantification. We can also account for Monte Carlo errors in approximating subset posteriors using samples under mild mixing conditions on the subset Markov chains; see Theorem 3 in the Appendix.

Experiments

We applied the proposed algorithm in a variety of cases, using consensus Monte Carlo , Wasserstein posterior , semiparametric density product , and variational Bayes as our competitors. Posterior summaries from Markov chain Monte Carlo applied to the full data served as the benchmark for all the comparisons. As our theory guarantees good performance for very large samples, we focused on simulations with moderate sample sizes. All Markov chain Monte Carlo algorithms were run for 10,000 iterations. After discarding the first 5000 samples as burn-in, we retained every fifth sample in all the chains; convergence diagnostics suggested that every chain had converged to its stationary distribution. We used the combination step implemented in R package parallelMCMCcombine for consensus Monte Carlo and semiparametric density product methods. We implemented the combination step of our algorithm in R and of Srivastava et al.’s algorithm in Matlab. All experiments were run on an Oracle Grid Engine cluster with 2.6GHz 16 core compute nodes. Memory resources were capped at 8GB for all the methods, except for Markov chain Monte Carlo based on the full data, which had a maximum memory limit of 16GB.

The accuracy of a density $q(\theta\mid X)$ approximating $\pi_{n}(\theta\mid X)$ was evaluated using the metric

This accuracy metric lies in $$, with larger values indicating better performance of$q$in approximating$\pi_{n}$. In our experiments, we first estimated$q(\theta\mid X)$and$\pi_{n}(\theta\mid X)$ based on the posterior samples using the bkde or bkde2D functions in R package KernSmooth, with automatic bandwidth selection via dpik . The density estimates were used to compute a numerical approximation of the integral in (4).

We first evaluated the performance of our proposed algorithm under varying sample size, dimension, and number of subsets in Bayesian linear models. Let the response, design matrix, regression coefficients, and random error be denoted as $y$, $X$, $\beta$, and $\epsilon$, where $y,\epsilon\in\mathcal{R}^{n}$, $\beta\in\mathcal{R}^{p\times 1}$, and $X\in\mathcal{R}^{n\times p}$. The model assumes that

where gdP denotes the generalized double Pareto shrinkage prior of and Half-$t$ is chosen to be weakly-informative . See Section D.1 in the Appendix for detailed specifications. The priors on $\beta$ and $\sigma$ in (5) are both heavy-tailed with infinite second moments, and therefore do not satisfy Assumption 6. However, one can verify that conditional on the initial $m_{0}$ observations with $m_{0}\geq p+4$, every subset posterior has finite second moments for both $\beta$ and $\sigma$. The result is summarized in Proposition 3 in the Appendix.

We applied our approach for inference on $\beta$ in (5) compared with an asymptotic normal approximation. We calculated the accuracy of approximations using a full data Gibbs sampler as the benchmark (Table 1). The first $10\%$ of entries of $\beta$ were set to $\pm 1$ with the remaining 0. The entries of $X$ were randomly set to $\pm 1$ and $\sigma^{2}$ was fixed at 1. We ran 10 replications for $n\in\{10^{4},10^{5}\}$ and $p\in\{10,100,200,300,400,500\}$. We varied $K\in\{10,20\}$ and applied Algorithm 1 after running a modification of the Gibbs sampler in (2.2) for each subset. We considered two versions of normal approximations for the full posterior. The first version used $\mathcal{N}(\widehat{m},\widehat{V})$ to approximate the posterior of $\beta$, where $\widehat{m}$ and $\widehat{V}$ are the maximum likelihood estimates of $\beta$ and its estimated asymptotic covariance matrix in (5). For the second version, we first obtained the asymptotic normal approximation of the $j$th subset posterior as $\mathcal{N}$($\widehat{m}_{j}$, $\widehat{V}_{j}$), where $\widehat{m}_{j}$ and $\widehat{V}_{j}$ ($j=1,\ldots,K$) are the maximum likelihood estimates of $\beta$ and its estimated asymptotic covariance matrix for the $j$th subset. Then we found the $W_{2}$ barycenter of the $K$ subset normal approximations, which is again a normal distribution $\mathcal{N}(m^{*},V^{*})$ . This provides an empirical verification of Theorem 1. See Section D.1 in the Appendix for details of the Gibbs sampler and the form of $m^{*}$ and $V^{*}$.

The performance of all the approaches was fairly similar across all simulations and agreed with our asymptotic theory (Table 1). The results in Table 1 show that the proposed algorithm closely matched the Gibbs sampling results for the full data in terms of uncertainty quantification. It also performed better than the asymptotic normal approximations in some cases. When the subset sample size was too small compared to the dimension, such as when $n=10^{4},p=400,K=20$ which has a subset size of only $m=500$, we observe poor performance for both the asymptotic approximations and the proposed approach.

2 Linear mixed effects model

Linear mixed effects models are widely used to characterize dependence in longitudinal and nested data structures. Let $n_{i}$ be the number of observations associated with the $i$th individual, for $i=1,\ldots,s$. Let $y_{i}\in\mathcal{R}^{n_{i}}$ be the responses of the $i$th individual, $X_{i}\in\mathcal{R}^{n_{i}\times p}$ and $Z_{i}\in\mathcal{R}^{n_{i}\times q}$ be matrices including predictors having coefficients that are fixed across individuals and varying across individuals, respectively. Let $\beta\in\mathcal{R}^{p}$ and $u_{i}\in\mathcal{R}^{q}$, respectively, represent the fixed effects and $i$th random effect. The linear mixed effects model lets

Many software packages are available for Markov chain Monte Carlo-based Bayesian inference in (6), but current implementations become intractable for data with large $s$ and $n=\sum_{i=1}^{s}n_{i}$.

We applied our algorithm for inference on $\beta$ and $\Sigma$ in (6) and compared its performance with maximum likelihood, consensus Monte Carlo, semiparametric density product, Wasserstein posterior, and variational Bayes. We set $s=5000$, $n_{i}=20$ for $i=1,\ldots,s$, $n=10^{5}$, $p=4$, $q=3$, $\beta=(-1,1,-1,1)^{\top}$, and $\sigma=1$. The random effects covariance $\Sigma$ had $\Sigma_{ii}=i,i=1,2,3$, $\Sigma_{12}=-$ 0.56, $\Sigma_{31}=$ 0.52, and $\Sigma_{23}=$ 0.0025. This matrix included negative, positive, and small to moderate strength correlations . The simulation was replicated 10 times. The approximate posterior distributions were obtained using consensus Monte Carlo, semiparametric density product, Wasserstein posterior, and our algorithm in three steps. First, full data were randomly partitioned into 20 subsets such that data for each individual were in the same subset. Second, the Markov chain Monte Carlo sampler for $\beta$ and $\Sigma$ in (6) was modified following (2.2) and Equation (2) in and implemented in Stan (Version 2.5.0). Finally, the posterior samples from all the subsets were combined. We used the streamlined algorithm for variational Bayes . Maximum likelihood produced a point estimate and asymptotic covariance for $\beta$, and only a point estimate for $\Sigma$.

We compared the performance of the seven methods for inference on the fixed effects $\beta$, the variances of random effects $\Sigma_{ii}$ ($i=1,2,3$), and the correlations of random effects $\rho_{ij}=\Sigma_{ij}/(\Sigma_{ii}\Sigma_{jj})^{1/2}$ ($1\leq i<j\leq 3$). The correlations are nonlinear functionals of the model parameters $\Sigma$. Maximum likelihood estimator, consensus Monte Carlo, semiparametric density product, Wasserstein posterior, and our algorithm had excellent performance in estimation of $\beta$ (Figure 1), as well as the variances and correlations (Tables 2 and 3 and Figure 2). Uncertainty quantification using consensus Monte Carlo, semiparametric density product, Wasserstein posterior, and our algorithm closely agreed with Markov chain Monte Carlo based on the full data. As shown in Figure 1 of the Appendix, variational Bayes was computationally most efficient, but it showed poor accuracy in approximating the posterior of $\beta$ and the variances, with underestimation of posterior uncertainty.

3 United States natality data

We applied our algorithm to United States natality data on birth weight of infants and variables related to their mothers’ health . Linear mixed effects models were used for the covariance in birth weights among siblings. Following the example in , we selected the data for mothers who smoked, had two infants, and had some college education but not a college degree. Detailed information about the variables are in the Appendix. The data set contained $s=3809$ mothers and $n=7618$ births. There were 13 variables related to mother’s health. All these covariates and an intercept were used as fixed effects in (6), so $p=14$. The random effects included mother’s age, gestation period, and number of living infants, so $q=3$. We performed 10 fold cross-validation and randomly split the data into 10 data sets such that data for siblings belonged to the same training data. We estimated fixed effects and covariance matrix for random effects as in Section 4.2 using $K=20$.

The seven methods in the previous section generally agreed in the inference on fixed effects (Figure 3), with variational Bayes deviating the furthest. Our algorithm and the algorithm of differed significantly from variational Bayes, consensus Monte Carlo, and semiparametric density product in the inference on variances and correlations of random effects (Tables 4 and 5 and Figure 4). Our algorithm and the algorithm of showed better agreement with Markov chain Monte Carlo based on the full data in estimating the correlations. The 90% credible intervals from our algorithm included the maximum likelihood estimates of correlations. Variational Bayes posterior concentrated very close to 0 for every element of the covariance matrix and significantly underestimated posterior uncertainty. Consensus Monte Carlo and semiparametric density product methods performed poorly in the inference on random effects but were better than variational Bayes. Markov chain Monte Carlo based on the full data was extremely slow compared to the other methods (see Figure 1 in the Appendix). Taking into account both the approximation accuracy and the computational efficiency, we concluded that our proposed algorithm performs better than the competing algorithms in estimating the covariance matrix of random effects.

4 Extension to multi-dimensional parameters

Although Algorithm 1 only applies to one-dimensional functionals, we provide a simple extension to the multi-dimensional case with a numerical illustration. Suppose our goal is to find the joint posterior of the $d$-dimensional parameter $\theta$. First, we center and scale the posterior samples of $\theta$ in every subset. Let $\widehat{m}_{j}$ and $\widehat{V}_{j}$ be the empirical mean and covariance matrix for the $j$th subset posterior samples $\{\theta_{1j},\ldots,\theta_{Tj}\}$. Let $\widehat{m}=K^{-1}\sum_{j=1}^{K}\widehat{m}_{j}$, $\widehat{V}^{-1}=K^{-1}\sum_{j=1}^{K}\widehat{V}_{j}^{-1}$. We transform every subset draw $\theta_{ij}$ to $\theta^{\prime}_{ij}=\widehat{V}^{-1/2}(\theta_{ij}-\widehat{m})$. If every subset posterior of $\theta$ is asymptotically normal, then the centered and rescaled version $\theta^{\prime}$ will be asymptotically standard normal with approximately independent components, since $T$ is large in practice. For every component of $\theta^{\prime}$, we apply Algorithm 1 to combine its $K$ subset posterior samples and obtain approximations of posterior quantiles for a fine grid of $$. This leads to accurate approximations of the marginal posteriors of$\theta^{\prime}$; we repeatedly draw samples from these marginals, and then transform back to the original parameter using$\theta=\widehat{V}^{1/2}\theta^{\prime}+\widehat{m}$. This yields approximate samples from the full posterior of$\theta$, and credible regions can be estimated based on these samples.

We implemented this generalized algorithm for combining subset posterior samples of all pairs of variances and covariances in the simulation from Section D.2, and compared the results with consensus Monte Carlo, semiparametric density product, Wasserstein posterior, and variational Bayes. The accuracies of our algorithm and the algorithm of were higher than the accuracies of the other three methods for all pairs of variances and covariances (Table 6). Variational Bayes performed poorly in the estimation of posterior distributions for all the pairs of variances. We obtained kernel density estimates of the three pairs of covariances in (6) using the combined posterior samples and the bkde2D function in the KernSmooth R package with a bandwidth of 0.01 (Figure 5). The kernel density estimates centered very close to the true values of the covariance pairs. Compared to the algorithm of , our algorithm was more efficient, easier to implement, and robust to the grid-size of quantiles, while having similar accuracy and stability across all simulation replications.

Conclusion

We have proposed a simple posterior interval estimation algorithm to rapidly and accurately estimate quantiles of the posterior distributions for different one-dimensional functionals of interest. The algorithm is simple and efficient relative to existing competitors, just averaging quantile estimates for each subset posterior based on applying existing sampling algorithms in an embarrassingly parallel manner. There is a fascinating mathematical relationship with the Wasserstein barycenter of subset posteriors: our algorithm calculates quantiles of the Wasserstein posterior without the optimization step in . The credible intervals from our algorithm asymptotically approximate those from the full posterior in the leading parametric order. The quality of approximation is the same even if the subset sample size increases slowly and the number of subsets increases polynomially fast. Our experiments have demonstrated excellent performance for linear mixed effects models and linear models with varying dimension.

Although our current theory focuses on parametric models and one-dimensional linear functionals, the proposed algorithm can be practically implemented for general one-dimensional functionals for semiparametric and nonparametric models. For example, in simulations not shown in the paper, we found that our algorithm shows excellent performance for Dirichlet process mixture models for multivariate categorical data , and Gaussian process nonparametric regression. Furthermore, we have provided an extension to the multi-dimensional case. It would be appealing to develop theory justification in these more complex settings, and to develop guarantees on approximation accuracy for fixed subset sizes and growing numbers of subsets. Also of interest in future work is to consider algorithms that do not require non-overlapping subsets, potentially relying on subsampling. Although such modifications can be implemented trivially, our proof techniques for the combining step in Theorem 1 do not apply directly. Other important extensions include optimal design of subsampling algorithms and extensions beyond product likelihoods.

In Section A we provide the detailed technical proofs of Theorem 1 and Theorem 2 in the main paper. In Section B, we present a theorem that quantifies the Monte Carlo errors in subset posterior sampling. In Section C, we verify Assumption 7 in the main paper for the normal linear model and some exponential family distributions. Section D includes further details about the data analysis in the main paper. In particular, for the heavy tailed priors used in Example 1 in the main paper, we verify a relaxed version of Assumption 6 in Section D.1.

Appendix A Proofs of Theorem 1 and Theorem 2

(Villani Theorem 6.15) For two measures $P_{1},P_{2}\in\mathcal{P}_{2}(\Theta)$, or similarly $\mathcal{P}_{2}(\Xi)$,

where the total variation of moments distance is defined as

In comparison, the usual parametric Bernstein-von Mises theorem on the subset $X_{j}$ without raising the likelihood to the $K$th power gives

in $P_{\theta_{0}}$-probability, where $z=m^{1/2}(\theta-\hat{\theta}_{j})$. See, for example, Theorem 8.2 in and Theorem 4.2 in .

Proof of Lemma 2: The relation (A.2) in Lemma 2 is the usual Bernstein-von Mises theorem for the overall posterior $\Pi_{n}(\theta\mid X)$. The proof of (A.2) follows a related line to the proof of Theorem 4.2 in , and can be treated as a special case of (A.1) with $m=n$ and $K=1$. In the following we focus on the proof of (A.1) in Lemma 2.

Given the independent and identically distributed assumption, we only need to show the result for a fixed index $j$. To emphasize the different roles played by the subset sample size $m$ and the number of subsets $K$, in the following proofs we will write the total sample size $n$ as $Km$. We complete the proof in 3 steps. For a generic matrix $A$ or a 3-dimensional array $A$, we use $\|A\|$ to denote its Frobenius norm.

Step 2: Show the following convergence as $m\to\infty$ in $P_{\theta_{0}}$-probability:

We prove this result for a fixed subset $X_{j}$, since the data are independent and identically distributed, and the conclusion is identical for any $j=1,\ldots,K$. Define the following quantities

Then based on the expression of $\Pi_{m}(\theta\mid X_{j})$, with the likelihood raised to the $K$th power,

The induced posterior density on $t=(Km)^{1/2}(\theta-\hat{\theta}_{j})$ can be written as

Let $\mathcal{T}=\{t=(Km)^{1/2}(\theta-\hat{\theta}_{j}):\theta\in\Theta\}$. Define

as $m\to\infty$ in $P_{\theta_{0}}$-probability. Hence, for the difference in (A.3), we obtain that

Divide the domain of the integral into 3 parts: $A_{1}=\{z:\|z\|\geq\delta_{1}(Km)^{1/2}\}$, $A_{2}=\{z:\delta_{2}\leq\|z\|<\delta_{1}(Km)^{1/2}\}$, $A_{3}=\{z:\|z\|<\delta_{2}\}$, where the constants $\delta_{1},\delta_{2}$ will be chosen later. Then

as $m\to\infty$, because the integral on the whole $z\in\mathcal{R}^{d}$ is finite, $\pi(\theta_{0})$ is bounded from above according to Assumption 6, and $K\geq 1$.

Next we bound the first term in (A.1). By Assumption 5 and the weak consistency of $\hat{\theta}_{j}$, there exists a constant $\epsilon_{1}$ that depends on $\delta_{1}$, such that for any $z\in A_{1}$ and all sufficiently large $m$, with $P_{\theta_{0}}$-probability approaching 1,

Furthermore, the weak consistency of $\hat{\theta}_{j}$ implies that for all sufficiently large $m$, with $P_{\theta_{0}}$-probability approaching 1, $\|\hat{\theta}_{j}\|\leq\|\hat{\theta}_{j}-\theta_{0}\|+\|\theta_{0}\|\leq\delta_{0}+\|\theta_{0}\|$. Therefore, as $m\to\infty$ in $P_{\theta_{0}}$-probability,

where we have used the finite second moment of $\pi(\theta)$ from Assumption 6 in the last step. Hence, we have proved that the first integral in (A.5) goes to zero in $P_{\theta_{0}}$-probability.

where the last convergence is almost surely in $P_{\theta_{0}}$ by the strong law of large numbers. Therefore, we can choose $\delta_{1}$ as

where $\lambda_{1}(A)$ denotes the smallest eigenvalue of a generic matrix $A$. Assumption 4 indicates that $\min_{\theta\in B_{\delta_{0}}(\theta_{0})}\lambda_{1}\{I(\theta)\}$ is bounded below by a constant. Thus, in (A.8), the choice of $\delta_{1}$ implies that for every $z\in A_{2}$, for all large $m$ with $P_{\theta_{0}}$-probability approaching 1,

Therefore for $z\in A_{2}$, for all large $m$ with $P_{\theta_{0}}$-probability approaching 1,

For the third integral in (A.5), we fix a constant $\delta_{2}>0$ and can use the similar Taylor series expansion above, and notice that when $\|z\|<\delta_{2}$, as $m\to\infty$,

as $m\to\infty$ in $P_{\theta_{0}}$-probability. Therefore, (A.10) and (A.11) together imply that as $m\to\infty$ in $P_{\theta_{0}}$-probability,

Hence by the definition of $g_{m}(z)$, as $m\to\infty$ in $P_{\theta_{0}}$-probability,

This has proved that the right-hand side of (A.5) converges to zero in $P_{\theta_{0}}$-probability, and also completes the proof of (A.3).

Step 3: Show the convergence in $L_{1}$ as $m\to\infty$. It is clear from the derivation of (A.1) that

In this display, the last term is a finite constant. The middle term is $\psi(X_{j})$ defined in Assumption 7. According to Assumption 7, for any fixed $j$, $\left\{\psi(X_{j}):m\geq m_{0},K\geq 1\right\}$ is uniformly integrable under $P_{\theta_{0}}$. Now since $TV_{2}\left[\Pi_{m,t}(t\mid X_{j}),\Phi\left\{t;0,I^{-1}(\theta_{0})\right\}\right]$ is upper bounded by $\psi(X_{j})+C$ for all $m,K$ and some constant $C>0$, we obtain that $\left\{TV_{2}\left[\Pi_{m,t}(t\mid X_{j}),\Phi\left\{t;0,I^{-1}(\theta_{0})\right\}\right]:m\geq m_{0},K\geq 1\right\}$ is also uniformly integrable. This uniform integrability together with the convergence in $P_{\theta_{0}}$-probability from Step 2 implies the $L_{1}$ convergence of $TV_{2}\left[\Pi_{m,t}(t\mid X_{j}),\Phi\left\{t;0,I^{-1}(\theta_{0})\right\}\right]$ to zero. $\blacksquare$

Similar to the $W_{2}$ distance, for any $l\geq 1$, we can define the Wasserstein-$l$ ($W_{l}$) distance: for any two measures $\nu_{1},\nu_{2}$ on $\Theta$, their $W_{l}$ distance is defined as

where $\Gamma(\nu_{1},\nu_{2})$ is the set of all probability measures on $\Theta\times\Theta$ with marginals $\nu_{1}$ and $\nu_{2}$, respectively. The $W_{l}$ distance on the space $\Xi$ can be similarly defined. The $W_{l}$ distance between two univariate distributions $F_{1}$ and $F_{2}$ is the same as the $L_{l}$ distance between their quantile functions (see Lemma 8.2 of ):

Let $\hat{\xi}_{j}=a^{\top}\hat{\theta}_{j}+b$. Then for any $l\geq 1$,

Proof of Lemma 3: We use $\Phi(\cdot)$ and $\Phi^{-1}(\cdot)$ to denote the cumulative distribution function and the quantile function of standard normal distribution $\mathcal{N}(0,1)$. From , the univariate Wasserstein-2 barycenter satisfies that for any $u\in(0,1)$,

Since $l\geq 1$, we apply Minkowski inequality to the right-hand side of (A.1) and obtain that

which concludes the proof. $\blacksquare$

where $o_{p}$ is in $P_{\theta_{0}}$ probability. Furthermore, if $\hat{\theta}_{1}$ is an unbiased estimator of $\theta_{0}$, then

Proof of Lemma 4: Because of the linearity $\xi=a^{\top}\theta+b$, it suffices to show

with the further assumption that $\hat{\theta}_{1}$ is an unbiased estimator of $\theta_{0}$.

Next we show the first term in (A.18) is of order $o_{p}(m^{-1/2})$ under Assumptions 1–7, and is of order $o_{p}(n^{-1/2})$ if furthermore $E_{P_{\theta_{0}}}\hat{\theta}_{1}=\theta_{0}$.

Hence, $\left\|\sum_{j=1}^{K}W_{j}/(Km^{1/2})\right\|=o_{p}\left(m^{-1/2}\right)$. This together with (A.18) and (A.19) leads to (A.15).

If we further assume unbiasedness $E_{P_{\theta_{0}}}\hat{\theta}_{1}=\theta_{0}$, then from (A.17) we can obtain that

for all $j=1,\ldots,K$. In other words, $W_{j}$’s are centered at zero. Since $X_{j}$’s ($j=1,\ldots,K$) are all independent and $W_{j}$ only depends on $X_{j}$, we have $E_{P_{\theta_{0}}}W_{j_{1}}^{\top}W_{j_{2}}=0$ for any $j_{1}\neq j_{2}$.

We can again apply Markov’s inequality to the first term in (A.18) and obtain that for any constant $c>0$,

Therefore, assuming that $E_{P_{\theta_{0}}}\hat{\theta}_{1}=\theta_{0}$ and $E_{P_{\theta_{0}}}(\|W_{1}\|^{2})\to 0$ as $m\to\infty$, which will be proven below, the display above implies that $\left\|\sum_{j=1}^{K}W_{j}/(Km^{1/2})\right\|=o_{p}\left(n^{-1/2}\right)$. This together with (A.18) and (A.19) leads to (A.16).

where $d$ is the dimension of $\theta$. We have used the property of the Frobenius norm: for a generic $d\times d$ symmetric positive definite matrix $A$, $\|A^{-1}\|\leq d^{1/2}\overline{\lambda}(A^{-1})=d^{1/2}\left\{\underline{\lambda}(A)\right\}^{-1}$, where $\overline{\lambda}(A)$ and $\underline{\lambda}(A)$ denotes the largest and the smallest eigenvalues of the matrix $A$, respectively. Furthermore, the envelop function condition in Assumption 3 implies that

It follows from (A.20) and (A.21) that for all large $m$,

where $c_{1},c_{2}$ are positive constants that only depend on $d,\underline{\lambda},\|I(\theta_{0})\|^{2}$.

Due to Assumption 3, the first term in (A.1) is bounded by

Thus we have shown that both terms on the right-hand side of (A.1) are finite. Therefore, $E_{P_{\theta_{0}}}(V_{1})<\infty$ and by the dominated convergence theorem, $E_{P_{\theta_{0}}}\|W_{1}\|^{2}\to 0$. $\blacksquare$

A.2 Proof of Theorem 1

Proof of Theorem 1(i): Since $\xi=a^{\top}\theta+b$, we can derive the following for subset posteriors in terms of $\xi$ using a change of variable from $\theta$ to $\xi$ in (A.1) of Lemma 2:

where $t=n^{1/2}(\xi-\hat{\xi}_{j})$ is now the local parameter for the $j$th subset. From the relation between norms $W_{2}$ and $TV_{2}$ in Lemma 1, this directly implies

We further use the rescaling property of the $W_{2}$ distance and obtain the equivalent form in terms of the original parameter $\xi$:

From Lemma 3, we have that for any constant $c>0$, as $m\to\infty$,

where (i) follows from Lemma 3 with $l=2$, (ii) uses Markov’s inequality, (iii) comes from the relation between $l_{1}$ norm and $l_{2}$ norm, and (iv) follows from (A.23). This result indicates that

which shows the first relation in Part (i) of Theorem 1. The second relation in Theorem 1 (i)

follows from a similar argument using (A.2) in Lemma 2 for the overall posterior.

From Lemma 4, we have $\left|\overline{\xi}-\hat{\xi}\right|=o_{p}\left(m^{-1/2}\right)$. Therefore,

where the first inequality follows because of the definition of $W_{2}$ distance and the same variance shared by the two normal distributions.

Finally, by (A.24), (A.25), (A.26) and the triangular inequality, we have

which is equivalent to the third relation in Part (i). $\blacksquare$

Proof of Theorem 1(ii): If $\hat{\theta}_{1}$ is an unbiased estimator of $\theta_{0}$, then by Lemma 4 and the definition of $W_{2}$ distance, it follows that

Applying the triangular inequality to (A.24), (A.25) and (A.27), we obtain that as $m\to\infty$,

Thus the conclusion of Part (ii) follows. $\blacksquare$

A.3 Proof of Theorem 2

Proof of Theorem 2(i): have shown that the barycenter $\overline{\Pi}_{n}(\xi\mid X)$ is related to the $K$ subset posteriors $\Pi_{m}(\xi\mid X_{j})$ ($j=1,\ldots,K$) through the quantile function:

To see why this is true, we notice that we have derived the following relation in the proof of Theorem 1:

But according to the definition of $r_{j}(u)$ in (A.13), by Cauchy-Schwarz inequality,

On the other hand, for the bias of the overall posterior $\Pi_{n}(\xi\mid X)$, we follow a similar argument as above and obtain that

where $r(u)=\Pi_{n}^{-1}(u|X)-\hat{\xi}-\left\{nI_{\xi}(\theta_{0})\right\}^{-1/2}\Phi^{-1}(u)$. Moreover we have

by Theorem 1. This completes the proof of Part (i). $\blacksquare$

Proof of Theorem 2(ii): Similar to the expectation, the variance of a generic univariate distribution $F$ can be calculated through its quantile functions: if $Y\sim F$,

From (A.14) (with $l=2$) and the conclusion of Theorem 1, we have

Again by Cauchy-Schwarz inequality, we have

Therefore, we have shown that $\var\left\{\overline{\Pi}_{n}(\xi\mid X)\right\}=\left\{nI_{\xi}(\theta_{0})\right\}^{-1}+o_{p}\left(n^{-1}\right)$.

For the variance of $\Pi_{n}(\xi\mid X)$, we use the same definition of $r(u)$ as in Part (i) and derive that

Based on the conclusion of Theorem 1 and Cauchy-Schwarz inequality, we have

which proves $\var\left\{\Pi_{n}(\xi\mid X)\right\}=\left\{nI_{\xi}(\theta_{0})\right\}^{-1}+o_{p}\left(n^{-1}\right)$. $\blacksquare$

Proof of Theorem 2(iii): The convergence in $W_{2}$ distance implies weak convergence. Therefore, it follows from Theorem 1 that in $P_{\theta_{0}}$ probability, both $\overline{\Pi}_{n,s}(s\mid X)$ and $\Pi_{n,s}(s\mid X)$ converge in distribution to normal distributions as $m\to\infty$. The weak convergence also implies the convergence of quantile functions at any continuous point. Since both $\overline{\Pi}_{n}\left(\xi\mid X\right)$ and $\Pi_{n}\left(\xi\mid X\right)$ are continuous distributions with posterior densities, their quantiles also converge pointwise to the quantiles of their limiting normal distributions. For any fixed $u\in(0,1)$, as $m\to\infty$, Theorem 1 implies that for $s=n^{1/2}(\xi-\overline{\xi})$,

We can make this convergence uniform over all quantiles $u\in[u_{1},u_{2}]\subset(0,1)$. Divide $[u_{1},u_{2}]$ into $L$ equally spaced subintervals $[u_{(j)},u_{(j+1)}]$ for $j=0,\ldots,L-1$ and $u_{(j)}=u_{1}+j(u_{2}-u_{1})/L$. For any $\epsilon>0$, since $\Phi^{-1}\{u;0,I_{\xi}^{-1}(\theta_{0})\}$ is uniformly continuous on $[u_{1},u_{2}]$, we can pick $L$ sufficiently large such that

for all $j=0,\ldots,L-1$. Furthermore, because $\Phi^{-1}(\cdot)$ is continuous everywhere, we can find a sufficiently large $n_{0}$, such that for all $n>n_{0}$, all $j=0,\ldots,L-1$ with the $L$ chosen above,

For any $u\in[u_{1},u_{2}]$, we can find a $j_{0}\in\{0,\ldots,L-1\}$ such that $u\in[u_{(j_{0})},u_{(j_{0}+1)}]$. Therefore using the monotonicity of quantile functions,

which implies that for the quantiles in terms of $\xi$,

By plugging in the order $\left|\overline{\xi}-\hat{\xi}\right|=o_{p}(m^{-1/2})$ from the proof of Theorem 1, we have

If we further assume that $\hat{\theta}_{1}$ is an unbiased estimator of $\theta_{0}$, then Lemma 4 says that $\left|\overline{\xi}-\hat{\xi}\right|=o_{p}\left(n^{-1/2}\right)$. Therefore, using the results from Part (i), we have

which completes the proof. $\blacksquare$

Appendix B Theoretical Results for the Posterior Monte Carlo Errors

In practice, the credible intervals are calculated from the averages of empirical quantiles from subset posterior samples. In Algorithm 1, suppose that for each $j=1,\ldots,K$, $\Pi^{\circ}_{j}(\theta)$ and $\kappa_{j}(\theta,\theta^{\prime})$ for $\theta,\theta^{\prime}\in\Theta$ are the initial distribution and the transition kernel for the Markov chain of the $j$th subset posterior. $\{\theta_{1j},\ldots,\theta_{Tj}\}$ with sample size $T$ are drawn sequentially with $\theta_{1j}\sim\Pi^{\circ}_{j}(\cdot)$ and $\theta_{l+1,j}\sim\kappa_{j}(\theta_{lj},\cdot)$ for $l=1,\ldots,T-1$. $\xi_{lj}=a^{\top}\theta_{lj}+b$ for $l=1,\ldots,T$ and $j=1,\ldots,K$. Let $\widehat{\Pi}_{m}(\xi\mid X_{j})$ be the empirical distribution of $\{\xi_{1j},\ldots,\xi_{Tj}\}$ for $j=1,\ldots,K$. Let $\widehat{\Pi}_{n}(\xi\mid X)$ be the Wasserstein barycenter of $\widehat{\Pi}_{m}(\xi\mid X_{1}),\ldots,\widehat{\Pi}_{m}(\xi\mid X_{K})$, which can be calculated through its quantile function $\widehat{\Pi}^{-1}_{n}(u\mid X)=\sum_{j=1}^{K}\widehat{\Pi}_{m}^{-1}(u\mid X_{j})/K$ for all $u\in(0,1)$. Let $L_{2}\{\Pi_{m}(\cdot\mid X_{j})\}$ for $j=1,\ldots,K$ be the $L_{2}$ space of functions on $\Theta$ such that for any $f\in L_{2}\{\Pi_{m}(\cdot\mid X_{j})\}$, $\|f(\theta)\|^{2}_{L_{2},j}=E_{\Pi_{m}(\cdot\mid X_{j})}f^{2}(\theta)<\infty$ almost surely in $P_{\theta_{0}}$. We need three additional assumptions as follows.

$\max_{1\leq j\leq K}E_{\Pi_{m}(\cdot\mid X_{j})}\|\theta\|^{7}$ is upper bounded by a constant almost surely in $P_{\theta_{0}}$. $~\max_{1\leq j\leq K}E_{\Pi_{m}(\cdot\mid X_{j})}\{\pi^{\circ}_{j}(\theta)/\pi_{m}(\theta\mid X_{j})\}^{3}$ is upper bounded by a constant almost surely in $P_{\theta_{0}}$, where $\pi^{\circ}_{j}(\theta)$ is the density of $\Pi^{\circ}_{j}(\theta)$ for $j=1,\ldots,K$.

Every subset posterior $\Pi_{m}(\theta\mid X_{j})$ ($j=1,\ldots,K$) is $\rho$-mixing: there exists a nonnegative constant sequence $\{\rho_{l}\}_{l\geq 1}$ decreasing to zero and $\sum_{l=1}^{\infty}\rho_{l}<\infty$, such that almost surely in $P_{\theta_{0}}$, for any integer $l\geq 1$, any $f\in L_{2}\{\Pi_{m}(\cdot\mid X_{j})\}$ and all $j=1,\ldots,K$,