Robust Mean Estimation with the Bayesian Median of Means

TL;DR

This paper proposes the Bayesian median of means, a new aggregation method that reduces variance in estimates at the cost of bias, demonstrating improved performance in high-variance scenarios through theoretical analysis and empirical validation.

Contribution

Introduction of the Bayesian median of means, a non-parametric aggregation rule that interpolates between mean and median, with theoretical guarantees and practical algorithms.

Findings

Significantly reduces variance in high-variance settings.

Maintains consistency and asymptotically negligible bias.

Performs well in real-world applications like importance sampling and cross-validation.

Abstract

The sample mean is often used to aggregate different unbiased estimates of a parameter, producing a final estimate that is unbiased but possibly high-variance. This paper introduces the Bayesian median of means, an aggregation rule that roughly interpolates between the sample mean and median, resulting in estimates with much smaller variance at the expense of bias. While the procedure is non-parametric, its squared bias is asymptotically negligible relative to the variance, similar to maximum likelihood estimators. The Bayesian median of means is consistent, and concentration bounds for the estimator's bias and $L_1$ error are derived, as well as a fast non-randomized approximating algorithm. The performances of both the exact and the approximate procedures match that of the sample mean in low-variance settings, and exhibit much better results in high-variance scenarios. The empirical performances are examined in real and simulated data, and in applications such as importance sampling, cross-validation and bagging.

Tables

Table 1. Table 1 : Results for Advertising data

	Linear Regression			Random forest			$k$NN
	CV	BMM	error	CV	BMM	error	CV	BMM	error
mean	3.002	2.939	2.889	0.714	0.7	0.667	2.532	2.514	2.264
std	0.246	0.23	0.937	0.08	0.076	0.227	0.205	0.201	0.592
RMSE	1.184	1.16	-	0.263	0.257	-	0.737	0.725	-

Table 2. Table 2 : Results for Auto data

	Linear Regression			Random forest			$k$NN
	CV	BMM	error	CV	BMM	error	CV	BMM	error
mean	536.23	531.63	508.11	389.7	384.74	366.69	1345.65	1322.42	1235.03
std	27.89	27.19	105.39	36.69	34.28	120.37	100.9	95.69	339.19
RMSE	134.72	133.09	-	135.12	133.92	-	443.72	432.38	-

Table 3. Table 3 : Results for College data

	Linear Regression			Random forest			$k$NN
	CV	BMM	error	CV	BMM	error	CV	BMM	error
mean	172.15	171.802	168.24	180.363	179.883	176.662	212.989	212.687	210.43
std	6.601	6.682	23.478	8.332	8.37	24.348	8.303	8.285	24.506
RMSE	30.044	29.989	-	31.478	31.445	-	31.949	31.861	-

Table 4. Table 4 : Results for Credit data

	Linear Regression			Random forest			$k$NN
	CV	BMM	error	CV	BMM	error	CV	BMM	error
mean	10116.4	10080.4	10451.1	13662	13485.5	12966.4	48018.5	47829.4	47001.5
std	444.9	452.6	1640.4	1201.4	1159.8	3912	2548.6	2599.8	7793.7
RMSE	2098	2111.9	-	4492.9	4481.8	-	10004.7	10027.3	-

Table 5. Table 5 : Results for Heart data

	LDA			Neural network			$k$NN
	CV	BMM	error	CV	BMM	error	CV	BMM	error
mean	0.1651	0.1645	0.1527	0.1786	0.1781	0.1707	0.3541	0.354	0.349
std	0.0139	0.0139	0.0401	0.0147	0.0147	0.0408	0.0236	0.0237	0.0598
RMSE	0.0514	0.0512	-	0.0504	0.0502	-	0.0741	0.0739	-

Table 6. Table 6 : Reduction in MSE and MAD when using θ ^ BMM subscript ^ 𝜃 BMM \hat{\theta}_{\textnormal{BMM}} instead of θ ¯ ¯ 𝜃 \overline{\theta}

	bike	demand	fires	GPUs	news	superconductors
MSE reduction	1.9%	17.3%	7.0%	0.6%	3.3%	-0.8%
MAD reduction	-0.3%	11.7%	7.0%	0.8%	4.4%	0.3%

Equations

$${\hat{\theta }}_{\text{aBMM}}=\overline{\theta }-\frac{1}{3}\frac{\sqrt{s_{\hat{\mathbf{\theta }}}^2}}{n\alpha +2}\text{accentset}\text{widehatsym}\mathrm{skew}(\hat{\mathbf{\theta }}),$$

$${\hat{\theta }}_{\text{MM}}=\text{accentset}\text{widehatsym}\mathrm{median}\left(\frac{1}{k}\sum _{i=1}^k{\hat{\theta }}_i,\dots ,\frac{1}{k}\sum _{i=n-k+1}^n{\hat{\theta }}_i\right),$$

$$|\text{accentset}\text{widehatsym}\mathrm{median}({\overline{\theta }}_1,\dots ,{\overline{\theta }}_g)-\theta |\le 6\sigma \sqrt{\frac{\log (1/\delta )}{n}}.$$

$${\hat{\theta }}_{\text{HL}}=\text{accentset}\text{widehatsym}\mathrm{median}\left(\left\{\frac{{\hat{\theta }}_i+{\hat{\theta }}_j}{2},i,j=1,\dots ,n,\text{ and }i\ne j\right\}\right).$$

$${\mathbf{p}}^{(j)}\sim \mathbb{P},j=1,\dots ,J,$$

$${\hat{\theta }}_{\mathbf{p}}=\text{accentset}\text{widehatsym}\mathrm{median}\left(\sum _{i=1}^n{p}_i^{(1)}{\hat{\theta }}_i,\dots ,\sum _{i=1}^n{p}_i^{(J)}{\hat{\theta }}_i\right).$$

$$\mathbb{E}[({\hat{\theta }}_{\mathbf{p}}-\theta {)}^2]$$

$$\mathbb{E}[({\hat{\theta }}_{\mathbf{p}}-\theta {)}^2]=\mathbb{E}[({\hat{\theta }}_{\mathbf{p}}-\overline{\theta }+\overline{\theta }-\theta {)}^2]=\mathbb{E}\left[({\hat{\theta }}_{\mathbf{p}}-\overline{\theta }{)}^2\right]+\mathbb{E}\left[(\overline{\theta }-\theta {)}^2\right]+2\mathbb{E}\left[({\hat{\theta }}_{\mathbf{p}}-\overline{\theta })(\overline{\theta }-\theta )\right],$$

$$\mathbb{E}\left[({\hat{\theta }}_{\mathbf{p}}-\overline{\theta })(\overline{\theta }-\theta )\right]$$

$$=\mathbb{E}\left[({\hat{\theta }}_{\mathbf{p}}-\mathbb{E}[{\hat{\theta }}_{\mathbf{p}}])(\overline{\theta }-\theta )\right]-\mathbb{E}\left[(\overline{\theta }-\theta {)}^2\right]$$

$$=\mathrm{Cov}({\hat{\theta }}_{\mathbf{p}},\overline{\theta })-\mathbb{V}[\overline{\theta }]$$

$$=\mathbb{V}[\overline{\theta }]\left({\rho }_{{\hat{\theta }}_{\mathbf{p}},\overline{\theta }}\sqrt{\frac{\mathbb{V}[{\hat{\theta }}_{\mathbf{p}}]}{\mathbb{V}[\overline{\theta }]}}-1\right),$$

$$\mathbb{E}\left[({\hat{\theta }}_{\mathbf{p}}-\theta {)}^2\right]=\mathbb{E}[(\overline{\theta }-\theta {)}^2]+\mathbb{E}\left[({\hat{\theta }}_{\mathbf{p}}-\overline{\theta }{)}^2\right]-2\mathbb{V}[\overline{\theta }]\left(1-{\rho }_{{\hat{\theta }}_{\mathbf{p}},\overline{\theta }}\sqrt{\frac{\mathbb{V}[{\hat{\theta }}_{\mathbf{p}}]}{\mathbb{V}[\overline{\theta }]}}\right),$$

$$\mathbb{E}\left[({\hat{\theta }}_{\mathbf{p}}-\theta {)}^2\right]$$

$$\le \mathbb{E}[({\hat{\theta }}_{\mathbf{p}}-\overline{\theta }{)}^2]-\mathbb{V}[\overline{\theta }]\left(1-2\sqrt{\frac{\mathbb{V}[{\hat{\theta }}_{\mathbf{p}}]}{\mathbb{V}[\overline{\theta }]}}\right).$$

$$\mathbf{p}$$

$$\mathbf{w}\mid \mathbf{p}$$

$$\mathbf{p}\mid \mathbf{w}$$

$${\hat{\theta }}_{\text{BMM}}=\overline{\theta }-\frac{1}{3}\frac{\sqrt{s_{\hat{\mathbf{\theta }}}^2}}{n\alpha +2}\text{accentset}\text{widehatsym}\mathrm{skew}(\hat{\mathbf{\theta }})+o\left(\frac{1}{n\alpha }\right),$$

$$\mathbb{E}[Y|\hat{\mathbf{\theta }}]=\mathbb{E}\left[\sum _{i=1}^n{p}_i{\hat{\theta }}_i\mid \hat{\mathbf{\theta }}\right]=\sum _{i=1}^n\mathbb{E}\left[p_i\right]{\hat{\theta }}_i=\sum _{i=1}^n\left(\frac{\alpha }{n\alpha }\right){\hat{\theta }}_i=\frac{1}{n}\sum _{i=1}^n{\hat{\theta }}_i=\overline{\theta },$$

$$\mathbb{V}$$

$$=\frac{1}{(n\alpha +1)}\left(\frac{1}{n}\sum _{i=1}^n{\hat{\theta }}_i-\left(\frac{1}{n}\sum _{i=1}^n{\hat{\theta }}_i\right)\left(\frac{1}{n}\sum _{i=1}^n{\hat{\theta }}_i\right)\right)=\frac{1}{n\alpha +1}\left(\overline{{\theta }^2}-{\left(\overline{\theta }\right)}^2\right)$$

$$=\frac{1}{n\alpha +1}{s}_{\hat{\mathbf{\theta }}}^2,$$

$$\mathbb{E}[Y^2|\hat{\mathbf{\theta }}]=\mathbb{V}[Y|\hat{\mathbf{\theta }}]-{\mathbb{E}}^2[Y|\hat{\mathbf{\theta }}]=\frac{\overline{{\theta }^2}+n\alpha (\overline{\theta }{)}^2}{n\alpha +1}.$$

$$\mathbb{E}[Y^m|\hat{\mathbf{\theta }}]=\sum _{k=0}^{m-1}\left(\frac{(n\alpha +k-1)!}{(n\alpha +m-1)!}\frac{(m-1)!}{k!}\mathbb{E}[Y^k|\hat{\mathbf{\theta }}]\cdot \sum _{i=1}^n\alpha {\hat{\theta }}_i^{m-k}\right).$$

$$\mathbb{E}[p_1^{{\beta }_1}{p}_2^{{\beta }_2}\cdots p_n^{{\beta }_n}]$$

$$=\frac{\Gamma (n\alpha )}{\Gamma (n\alpha +{\sum }_{i=1}^n{\beta }_i)}\prod _{i=1}^n\frac{\Gamma ({\beta }_i+\alpha )}{\Gamma (\alpha )}.$$

$$\mathbb{E}\left[{\left(\sum _{i=1}^n{p}_i{\hat{\theta }}_i\right)}^m\right]$$

$$=\sum _j{\hat{\theta }}_j\mathbb{E}\left[\sum _{k_1+\cdots +k_n=m-1}\left(\genfrac{}{}{0ex}{}{m-1}{k_1,\dots ,k_n}\right)\left(\prod _{i=1}^n{p}_i^{k_i+{\mathbb{I}}_{[i=j]}}\right)\right]$$

$$=\sum _j{\hat{\theta }}_j\left(\sum _{k_1+\cdots +k_n=m-1}\left(\genfrac{}{}{0ex}{}{m-1}{k_1,\dots ,k_n}\right)\left(\prod _{i=1}^n{\hat{\theta }}_i^{k_i}\right)\mathbb{E}\left[\prod _{i=1}^n{p}_i^{k_i+{\mathbb{I}}_{[i=j]}}\right]\right).$$

$$\mathbb{E}\left[{\left(\sum _{i=1}^n{p}_i{\hat{\theta }}_i\right)}^m\right]$$

$$=\sum _j{\hat{\theta }}_j\left(\sum _{k_1+\cdots +k_n=m-1}\left(\genfrac{}{}{0ex}{}{m-1}{k_1,\dots ,k_n}\right)\left(\prod _{i=1}^n{\hat{\theta }}_i^{k_i}\right)\frac{\alpha +k_j}{n\alpha +m-1}\mathbb{E}\left[\prod _{i=1}^n{p}_i^{k_i}\right]\right)$$

$$=\sum _j{\hat{\theta }}_j\alpha \frac{1}{n\alpha +m-1}\mathbb{E}\left[\sum _{k_1+\cdots +k_n=m-1}\left(\genfrac{}{}{0ex}{}{m-1}{k_1,\dots ,k_n}\right)\prod _{i=1}^n{\hat{\theta }}_i^{k_i}{p}_i^k\right]$$

$$+\sum _j{\hat{\theta }}_j\frac{1}{n\alpha +m-1}\mathbb{E}\left[\sum _{k_1+\cdots +k_n=m-1}{k}_j\left(\genfrac{}{}{0ex}{}{m-1}{k_1,\dots ,k_n}\right)\prod _{i=1}^n{\hat{\theta }}^{k_i}{p}_i\right].$$

$$\sum _j{\hat{\theta }}_j\alpha \frac{1}{n\alpha +m-1}\mathbb{E}\left[\sum _{k_1+\cdots +k_n=m-1}\left(\genfrac{}{}{0ex}{}{m-1}{k_1,\dots ,k_n}\right)\prod _{i=1}^n{\hat{\theta }}_i^{k_i}{p}_i^k\right]=\sum _j\alpha {\hat{\theta }}_j\frac{1}{n\alpha +m-1}\mathbb{E}[Y^{m-1}],$$

$$\mathbb{E}\left[\sum _{k_1+\cdots +k_n=m-1}{k}_j\left(\genfrac{}{}{0ex}{}{m-1}{k_1,\dots ,k_n}\right)\prod _{i=1}^n{\hat{\theta }}^{k_i}{p}_i\right]$$

$$\sum _j{\hat{\theta }}_j\frac{1}{n\alpha +m-1}\mathbb{E}\left[\sum _{k_1+\cdots +k_n=m-1}{k}_j\left(\genfrac{}{}{0ex}{}{m-1}{k_1,\dots ,k_n}\right)\prod _{i=1}^n{\hat{\theta }}^{k_i}{p}_i\right]$$

Full text

1 Introduction

The problem of combining many unbiased estimates $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$ of a parameter $\theta\in\mathbb{R}$ to form a single estimate $\hat{\theta}$ arises throughout statistics. A common way to perform aggregation when the distribution of $\hat{\theta}_{i}$ is unknown is via the sample mean, $\overline{\theta}=\frac{1}{n}\sum_{i=1}^{n}\hat{\theta}_{i}$. However, this can often lead to poor performance if the underlying distribution is very skewed or heavy-tailed.

For example, the importance sampling estimate of an integral $\theta=\int f(x)p(x)dx$ is formed by taking samples $X_{i}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}q$, letting $\hat{\theta}_{i}=f(X_{i})p(X_{i})/q(X_{i})$ and combining the estimates using the sample mean, $\frac{1}{n}\sum_{i=1}^{n}\hat{\theta}_{i}\approx\theta$. The estimator is unbiased and guaranteed to converge almost surely to $\theta$, as long as $\theta$ exists. Still, in many cases this ratio estimator has a complicated distribution with extremely high or infinite variance, making the estimates unreliable. Similar situations arise in various modern statistical procedures such as cross-validation, bagging and random forests.

If the underlying distribution of $\hat{\theta}_{i}$ were known, aggregation could be performed using maximum likelihood, which enjoys great theoretical properties. In general, both its bias and variance decrease as $O(1/n)$, so some amount of bias is introduced at the expense of variance reduction. However, since mean squared error can be decomposed as the sum of bias squared and variance, the bias component is asymptotically negligible compared to the variance. In the problem at hand, the distribution of $\hat{\theta}_{i}$ is not known, but one can still look for a non-parametric estimator that similarly introduces a small, asymptotically negligible bias, to obtain a significant variance reduction relative to the sample mean.

Given many unbiased estimators $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}\in\mathbb{R}$, consider the following procedure:

draw $\mathbf{p}^{(j)}\sim\operatorname{Dir}_{n}(\alpha,\ldots,\alpha)$ for $j=1,\ldots,J$; 2. 2.

compute $Y_{j}=\sum_{i=1}^{n}p_{i}^{(j)}\hat{\theta}_{i}$, for $j=1,\ldots,J$; 3. 3.

estimate \hat{\theta}_{\textnormal{BMM}}=\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{median}}}(Y_{1},\ldots,Y_{j}).

Here, $\operatorname{Dir}_{n}(\alpha,\ldots,\alpha)$ denotes the Dirichlet distribution with parameter vector $(\alpha,\ldots,\alpha)\in\mathbb{R}^{n}$.

This procedure, referred to as the Bayesian median of means, is detailed in Section 4. Since the median and mean of a symmetric distribution coincide, one can think of this estimator as symmetrizing $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$ via averaging before applying the median for added robustness. Section 4.2 studies the theoretical properties of this estimator, in particular proving that, when $J=O(n)$, the squared bias grows as $O(1/n^{2})$ while the variance only grows as $O(1/n)$, so asymptotically the amount of bias introduced is insignificant; finite-sample guarantees are also provided, as well as consistency results. Section 5 investigates the empirical performance of this estimator in many different settings, showing that in general the trade-off is well worth it.

The procedure only has one hyperparameter, the Dirichlet concentration level $\alpha$. When $\alpha\to 0$ and $J$ is large the Bayesian median of means approximates the sample median, whereas when $\alpha\to\infty$ it approximates the sample mean. Thus, $\alpha$ can be thought of as controlling where the procedure stands between mean and median. In Proposition 4.37, it is shown that the Bayesian median of means is first-order equivalent to a skewness-corrected sample mean, with the correction controlled by $\alpha$, suggesting an approximate, non-randomized version of the procedure above, dubbed the approximate Bayesian median of means:

[TABLE]

where $s^{2}_{\hat{\bm{\theta}}}$ is the sample variance of $\hat{\bm{\theta}}=(\hat{\theta}_{1},\ldots,\hat{\theta}_{n})$ and \mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{skew}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{skew}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{skew}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{skew}}}(\hat{\bm{\theta}}) is the sample skewness. Section 4.3 argues that $\alpha=1$ gives a good compromise between robustness and efficiency.

In many of the results below, $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$ are regarded as fixed, and the statistical procedures are analyzed conditionally. This is intended to reproduce the common setting in which the unbiased estimators $\hat{\theta}_{i}$ have been obtained before the Bayesian median of means is used. It also implies some of the analyses hold even if $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$ are not independent, which is the case for many important applications, such as cross-validation. Regarding $\hat{\theta}_{i}$ as fixed also bypasses the issue of computation, since in some cases the computational effort of drawing Dirichlet samples could be directed to obtaining further data points, as in Example 1.1 below. When obtaining more points is relatively cheap, the approximate Bayesian median of means remains an attractive alternative that does not require further sampling.

Example 1.1** (Importance sampling).**

As a first example, consider the importance sampling setting discussed above. Let $X_{1},\ldots,X_{n}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}q=\text{Expo}(1)$ be independent and identically distributed Exponential random variables, to be used in estimating the mean of $p=\text{Expo}(\lambda)$, that is, $1/\lambda$. Given a sample $X_{1},\ldots,X_{n}$, one forms the importance sampling terms $\hat{\theta}_{i}=X_{i}\cdot p(X_{i})/q(X_{i})=\lambda X_{i}e^{-(\lambda-1)X_{i}}$, and estimates $\hat{\theta}_{\text{IS}}=\frac{1}{n}\sum_{i=1}^{n}\hat{\theta}_{i}$. This sample mean will be compared to the Bayesian median of means estimator, $\hat{\theta}_{\textnormal{BMM}}$, and its approximation, $\hat{\theta}_{\textnormal{aBMM}}$, formed using the procedures outlined above.

Figure 1 shows the mean squared error (MSE) and mean absolute deviation (MAD) of both procedures, with $30$ values of $1/\lambda$, equispaced between $1.25$ and $2$. This ranges from easy settings (when $1/\lambda\approx 1.25$) to complicated ones ($1/\lambda\approx 2$). In particular, when $1/\lambda\geq 2$ the variance of $\hat{\theta}_{\text{IS}}$ becomes infinite. Here, the number of importance sampling draws is $n=1000$, the number of Dirichlet draws is $J=1000$, and $\alpha=1$. For each $\lambda$ the estimator is computed $5000$ times to estimate the errors. Note the Bayesian median of means requires sampling $J=1000$ additional Dirichlet random variables, and this extra computational allowance could alternatively have been used to generate more points $\hat{\theta}_{i}$; the approximate algorithm, on the other hand, requires no extra sampling and exhibits similar performance.

In this simulation, the MSE of the Bayesian median of means and its approximation were smaller than the importance sampling estimator’s 30 out of 30 times. The same was true for mean absolute deviation (MAD). For small values of $\lambda$ both estimators perform comparably. As $\lambda$ increases, it is easier to see the stability brought about by the Bayesian median of means makes it significantly better. In all 30 instances, the usual importance sampling estimator had larger standard deviation but smaller bias than both procedures, displaying the bias-variance trade-off involved. While considering values of $1/\lambda$ closer to 1 should favor the usual estimator, taking 30 equispaced values of $\lambda$ between $1$ and $2$ still sees better MSE and MAD performance by the Bayesian median of means 29 and 28 times; the approximate version is still better 30 and 29 times. Section 5.3 considers more extensive simulations and real data examples, reaching similar conclusions.

Paper structure. Section 2 reviews the literature and similar attempts in using the median to add robustness to mean estimates. Section 3 introduces the general idea of symmetrizing estimates before using the median to estimate location parameters, and investigates to what degree that can help decrease mean squared error. Section 4 settles on a particular distribution for symmetrizing estimates, giving rise to the Bayesian median of means. It also analyzes its theoretical properties and gives both asymptotic and finite-sample guarantees, and presents the (non-randomized) approximate Bayesian median of means. Finally, it also considers different ways of setting $\alpha$, the only hyperparameter in the procedure. Section 5 looks at the performance of these estimators in a myriad settings, including both real and simulated data, in particular comparing them to the sample mean. Finally, Section 6 concludes by giving further research directions.

2 Related Work

The idea of combining mean and median estimators has been visited several times in the statistical robustness literature, particularly for the estimation of location parameters in symmetric distributions. For instance, [Lai et al., 1983] propose an adaptive estimator that picks either the sample mean or median to estimate the center of a symmetric distribution, while [Damilano and Puig, 2004] and [Chan and He, 1994] investigate using linear combinations of mean and medians, with weights picked according to asymptotic criteria.

More recently, there has been intense work on the so-called median of means estimator (see [Alon et al., 1999], [Jerrum et al., 1986], [Nemirovsky and Yudin, 1983]). Given independent and identically distributed random variables $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$, the median of means estimator for $\mathbb{E}[\hat{\theta}_{1}]$ is given by dividing the data into blocks with $k$ elements and estimating

[TABLE]

with minor adjustments if $n/k$ is not an integer. For instance, [Devroye et al., 2016] discuss the mean estimation problem from a non-asymptotic perspective and show that the median of means estimator, among others, can outperform the sample mean in terms of concentration bounds, via the following proposition, proved in the same paper.

Proposition 2.2**.**

Consider an iid sample $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$ with mean $\theta$ and variance $\sigma^{2}$. Define the median of means estimator $\hat{\theta}_{\textnormal{MM}}$ as in (1), with $k$ elements in each of the $g=n/k$ groups. Then, with probability $1-\delta$,

[TABLE]

Such concentration cannot be achieved by the sample mean in general, unless stronger hypotheses, such as sub-Gaussianity of $\hat{\theta}_{i}$, are assumed. Variants of this estimator are further analyzed in the heavy-tailed setting of [Brownlees et al., 2015], [Hsu and Sabato, 2016] and [Bubeck et al., 2013]. The estimator was also used to combine Bayesian posterior updates in split datasets in [Minsker et al., 2014].

Unfortunately, however, there are practical challenges in using the median of means estimator. First, there is little guidance in how to pick the number of groups, which essentially amounts to the estimator’s willingness to trade off bias for variance. Furthermore, in spite of its theoretical properties, the median of means underutilizes the data available by only using each datapoint once. This guarantees independence between blocks, but limits the number of means one can obtain in a given dataset. This requirement can be relaxed to a degree, but not completely (see [Joly and Lugosi, 2016]). The estimator considered here has no such restrictions, and can be viewed as a smoothed version of the median of means. Besides, the randomness introduced in sampling the blocks allow for probabilistic analyses and parameter choices conditional on the realized datapoints. A further benefit is that, unlike the median of means, its smoothed counterpart does not depend on the order of the datapoints while still being computationally tractable.

In fact, the median of means itself can be cast as a computation compromise on the celebrated Hodges-Lehmann estimator proposed in [Hodges Jr and Lehmann, 1963]:

[TABLE]

Many theoretical properties are known about it; for instance, it has an asymptotic breakdown point of $0.29$, meaning a contamination of up to $29\%$ of the data taking arbitrary value still leaves the estimator bounded (unlike the sample mean, which has an asymptotic breakdown point of 0). However, generalizing it to $m$-averages in a straightforward way, with $m\geq 2$, is harder both in theoretical and computational terms.

Closer in spirit to the present work is the suggestion in [Bühlmann, 2003] to compute sample means using bootstrap samples and then aggregate them via the median. This is similar to bagging estimators $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$, but using the median instead of the mean in assembling the averages. In particular, [Bühlmann, 2003] empirically observed that the median helped make the final estimators better in terms of MSE in non-convex problems. This work generalizes [Bühlmann, 2003] by resampling the weights used in the averaging using a Dirichlet distribution and by providing a more extensive theoretical analysis, in particular, deriving a non-randomized approximation. This leads to a smoother estimator, particularly relevant in the heavy-tailed or skewed setting, and also allows for theoretical results and hyperparameter recommendations based on the literature on Dirichlet means (see [Cifarelli and Regazzini, 1990], [Cifarelli and Melilli, 2000], [Regazzini et al., 2002] and [Diaconis and Kemperman, 1996]) and, more generally, on $P$-means ([Pitman, 2018]).

3 Median of Weighted Means

Given a probability space $(\Omega,\mathcal{F},\mathbb{P}_{\hat{\theta}})$, let $\hat{\theta}_{i}:\Omega\to\mathbb{R}$, $i=1,\ldots,n$ be a collection of independent and identically distributed random variables, with $\mathbb{E}[\hat{\theta}_{1}]=\theta$ and $\mathbb{V}[\hat{\theta}_{1}]=\sigma^{2}<\infty$. This will be denoted below by $\hat{\theta}_{i}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}[\theta,\sigma^{2}]$. Consider the problem of obtaining an estimator for $\theta$ given observations $\hat{\bm{\theta}}=(\hat{\theta}_{1},\ldots,\hat{\theta}_{n})$.

A common aggregation procedure is the sample mean, $\overline{\theta}=\frac{1}{n}\sum_{i=1}^{n}\hat{\theta}_{i}$. Besides retaining unbiasedness, it also possesses many satisfying theoretical properties; for example, the sample mean is the best linear unbiased estimator for the population mean, minimizes the maximum expected mean squared error populations with bounded variance (see [Bickel and Lehmann, 1981]), is an admissible procedure in one dimension ([Joshi, 1968]), and is the maximum likelihood estimate in exponential families under independent sampling. For these reasons, as well as computationally simplicity, the sample mean is a widely adopted non-parametric procedure for aggregating one-dimensional estimators.

However, in many settings the underlying distribution of $\hat{\theta_{1}}$ is very skewed or heavy tailed, in which case the sample mean becomes highly unstable. A robust solution to the aggregation problem is to use the sample median. If the distribution is symmetric, then mean and median coincide, the sample median is still an unbiased estimator, and the median can have better asymptotic mean squared error than the bias, as in Example 3.3.

Example 3.3** (Mean and median).**

Let $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}N(\theta,\sigma^{2})$, and consider estimators $\overline{\theta}=\frac{1}{n}\sum_{i=1}^{n}\hat{\theta}_{i}$ and \mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{median}}}(\hat{\theta}_{1},\ldots,\hat{\theta}_{n}), which are unbiased for estimating $\theta$. In this case, the asymptotic mean squared error is given by the asymptotic variance. Since $n\cdot\mathbb{V}[\overline{\theta}]\to\sigma^{2}$ as $n\to\infty$, and also n\cdot\mathbb{V}[\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{median}}}(\hat{\theta}_{1},\ldots,\hat{\theta}_{n})]\to\pi\sigma^{2}/2 (see Proposition 4.16 below), it holds that, asymptotically, \operatorname{MSE}(\overline{\theta})=(2/\pi)\cdot\operatorname{MSE}(\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{median}}}(\hat{\theta}_{1},\ldots,\hat{\theta}_{n})), so the sample mean exhibits better asymptotic performance.

On the other hand, consider $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}\text{Laplace}(\theta,\sigma^{2})$. As before, both $\overline{\theta}=\frac{1}{n}\sum_{i=1}^{n}\hat{\theta}_{i}$ and \mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{median}}}(\hat{\theta}_{1},\ldots,\hat{\theta}_{n}) are unbiased for estimating $\theta$. In particular, note \mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{median}}}(\hat{\theta}_{1},\ldots,\hat{\theta}_{n}) is the maximum likelihood estimator. When $n\to\infty$, $n\cdot\mathbb{V}[\overline{\theta}]\to\sigma^{2}$ while n\cdot\mathbb{V}[\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{median}}}(\hat{\theta}_{1},\ldots,\hat{\theta}_{n})]\to\sigma^{2}/2, so \operatorname{MSE}(\overline{\theta})=2\cdot\operatorname{MSE}(\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{median}}}(\hat{\theta}_{1},\ldots,\hat{\theta}_{n})), and the median is asymptotically better than the sample mean. In more extreme cases, say if $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}\text{Cauchy}(0,1)$, the median can be asymptotically infinitely better than the sample mean.

One can try to use a compromise between the two estimators, for instance the Hodges-Lehmann estimator, $\hat{\theta}_{\textnormal{HL}}$, defined in (2), which is still unbiased for symmetric distributions. Asymptotically, with a Normal sample, $\operatorname{MSE}(\overline{\theta})=(3/\pi)\cdot\operatorname{MSE}(\hat{\theta}_{\textnormal{HL}})$, and with a Laplace sample, $\operatorname{MSE}(\overline{\theta})=1.5\cdot\operatorname{MSE}(\hat{\theta}_{\textnormal{HL}})$. In fact, if $\mathcal{S}$ is the set of symmetric distributions centered at $\theta$, it can be shown that asymptotically $\inf_{\{F_{\hat{\bm{\theta}}}:F_{\hat{\bm{\theta}}}\in\mathcal{S}\}}\operatorname{MSE}(\overline{\theta})/\operatorname{MSE}(\hat{\theta}_{\textnormal{HL}})=0.864$ (see [Hodges et al., 1956]), so $\hat{\theta}_{\textnormal{HL}}$ never fares much worse than $\overline{\theta}$ but can sometimes do much better.

In general, however, the median will be a biased estimator. For instance, if the distribution of $\hat{\theta}_{i}$ is highly skewed, \mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{median}}}(\hat{\theta}_{i}) might be very different from $\theta=\mathbb{E}[\hat{\theta}_{i}]$. One way to soften the bias is to take the median of weighted averages of $\hat{\theta}_{i}$, which are more symmetric around $\theta$ by the Central Limit Theorem. Hence, the application of the median becomes less costly in terms of bias, while still guaranteeing increased robustness.

A general scheme for aggregating many unbiased estimators $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$ using a vector of probabilities $\mathbf{p}=(p_{1},\ldots,p_{n})$, with $\sum_{i=1}^{n}p_{i}=1$, $p_{i}\geq 0$, is as follows:

Sample:

[TABLE]

where $\mathbb{P}$ is a probability measure such that $\sum_{i=1}^{n}p_{i}^{(j)}=1$. 2. 2.

Estimate:

[TABLE]

Consider the following choices for $\mathbb{P}$:

• •

if $\mathbb{P}$ sets $p_{i}=1$ for $i$ chosen uniformly at random and $J$ is sufficiently large, then $\hat{\theta}_{\mathbf{p}}$ is essentially the sample median;

• •

if $\mathbb{P}$ is a delta mass at $(\frac{1}{n},\ldots,\frac{1}{n})$, then $\hat{\theta}_{\mathbf{p}}$ is the sample mean;

• •

if $\mathbb{P}$ selects sets $S$ of size $k$ in a uniformly chosen partition of $\{1,\ldots,n\}$ and sets $p_{i}=\frac{1}{k}$ for $i\in S$ and [math] otherwise, then $\hat{\theta}_{\mathbf{p}}$ is a randomly-partitioned median of means estimator;

• •

if $\mathbb{P}$ is $\operatorname{Dir}_{n}(\alpha,\ldots,\alpha)$, then the resulting estimator is the Bayesian median of means, $\hat{\theta}_{\textnormal{BMM}}$.

Hence, estimators with randomized weights give a way to interpolate between the sample mean, with low bias and possibly high variance, and the sample median, with low variance but possibly high bias. See Figure 2. The two extreme choices of $\mathbb{P}$ leading to the sample mean and sample median come from degenerate Dirichlet distributions placing all mass at the centroid of the simplex ($\alpha\to 0$), or splitting the mass equally among the vertices of the simplex ($\alpha\to\infty$). More general distributions over the simplex are possible, and many of the results in this paper can readily be extended to that case (see [Newton and Raftery, 1994] and [Pitman, 2018]). In particular, distributions can be picked to encode any prior information available about the sample $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$ such as skewness or symmetry.

The reduction in variance achieved by the median of weighted means can have a drastic effect on mean squared error when $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$ come from a distribution with high variance. This can be understood through the following proposition.

Proposition 3.4**.**

Let $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$ be iid, unbiased estimates of a parameter $\theta\in\mathbb{R}$. Let $\overline{\theta}=\frac{1}{n}\sum_{i=1}^{n}\hat{\theta}_{i}$ be the sample mean, and $\hat{\theta}_{\mathbf{p}}$ be any median of weighted means estimator (4). Then the mean squared error of $\hat{\theta}_{\mathbf{p}}$ can be bounded by

[TABLE]

where the expectation is over both the data $\{\hat{\theta}_{i}\}_{i=1}^{n}$ and the weights $\{\mathbf{p}^{(j)}\}_{j=1}^{J}$.

Proof 3.5**.**

Decompose the expectation as

[TABLE]

and note the cross-term can be written

[TABLE]

where $\rho_{\hat{\theta}_{\mathbf{p}},\overline{\theta}}$ is the correlation between $\overline{\theta}$ and $\hat{\theta}_{\mathbf{p}}$. Putting it together,

[TABLE]

and since $\overline{\theta}$ is unbiased, $\mathbb{E}[(\overline{\theta}-\theta)^{2}]=\mathbb{V}[\overline{\theta}]$, so

[TABLE]

Thus, the mean squared error of $\hat{\theta}_{\mathbf{p}}$ is given by a term measuring the discrepancy between $\hat{\theta}_{\mathbf{p}}$ and $\overline{\theta}$, minus a term measuring the variance reduction achieved by $\hat{\theta}_{\mathbf{p}}$ with respect to $\overline{\theta}$. If $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$ are coming from a distribution with high variance, then generally $\mathbb{V}[\hat{\theta}_{\mathbf{p}}]\ll\mathbb{V}[\overline{\theta}]$, and the second term in the right-hand side of the bound is negative and very large. If, on the other hand, $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$ are coming from a distribution with low variance, then $\hat{\theta}_{\mathbf{p}}\approx\overline{\theta}$, and there should be no significant differences between the mean squared errors of the two estimators.

4 Bayesian Median of Means

Now consider the median of weighted means estimator obtained by sampling the probabilities in (3) from a $\operatorname{Dir}_{n}(\alpha,\ldots,\alpha)$ distribution. Recall this scheme is broad enough to interpolate between median and mean, while still being analytically tractable, and is given by

draw $\mathbf{p}^{(j)}\sim\operatorname{Dir}_{n}(\alpha,\ldots,\alpha)$ for $j=1,\ldots,J$; 2. 2.

compute $Y_{j}=\sum_{i=1}^{n}p_{i}^{(j)}\hat{\theta}_{i}$, for $j=1,\ldots,J$; 3. 3.

estimate \hat{\theta}_{\textnormal{BMM}}=\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{median}}}(Y_{1},\ldots,Y_{j}).

This estimator is ‘Bayesian’ in the sense that the probabilities $\mathbf{p}^{(j)}$ are being generated according to the Bayesian bootstrap. Indeed, consider weights $\mathbf{w}=(w_{1},\ldots,w_{n})$ with $\sum_{i=1}^{n}w_{i}=n$, $w_{i}\geq 0$ and $\mathbf{p}=\frac{1}{n}\mathbf{w}$, and assume the following underlying model

[TABLE]

so the posterior distribution is

[TABLE]

With a non-informative prior, $\gamma\to 0$, the posterior distribution becomes $\operatorname{Dir}_{n}(m/n,\ldots,m/n)$, which amounts to step 1 above with $\alpha=m/n$. In particular, if $m=n$, it is $\operatorname{Dir}_{n}(1,\ldots,1)$.

This gives a posterior on the sums $\sum_{i=1}^{n}p_{i}^{(j)}\hat{\theta}_{i}$ for a randomly sampled probability vector $\mathbf{p}^{(j)}$ and fixed $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$. Summarizing the posterior distribution by minimizing the $l_{1}$ loss for robustness yields the Bayesian median of means.

Note the usual bootstrap would sample $\mathbf{p}^{(j)}\sim\text{Mult}_{n}(m,(1/n,\ldots,1/n))$. This has the same mean as $\operatorname{Dir}_{n}(m/n,\ldots,m/n)$, and nearly the same variance. The main reason for using the Dirichlet distribution is that it confers additional smoothness to the estimator that are important in establishing theoretical results, in particular asymptotic expansions (see Section 4.2.3).

How can this estimator improve on the sample mean? Proposition 4.37 below shows that, under some regularity assumption and assuming $J=O(n)$, one can approximate,

[TABLE]

where $s^{2}_{\hat{\bm{\theta}}}=\frac{1}{n}\sum_{i=1}^{n}(\hat{\theta}_{i}-\overline{\theta})^{2}$ and \mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{skew}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{skew}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{skew}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{skew}}}(\hat{\bm{\theta}})=\frac{1}{(s^{2}_{\hat{\bm{\theta}}})^{3/2}}\cdot\frac{1}{n}\sum_{i=1}^{n}(\hat{\theta_{i}}-\overline{\theta})^{3}, so to first-order the Bayesian median of means is a corrected sample mean, with the correction proportional to the sample standard deviation and skewness. This is reminiscent of a shrinkage-type estimator.

Indeed, if $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$ are coming from a symmetric distribution, then $\hat{\theta}_{\textnormal{BMM}}\approx\overline{\theta}$, so consider the case where the underlying distribution is very right-skewed, as in the first plot in Figure 3. Assume two samples of $n$ points are obtained, as represented in red and blue in the figure. Note the blue sample happens to include a very large, but unlikely, sample point.

The second plot in Figure 3 shows the result of applying both the sample mean and the Bayesian median of means to the red and blue samples. Since the red sample is virtually symmetric, \mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{skew}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{skew}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{skew}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{skew}}}(\hat{\bm{\theta}})\approx 0 and so $\overline{\theta}\approx\hat{\theta}_{\textnormal{BMM}}$. For the blue sample, the large sample point means the sample mean overestimates $\theta$ by a lot; in this case, both $s^{2}_{\hat{\bm{\theta}}}$ and \mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{skew}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{skew}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{skew}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{skew}}}(\hat{\bm{\theta}}) are large and positive, so $\hat{\theta}_{\textnormal{BMM}}<\overline{\theta}$, as represented in Figure 3.

This schematic example explains the mechanics behind the Bayesian median of means: it introduces some bias for samples with large sample variance and skewness, and in doing so greatly reduces variance. When the underlying distribution of $\hat{\theta}_{i}$ has small variance or is symmetric, $\overline{\theta}$ and $\hat{\theta}_{\textnormal{BMM}}$ are virtually indistinguishable; however, when the underlying distribution has heavy tails or is very skewed, the estimators disagree (this phenomena was also observed in Example 1.1). Hence, the Bayesian median of means matches the sample mean performance for samples coming from approximately Normal distributions, but incurs in corrections once this is no longer the case.

To establish the theoretical properties of this estimator, Section 4.1 first looks at the distribution of the resampled averages $Y_{j}$, $j=1,\ldots,J$. Section 4.2 establishes finite-sample guarantees for the Bayesian median of means, as well as asymptotic approximations. Finally, Section 4.3 considers the issue of picking a value for the hyperparameter $\alpha$.

4.1 Conditional Moments and Density for $Y$

To understand the behavior of the Bayesian median of means, it is first necessary to study $Y_{j}$. In particular, it will be important to characterize the moments and distribution of $Y$ (the subscript $j$ will be dropped when the meaning is clear). Calculating the first two conditional moments of $Y$ is straightforward:

[TABLE]

and

[TABLE]

where $s_{\hat{\bm{\theta}}}^{2}=\frac{1}{n}\sum_{i=1}^{n}(\hat{\theta}_{i}-\overline{\theta})^{2}$, so

[TABLE]

Higher conditional moments can be found in a recursive fashion.

Lemma 4.6**.**

Let $\mathbf{p}\sim\operatorname{Dir}_{n}(\alpha,\ldots,\alpha)$ and $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$ be a set of fixed of estimates. The linear combination of Dirichlet components $Y=\sum_{i=1}^{n}p_{i}\hat{\theta}_{i}$ has moments given by the recursion

[TABLE]

Proof 4.7**.**

First, recall that if $\mathbf{p}\sim\operatorname{Dir}_{n}(\alpha,\ldots,\alpha)$, then

[TABLE]

Using $\mathbb{E}[\cdot]$ below to denote the expectation conditional on $\hat{\bm{\theta}}$, write

[TABLE]

Using (12), this becomes

[TABLE]

The first term in the last equality above is just

[TABLE]

while the expectation in the second term is

[TABLE]

so the second term becomes

[TABLE]

where the last equality follows by applying the argument above with $m-1$ instead of $m$. Putting it all together,

[TABLE]

Proceeding with the inductive argument, this gives

[TABLE]

as desired.

For example, the lemma above gives

[TABLE]

and it is not hard to compute higher-order moments as needed. From this, one can also obtain the unconditional moments: using the Law of Iterated Expectations,

[TABLE]

so $Y$ is unbiased, and by the Law of Total Variance,

[TABLE]

Note in particular that $\mathbb{V}[Y]\stackrel{{\scriptstyle\alpha\to 0}}{{\longrightarrow}}\sigma^{2}=\mathbb{V}[\hat{\theta}_{1}]$ and $\mathbb{V}[Y]\stackrel{{\scriptstyle\alpha\to\infty}}{{\longrightarrow}}\sigma^{2}/n=\mathbb{V}[\overline{\theta}]$, as expected.

Consider now the conditional distribution of $Y$. Since it is a random mean of the $\hat{\theta}_{i}$, it admits an asymptotic Normal distribution concentrated around $\overline{\theta}$.

Proposition 4.8**.**

Let $\hat{\theta}_{1},\ldots,\hat{\theta}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}[\theta,\sigma^{2}]$, suppose $\mathbf{p}\sim\operatorname{Dir}_{n}(\alpha,\ldots,\alpha)$, with $\alpha>0$, and $Y=\sum_{i=1}^{n}p_{i}\hat{\theta}_{i}$. If $\mathbb{V}[\hat{\theta}_{i}]=\sigma^{2}<\infty$, then, for almost all sequences $\hat{\theta}_{1},\hat{\theta}_{2},\ldots$, the conditional distribution of $\sqrt{n}(Y-\overline{\theta})$ converges to $N(0,\sigma^{2}/\alpha)$.

Proof 4.9**.**

Since $\mathbf{p}\sim\operatorname{Dir}_{n}(\alpha,\ldots,\alpha)$, each coordinate $p_{i}$ can be written $p_{i}\stackrel{{\scriptstyle d}}{{=}}G_{i}/\sum_{j=1}^{n}G_{j}$, where $G_{j}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}\text{Gamma}(\alpha,1)$. Thus,

[TABLE]

By the Strong Law of Large Numbers, $\frac{1}{\alpha n}\sum_{j=1}^{n}G_{j}\stackrel{{\scriptstyle as}}{{\longrightarrow}}1$. Conditioned on $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$, $\sum_{i=1}^{n}G_{i}(\hat{\theta}_{i}-\overline{\theta})$ is a weighted sum of the $G_{i}$. Note

[TABLE]

To obtain a Central Limit Theorem, it suffices to see Lindeberg’s condition is implied by

[TABLE]

and so $\frac{1}{\alpha\sqrt{n}}\sum_{i=1}^{n}G_{i}(\hat{\theta_{i}}-\overline{\theta})\Rightarrow N(0,\frac{\sigma^{2}}{\alpha})$ conditioned on $\hat{\theta}_{1},\hat{\theta}_{2},\ldots$, almost surely. An application of Slutsky’s theorem then yields the proposition.

Computing the exact distribution of $Y_{j}$, however, is more intricate. It can be done using the theory of Dirichlet means first developed by Cifarelli and Regazzini ([Cifarelli and Regazzini, 1990], [Cifarelli and Regazzini, 1994]). Through a connection between the Stieltjes transform of the distribution function of $Y$ and the Laplace transform of a related Gamma process, known as the Markov-Krein identity, they were able to obtain the cumulative distribution function for $Y$ in [Cifarelli and Regazzini, 1994].

Proposition 4.10**.**

Consider $\mathbf{p}\sim\operatorname{Dir}_{n}(\alpha,\ldots,\alpha)$, and take $Y=\sum_{i=1}^{n}p_{i}\hat{\theta}_{i}$, with $\hat{\theta}_{i}$ fixed for $i=1,\ldots,n$. The cumulative distribution function of $Y$, denoted $F_{Y|\hat{\bm{\theta}}}$, is supported on $[\min_{i}\hat{\theta}_{i},\max_{i}\hat{\theta}_{i}]$, degenerate at points $y=\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$, and absolutely continuous with respect to the Lebesgue measure as long as $\min_{i}\hat{\theta}_{i}<\max_{i}\hat{\theta}_{i}$. For $y\neq\hat{\theta}_{i}$, the probability density function of $Y$ given $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$, denoted $f_{Y|\hat{\bm{\theta}}}(y)$, can be written as:

1. (i)

if $\alpha>1/n$,

[TABLE]

furthermore, if $\alpha<1$,

[TABLE]

where $\operatorname{Im}\left\{\cdot\right\}$ denotes the imaginary part of the complex number, and $\hat{\theta}_{(i)}$ the $i$-th smallest value in $\{\hat{\theta}_{1},\ldots,\hat{\theta}_{n}\}$; 2. (ii)

if $\alpha=1/n$,

[TABLE] 3. (iii)

if $1/n>\alpha>0$,

[TABLE]

Proof 4.11**.**

Degeneracy of the cumulative distribution function is considered in Theorem 1 of [Cifarelli and Regazzini, 1994], while absolute continuity with respect to Lebesgue measure is established in Proposition 4 of [Regazzini et al., 2002]. The density formulas are given as Proposition 9 of [Regazzini et al., 2002], specialized here to the case where the underlying measure of the Dirichlet process $\tilde{\alpha}$ is supported on $\{\hat{\theta}_{1},\ldots,\hat{\theta}_{n}\}$ with $\tilde{\alpha}(\hat{\theta}_{i})=\alpha$ for $i=1,\ldots,n$. The general proof can be found in [Regazzini et al., 2000].

Example 4.12** (Behavior at the median).**

Consider the behavior of $Y|\hat{\bm{\theta}}$ near the median, $m=\operatorname{median}(Y|\hat{\bm{\theta}})$, when $\alpha=1/n$, in which case its density can be given in explicit form. From Proposition 4.10,

[TABLE]

where $C>0$ is a positive constant. Using the Strong Law of Large Numbers, as $n\to\infty$,

[TABLE]

where $F_{\hat{\bm{\theta}}}$ denotes the law of $\hat{\theta}$, and the last equality uses the fact that $m\stackrel{{\scriptstyle n\to\infty}}{{\longrightarrow}}\theta$, shown in Proposition 4.29. Thus,

[TABLE]

giving a way to specify, at least asymptotically, $f_{Y|\hat{\bm{\theta}}}(m)$.

While Proposition 4.10 determines the conditional distribution $Y|\hat{\bm{\theta}}$, directly relating the unconditional distribution of $Y_{j}$ to $\hat{\bm{\theta}}$ requires different techniques. It can be done, for example, via the Laplace transform of $G_{n\alpha}Y_{j}$, with $G_{n\alpha}\sim\text{Gamma}(n\alpha)$ independent of $Y_{j}$, which also uniquely determines the distribution of $Y_{j}$:

[TABLE]

This is proved as Proposition 29 in [Pitman, 2018], dating back to [Von Neumann, 1941] and [Watson, 1956]. The following is one of the few examples where the distribution of both $\hat{\theta}_{i}$ and $Y_{j}$ are known.

Example 4.13** (Beta distribution).**

Suppose $\hat{\theta}_{i}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}\text{Beta}(a,b)$, and $\mathbf{p}^{(j)}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}\operatorname{Dir}_{n}(a+b,\ldots,a+b)$, so $\alpha=a+b$. In this case, $\hat{\theta}_{i}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}\text{Beta}(a,b)$ is also defined by its generalized Stieltjes transform $\mathbb{E}[(1-\lambda\hat{\theta}_{i})^{-(a+b)}]=(1-\lambda)^{-a}$. Hence, (20) gives

[TABLE]

which in turn implies, by the Stieltjes transform characterization, that $Y_{j}\sim\text{Beta}(na,nb)$. In particular,

[TABLE]

agreeing with equations (16) and (19) for $\hat{\theta}_{i}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}\text{Beta}(a,b)$. Also,

[TABLE]

so the $Y_{j}$ are more concentrated and symmetrized around its mean than the original sample $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$. See Figure 4. Further, note the goal is to estimate $\theta=\mathbb{E}[\hat{\theta}_{i}]=a/(a+b)$, while the Bayesian median of means is estimating

[TABLE]

where this approximation for the median of a Beta distribution is valid for $a,b\geq 1/n$ (see [Kerman, 2011]). Hence, the bias introduced in using $\hat{\theta}_{\textnormal{BMM}}$ to estimate $\theta$ is fairly small.

The previous example also sheds light on the effects the Bayesian median of means procedure has on symmetric and unimodal distributions. In particular, if $\hat{\theta}_{i}$ comes from a symmetric distribution, then $Y$ is also symmetrically distributed (so the median adds no additional bias).

Proposition 4.14**.**

Let $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}F_{\hat{\bm{\theta}}}$, and assume $F_{\hat{\bm{\theta}}}$ is symmetric, with density $f_{\hat{\bm{\theta}}}$, mean $\theta$ and variance $\sigma^{2}$. If $\mathbf{p}\sim\text{Dir}(\alpha,\ldots,\alpha)$ and $Y=\sum_{i=1}^{n}p_{i}\hat{\theta}_{i}$, then the distribution of $Y$ is also symmetric. If in addition $F_{\hat{\bm{\theta}}}$ is unimodal, then so is the distribution of $Y$. The converse is not true.

Proof 4.15**.**

First, assume $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$ is coming from a symmetric distribution, and without loss of generality take the center of symmetry to be [math]. Then, for any constant $p_{i}$, $p_{i}\hat{\theta}_{i}\stackrel{{\scriptstyle d}}{{=}}-p_{i}\hat{\theta}_{i}$, and because $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$ are iid, $\sum_{i=1}^{n}p_{i}\hat{\theta}_{i}\stackrel{{\scriptstyle d}}{{=}}-\sum_{i=1}^{n}p_{i}\hat{\theta}_{i}$. Thus,

[TABLE]

so $Y$ is also symmetric.

Now, assume additionally that $F_{\hat{\theta}}(t)$ is unimodal. Conditional on $\mathbf{p}$, the sum of independent, symmetric and unimodal random variables is still symmetric and unimodal (see [Purkayastha, 1998]), so $Y$ itself is also symmetric and unimodal. The unconditional case follows by taking expectations as above.

On the other hand, note that by Example 4.13, a multimodal distribution such as $\hat{\theta}_{i}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}\text{Beta}(\frac{1}{2},\frac{1}{2})$ leads to $Y\sim\text{Beta}(\frac{n}{2},\frac{n}{2})$, which is unimodal if $n>2$.

4.2 Theoretical Guarantees

This section establishes the main theoretical results in the paper. Section 4.2.1 characterizes the convergence of $\hat{\theta}_{\textnormal{BMM}}$ to $\operatorname{median}(Y|\hat{\theta}_{1},\ldots,\hat{\theta}_{n})$, as $J\to\infty$ and for fixed $n$. Section 4.2.2 then bounds the error in estimating $\theta$ using $\operatorname{median}(Y|\hat{\theta}_{1},\ldots,\hat{\theta}_{n})$ in finite samples. Section 4.2.3 considers asymptotic guarantees.

4.2.1 Approximating $\widehat{\theta}_{\text{BMM}}$ with $\operatorname{median}(Y\mid\widehat{\theta})$

The limiting distribution of $\hat{\theta}_{\textnormal{BMM}}$ as $J\to\infty$ can be found via the Central Limit Theorem for medians and Proposition 4.10.

Proposition 4.16**.**

Let $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}\in\mathbb{R}$ be fixed, and consider $Y_{j}=\sum_{i=1}^{n}p^{(j)}_{i}\hat{\theta}_{i}$, where $\mathbf{p}^{(j)}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}\operatorname{Dir}_{n}(\alpha,\ldots,\alpha)$ for $j=1,\ldots,J$. If \hat{\theta}_{\textnormal{BMM}}^{(J)}=\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{median}}}(Y_{1},\ldots,Y_{J}|\hat{\bm{\theta}}) is the sample median and $m=\operatorname{median}(Y|\hat{\bm{\theta}})$ is the population median, with $m\neq\hat{\theta}_{i}$, $i=1,\ldots,n$ and $f_{Y|\hat{\bm{\theta}}}(m)>0$, then, as $J\to\infty$,

[TABLE]

with $f_{Y|\hat{\bm{\theta}}}(y)$ given in Proposition 4.10.

Proof 4.17**.**

Since conditional on $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$, $\left\{Y_{j}\right\}_{j=1}^{J}$ are iid, defining $W_{J}=\frac{1}{J}\sum_{j=1}^{J}\mathbb{I}_{[Y_{j}\leq m]}$, the Central Limit Theorem gives $\sqrt{J}(W_{J}-1/2)\Rightarrow N(0,1/4)$. Considering the generalized inverse $F_{Y|\hat{\bm{\theta}}}^{-1}(t)$ and recalling $\frac{d}{dt}F_{Y|\hat{\bm{\theta}}}^{-1}(t)=\frac{1}{f_{Y|\hat{\bm{\theta}}}(F^{-1}_{Y|\hat{\bm{\theta}}}(t))}$, a straightforward application of the Delta Method yields $\sqrt{J}(\hat{\theta}_{\textnormal{BMM}}^{(J)}-m)\Rightarrow N(0,1/(4f^{2}_{Y|\hat{\bm{\theta}}}(m)))$, where $f_{Y|\hat{\bm{\theta}}}(y)$ is defined in Proposition 4.10.

It is possible to obtain better control over the fluctuations between $\hat{\theta}_{\textnormal{BMM}}$ and $\operatorname{median}(Y|\hat{\bm{\theta}})$ via finite-sample concentration bounds. In particular, this suggests how to set $J$ as a function of $n$ to obtain concentration around the median for any small $t=O(1/n)$.

Proposition 4.18**.**

Consider $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}\in\mathbb{R}$, and let $Y_{j}=\sum_{i=1}^{n}p^{(j)}_{i}\hat{\theta}_{i}$, where $\mathbf{p}^{(j)}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}\operatorname{Dir}_{n}(\alpha,\ldots,\alpha)$ for $j=1,\ldots,J$. Let \hat{\theta}_{\textnormal{BMM}}=\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{median}}}(Y_{1},\ldots,Y_{J}|\hat{\bm{\theta}}) and assume $m=\operatorname{median}(Y|\hat{\bm{\theta}})$ is unique. For $t$ small enough that $\min_{[m,m+t]}f_{Y|\hat{\bm{\theta}}}(\xi)>C>0$ for some $C=C(\hat{\theta}_{1},\ldots,\hat{\theta}_{n})$, it holds that

[TABLE]

so $\hat{\theta}_{\textnormal{BMM}}$ concentrates exponentially fast.

Proof 4.19**.**

To establish concentration, split

[TABLE]

and bound each term in turn. For the first, note that conditionally on $\hat{\bm{\theta}}$ the $Y_{1},\ldots,Y_{J}$ are independent, so the sample median $\hat{\theta}_{\textnormal{BMM}}$ is bigger than $m+t$ if at least half of the $J$ points is so. Thus,

[TABLE]

where $q_{+}(t)=P_{\mathbf{p}}\left[Y\geq m+t\mid\hat{\bm{\theta}}\right]<\frac{1}{2}$. The Hoeffding bound gives

[TABLE]

Note the density of $Y$, $f_{Y|\hat{\bm{\theta}}}(y)$ is strictly positive in $[m,m+t]$ since the support of $Y|\hat{\bm{\theta}}$ is the convex hull of the $\hat{\theta}_{i}$. Thus, there exists $C>0$ such that $\min_{[m,m+t]}f_{Y|\hat{\bm{\theta}}}(\xi)>C$, and

[TABLE]

so it is possible to conclude

[TABLE]

Bounding the other term in (21) similarly yields the proposition.

In some cases, it is possible to get an asymptotic description of the constant $C$ above.

Example 4.20** (Value of $C$).**

Recall from Example 4.12 that when $\alpha=1/n$, if $F_{\hat{\theta}}$ denotes the cdf of $\hat{\theta}_{i}$ and $\tilde{C}=\frac{1}{\pi}\sin\left(\frac{\pi}{n}\sum_{i=1}^{n}\mathbb{I}_{[m>\hat{\theta}_{i}]}\right)$, then, as $n\to\infty$, the density at the median $m$ converges to

[TABLE]

so the bound in (22) becomes

[TABLE]

The bound in Proposition 4.18 works well for small values of $t$. For example, it ensures that $J=O(n^{2})$ suffices to guarantee that, with constant probability, $|\hat{\theta}_{\textnormal{BMM}}-m|$ is $O(\frac{1}{n})$ as $n\to\infty$. For large values of $t$, a different bound can yield better results.

Proposition 4.21**.**

Let $Y_{1},\ldots,Y_{J}$ be iid given $\hat{\bm{\theta}}=(\hat{\theta}_{1},\ldots,\hat{\theta}_{n})$, and assume $m=\operatorname{median}(Y|\hat{\bm{\theta}})$ is unique. For $t>2\sqrt{\mathbb{V}[Y|\hat{\bm{\theta}}]}$,

[TABLE]

Proof 4.22**.**

From the definition of conditional median and Jensen’s inequality,

[TABLE]

so, from Markov’s inequality,

[TABLE]

Now, if the sample median deviates from the conditional median, $m$, by more than $t$ then at least $J/2$ of the $Y_{j}$ are more that than $t$ apart from $m$. Since $Y_{1},\ldots,Y_{J}$ are conditionally iid, this implies

[TABLE]

To bound this probability, recall the Chernoff bound for $X\sim\text{Bin}(J,p)$ yields, for $0\leq a\leq J(1-p)$,

[TABLE]

Thus, for $t>2\sqrt{\mathbb{V}[Y|\hat{\bm{\theta}}]}$,

[TABLE]

which gives suitable control over $|\hat{\theta}_{\textnormal{BMM}}-m|$ for large $t$.

Both Propositions 4.18 and 4.21 are useful in quantifying how big $J$ should be to control $|\hat{\theta}_{\textnormal{BMM}}-m|$. An alternative perspective on how large $J$ should be comes from analyzing the asymptotic bias and variance incurred by $\hat{\theta}_{\textnormal{BMM}}$ in estimating $m$. Both decrease as $O(\frac{1}{J})$, as the next proposition shows.

Proposition 4.23**.**

Let $Y_{1},\ldots,Y_{J}$ be conditionally iid and $m=\operatorname{median}(Y|\hat{\bm{\theta}})$ unique, and define \hat{\theta}_{\textnormal{BMM}}=\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{median}}}(Y_{1},\ldots,Y_{J}|\hat{\bm{\theta}}). Then:

[TABLE]

In particular both the bias of $\hat{\theta}_{\textnormal{BMM}}$ (in estimating $m$) and the variance decrease as $O(1/J)$.

Proof 4.24**.**

Assume for simplicity $J=2k+1$ so $\hat{\theta}_{\textnormal{BMM}}=Y_{(k)}$. Thus,

[TABLE]

The idea will be to approximate this integral via a Laplace approximation. First, using Stirling’s formula,

[TABLE]

Since $f_{Y|\hat{\bm{\theta}}}(m)>0$, $\log F_{Y|\hat{\bm{\theta}}}(t)-\log\left(1-F_{Y|\hat{\bm{\theta}}}(t)\right)$ is maximized at the conditional median, $m$. Taking $r=1$, the Laplace approximation yields

[TABLE]

Putting it together,

[TABLE]

so $\mathbb{E}[\hat{\theta}_{\textnormal{BMM}}]-m=O(1/k)$, and the bias of $\hat{\theta}_{\textnormal{BMM}}$ in estimating $m$ decreases as $O(1/J)$. In fact, a more careful approximation yields

[TABLE]

Using $r=2$ in (23) leds to the following expression for the variance:

[TABLE]

Note this matches the asymptotic variance found in Proposition 4.16. Hence, the variance of $\hat{\theta}_{\textnormal{BMM}}$ also decreases as $O\left(\frac{1}{J}\right)$.

In terms of variance, note that $\mathbb{V}[\overline{\theta}]=\sigma^{2}/n$ while $\mathbb{V}[\hat{\theta}_{\textnormal{BMM}}]=1/(4J(f_{Y|\hat{\bm{\theta}}}(m))^{2})+O(1/J^{2})$. Suppose $J=n$ Dirichlets are sampled to form $\hat{\theta}_{\textnormal{BMM}}$. Then, as $n\to\infty$, $\mathbb{V}[\overline{\theta}]>\mathbb{V}[\hat{\theta}_{\textnormal{BMM}}]$ if $4\sigma^{2}f_{Y|\hat{\bm{\theta}}}^{2}(m)>1$. The larger the tails of the underlying distribution of $\hat{\theta}_{i}$ the more $\sigma^{2}$ grows while $f_{Y|\hat{\bm{\theta}}}(m)$ stays the same, thereby making $\hat{\theta}_{\textnormal{BMM}}$ more attractive. Of course, variance is only half of the picture. Since the bias in using $\hat{\theta}_{\textnormal{BMM}}$ to estimate $m$ is of order $O(1/J)$, the next subsection considers how far apart $m$ and $\theta$ can be.

4.2.2 Concentration bounds for $|m-\theta|$

Since $\hat{\theta}_{\textnormal{BMM}}$ concentrates around $m=\operatorname{median}(Y|\hat{\theta}_{1},\ldots,\hat{\theta}_{n})$, it is important to characterize how $m$ and $\theta$ differ. The results in this section show that when $\mathbb{V}[\hat{\theta}_{i}]=\sigma^{2}$ is small, $m$ and $\theta$ are close and so are $\hat{\theta}_{\textnormal{BMM}}$ and $\overline{\theta}$.

The following mean-median inequality will prove to be extremely valuable. It guarantees that the bias in estimating the median instead of the mean is bounded by the square root of the variance.

Proposition 4.25**.**

If $X$ is a random variable with finite variance $\sigma^{2}$, then

[TABLE]

so the distance between mean and median is at most a standard deviation. Similarly, conditional on random variables $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$, almost surely,

[TABLE]

Proof 4.26**.**

Using Jensen’s inequality and the fact that the median minimizes the $L_{1}$ loss,

[TABLE]

for any $c\in\mathbb{R}$. Taking $c=\mathbb{E}[X]$,

[TABLE]

The conditional result follows analogously by using the conditional Jensen’s inequality.

This inequality can often be strengthened. For example, one could also have taken $c=0$ in the proposition above to obtain $|\operatorname{median}(X)-\mathbb{E}[X]|\leq\mathbb{E}[|X|]$, which sometimes yields stronger results. For unimodal distributions, the upper bound can be tightened to $\sqrt{0.6}\sigma$ (see [Basu and DasGupta, 1997]). If $X$ concentrates exponentially around the mean or median, it is also possible to obtain better results, as the next proposition shows.

Proposition 4.27**.**

Let $X$ be a random variable with mean $\theta$ and median $m$. If there exist $a,b>0$ such that $P[|X-\operatorname{median}(X)|>t]\leq ae^{-\frac{t^{2}}{b}}$ or $P[|X-\mathbb{E}[X]|>t]\leq\frac{a}{2}e^{-\frac{4t^{2}}{b}}$, then

[TABLE]

Proof 4.28**.**

First, note $P[|X-\mathbb{E}[X]|]\leq\frac{a}{2}e^{-\frac{4t^{2}}{b}}$ implies $P[|X-\operatorname{median}(X)|>t]\leq ae^{-\frac{t^{2}}{b}}$. Indeed, consider two cases: (i) $t\geq 2|\mathbb{E}[X]-\operatorname{median}(X)|$, and (ii) $t<2|\mathbb{E}[X]-\operatorname{median}(X)|$. For (i), note

[TABLE]

using the fact that $|X-\operatorname{median}(X)|\leq|X-\mathbb{E}[X]|+|\mathbb{E}[X]-\operatorname{median}(X)|$. For (ii), by the definition of median,

[TABLE]

which implies $2ae^{-\frac{t^{2}}{4b}}\geq 1$, so the bound must hold. This proves exponential mean concentration implies exponential median concentration.

Now, proceeding as in Proposition 4.25,

[TABLE]

On the other hand,

[TABLE]

so, using Proposition 4.25,

[TABLE]

Hence, $|\operatorname{median}(X)-\mathbb{E}[X]|\leq\min(\sqrt{ab},a\sqrt{\pi b}/2)$.

Note that taking expectation on the result of Proposition 4.25 gives the following bound on the bias of the conditional median:

[TABLE]

These results provide useful bounds for the Bayesian median of means estimator.

Proposition 4.29**.**

Suppose $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}[\theta,\sigma^{2}]$, $\mathbf{p}\sim\operatorname{Dir}_{n}(\alpha,\ldots,\alpha)$ and take $Y=\sum_{i=1}^{n}p_{i}\hat{\theta}_{i}$. Then,

[TABLE]

and, conditioning on $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$, almost surely

[TABLE]

In particular, as $n\to\infty$, the unconditional median converges to $\theta$ at least as fast as $O(1/\sqrt{n})$ and the conditional median converges to $\overline{\theta}$ at the same speed.

Proof 4.30**.**

From Proposition 4.25 and (19),

[TABLE]

while Proposition 4.25 and (9) imply

[TABLE]

The results above highlights the impact $\alpha$ has on the difference between $\operatorname{median}(Y|\hat{\bm{\theta}})$ and $\theta$. For instance, if $\alpha=O(n)$, the bound implies a maximum bias of order $O(1/n)$. In fact, any value of $O(n^{\kappa})$, for $\kappa>0$, implies the bias in estimating the median instead of the mean is asymptotically negligible relative to the variance of $\hat{\theta}_{\textnormal{BMM}}$. This, of course, comes at the expense of variance reduction, since larger $\alpha$ mean $\hat{\theta}_{\textnormal{BMM}}$ becomes closer to $\overline{\theta}$. In any case, they are enough to establish $\operatorname{median}(Y|\hat{\bm{\theta}})$ and $\overline{\theta}$ are asymptotically equivalent.

Corollary 4.31**.**

If $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}[\theta,\sigma^{2}]$, $\mathbf{p}\sim\operatorname{Dir}_{n}(\alpha,\ldots,\alpha)$ and $Y=\sum_{i=1}^{n}p_{i}\hat{\theta}_{i}$, then $\operatorname{median}(Y|\hat{\bm{\theta}})-\overline{\theta}\stackrel{{\scriptstyle L_{2}}}{{\longrightarrow}}0$, and $\operatorname{median}(Y|\hat{\bm{\theta}})\stackrel{{\scriptstyle L_{2}}}{{\longrightarrow}}\theta$. In particular, almost surely, $\operatorname{median}(Y|\hat{\bm{\theta}})$ is asymptotically unbiased.

Proof 4.32**.**

Note from Proposition 4.29 that

[TABLE]

so in particular the $L_{2}$ norm of $\operatorname{median}(Y|\hat{\bm{\theta}})$ and $\overline{\theta}$ goes to zero as $n\to\infty$, the more so the smaller $\sigma^{2}$ is. Furthermore, by the triangle inequality,

[TABLE]

In particular, since $L_{2}$ convergence implies $L_{1}$ convergence,

[TABLE]

so $\operatorname{median}(Y|\hat{\bm{\theta}})$ is asymptotically unbiased.

Note Proposition 4.18 bounds the distance between $\hat{\theta}_{\textnormal{BMM}}$ and $m=\operatorname{median}(Y|\hat{\bm{\theta}})$ and Proposition 4.29 bounds the distance between $m$ and $\theta$. Putting them together, it is possible to upper bound the finite-sample (unconditional) bias of the Bayesian median of means.

Proposition 4.33**.**

Let $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}[\theta,\sigma^{2}]$, $\mathbf{p}^{(j)}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}\operatorname{Dir}_{n}(\alpha,\ldots,\alpha)$ and $Y_{j}=\sum_{i=1}^{n}p^{(j)}_{i}\hat{\theta}_{i}$ for $j=1,\ldots,J$, with $J>2$ and unique $m=\operatorname{median}(Y|\hat{\bm{\theta}})$. Consider $t_{0}=4\sqrt{s^{2}_{\hat{\bm{\theta}}}/\alpha}$, and assume there exists $C_{t_{0}}(\hat{\bm{\theta}})>0$ such that $\min_{[m,m+t_{0}]}f(\xi)>C_{t_{0}}(\hat{\bm{\theta}})$. Then

[TABLE]

where $\tilde{C}=\sqrt{\mathbb{E}[1/C_{t_{0}}(\hat{\bm{\theta}})]}$.

Proof 4.34**.**

First, decompose the bias as

[TABLE]

To bound the first term, use Propositions 4.18 and 4.21 to get

[TABLE]

Thus, by the Law of Iterated Expectations and Jensen,

[TABLE]

For the second term on the left-hand side of (25), note

[TABLE]

which yields the bound.

Note to have the bias go to zero it is necessary to take both $J,n\to\infty$. Increasing $J$ concentrates $\hat{\theta}_{\textnormal{BMM}}$ around the conditional median, and increasing $n$ (or $\alpha$) lessens the bias in estimating the median instead of the mean. Indeed, as $\alpha\to\infty$, $\hat{\theta}_{\textnormal{BMM}}$ converges to $\overline{\theta}$ and $\operatorname{median}(Y\mid\hat{\bm{\theta}})=\theta$, but at the expense of possibly higher variance.

The decomposition in (25) also readily gives a bound for the $L_{1}$ error of the procedure.

Corollary 4.35**.**

Consider the setting of Proposition 4.33. Then,

[TABLE]

In particular, if $J=O(n)$, then $\hat{\theta}_{\textnormal{BMM}}$ is consistent.

Proof 4.36**.**

Note that

[TABLE]

Proposition 4.33 gives a bound on the first two terms, while the second can be controlled via Chebyshev: $\mathbb{E}\left[|\overline{\theta}-\theta|\right]\leq\sqrt{\mathbb{E}(\overline{\theta}-\theta)}=\sqrt{\sigma^{2}/n}$. Consistency follows from $L_{1}$ convergence.

The bound on the $L_{1}$ error above, while holding quite generally, is not very useful, particularly in terms of comparing $\hat{\theta}_{\textnormal{BMM}}$ with $\overline{\theta}$, since $\mathbb{E}[|\overline{\theta}-\theta|]$ is needed to upper bound $\mathbb{E}[|\hat{\theta}_{\textnormal{BMM}}-\theta|]$. To get a better understanding of when and how $\hat{\theta}_{\textnormal{BMM}}$ outperforms the sample mean, it is necessary to investigate further asymptotic properties of this estimator using the results of Section 4.1.

4.2.3 Asymptotic approximation

Proposition 4.23 showed the variance of the Bayesian median of means is of order $O(1/J)$, and so is the bias in estimating $\operatorname{median}(Y|\hat{\bm{\theta}})$ with $\hat{\theta}_{\textnormal{BMM}}$. This section develops an asymptotic expansion of $\operatorname{median}(Y|\hat{\bm{\theta}})$ to show the bias in estimating $\theta$ with $\operatorname{median}(Y|\hat{\bm{\theta}})$ is of order $O(1/(n\alpha))$. This implies, when $J=n$ and $\alpha=1$, the variance of $\hat{\theta}_{\textnormal{BMM}}$ is of order $O(1/n)$ while the squared bias is only $O(1/n^{2})$. The asymptotic approximation also gives a deterministic, approximate algorithm for $\hat{\theta}_{\textnormal{BMM}}$.

Proposition 4.37**.**

Let $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}[\theta,\sigma^{2}]$ with $\mathbb{E}[|\hat{\theta}_{i}|^{3}]<\infty$, and denote by $\overline{\theta}$ the sample mean. If $Y=\sum_{i=1}^{n}p_{i}\hat{\theta}_{i}$ for $p\sim\operatorname{Dir}(\alpha,\ldots,\alpha)$, then almost surely with respect to the sampled $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$,

[TABLE]

Proof 4.38**.**

Recall from Proposition 4.8 that, conditional on $\hat{\bm{\theta}}=(\hat{\theta}_{1},\ldots,\hat{\theta}_{n})$,

[TABLE]

where $\tilde{Y}$ is the standardized version of $Y$. Since $\mathbb{E}[\hat{\theta}_{i}^{3}]<\infty$, the Edgeworth expansion of [Weng, 1989] implies that, almost surely with respect to the empirical distribution of $\hat{\theta}_{i}$,

[TABLE]

where $\Phi(\tilde{y})$ and $\varphi(\tilde{y})$ denote the cumulative distribution and probability density functions, respectively, of a standardized Normal random variable, and \mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{skew}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{skew}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{skew}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{skew}}}(\hat{\bm{\theta}}) is the sample skewness of $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$, that is,

[TABLE]

The median of a standardized Normal is zero, so, from Proposition 4.29, it suffices to consider the median of $\tilde{Y}$ to be $\tilde{m}=O(1/\sqrt{n\alpha})$. In this case, a Taylor expansion yields

[TABLE]

and so plugging in $\tilde{y}=\tilde{m}$ in (27) gives

[TABLE]

so

[TABLE]

Since $\operatorname{median}(Y|\hat{\bm{\theta}})=m=\overline{\theta}+\sqrt{\frac{1}{n\alpha+1}s^{2}_{\hat{\bm{\theta}}}}\tilde{m}$, this implies

[TABLE]

so the conditional median is, to first-order, a variance and skewness correction applied to the sample mean.

Recall that, as $J\to\infty$, $\operatorname{median}(Y|\hat{\bm{\theta}})$ becomes well-approximated by $\hat{\theta}_{\textnormal{BMM}}$, and increasing the value of $J$ is generally easy since it only involves sampling Dirichlets. In this sense, the previous proposition illuminates important aspects of the Bayesian median of means estimator. Several key remarks are collected below.

Relationsip between $\hat{\theta}_{\textnormal{BMM}}$ and $\overline{\theta}$. Formula (26) synthesizes many previous results. For example, it makes it clear that as $\alpha\to\infty$, $\hat{\theta}_{\textnormal{BMM}}$ is converging to $\overline{\theta}$ and so, in particular, $\hat{\theta}_{\textnormal{BMM}}$ is asymptotically unbiased. To first-order, the two estimators differ when the sample variance or the sample skewness are very large, in which case $\hat{\theta}_{\textnormal{BMM}}$ applies a correction that reduces variance at the expense of bias. In cases where $\overline{\theta}$ is known to be optimal, for instance when $\hat{\theta}_{i}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}N(\theta,\sigma^{2})$, $\hat{\theta}_{\textnormal{BMM}}$ and $\overline{\theta}$ are virtually the same. Also, while $\overline{\theta}$ performs worse the higher $\sigma^{2}$ is, $\hat{\theta}_{\textnormal{BMM}}$ does not have its asymptotic variance affected by $\sigma^{2}$, but its bias roughly depends on $\sigma^{2}\cdot\operatorname{skew}(\hat{\bm{\theta}})$. Thus, the ideal setting for $\hat{\theta}_{\textnormal{BMM}}$ relative to $\overline{\theta}$ is one with high variance, low skewness — precisely those of Example 3.3.

How $\alpha$ affects convergence. The values of $n$ and $\alpha$ are directly related. If $\alpha$ is much smaller than $1$, the bias of the estimator can be relatively large, even asymptotically. For instance, if $\alpha=1/n$ the Bayesian median of means ceases to be asymptotically unbiased, as there is a constant correction to the sample mean, possibly incurring high variance. Values of $\alpha$ much larger than 1 suggest more data is available than is actually the case, and thus the estimator relies more on the sample mean. Proposition 4.8 also holds when $\alpha\to\infty$, since a $\text{Gamma}(\alpha)$ can be thought of as the convolution of $\alpha$ Exponential random variables, in which case the Central Limit Theorem yields $\sqrt{\alpha}(Y-\overline{\theta})\mid\hat{\bm{\theta}}\Rightarrow N(0,s^{2}_{\hat{\bm{\theta}}}/n)$, so Proposition 4.37 still carries through for $\alpha\to\infty$ with fixed $n$.

Example 4.41 compares the behavior of $\hat{\theta}_{\textnormal{BMM}}$ and $\overline{\theta}$ for varying values of $n$ and $\alpha$ for a Skewnormal simulation. The larger $n$ or $\alpha$, the more the estimators look alike.

Cornish-Fisher expansion. The proof of Proposition 4.37 can be thought of as the development of a Cornish-Fisher expansion for the median of $Y|\hat{\bm{\theta}}$, which provides a full asymptotic expansion of the quantiles of $F_{Y|\hat{\bm{\theta}}}(y)$ in terms of its cumulants. For example, using (15), the skewness of $Y|\hat{\bm{\theta}}$ can be written

[TABLE]

which appears as part of the coefficient of the first factor in (27). Section 4.1 established the density of $Y|\hat{\bm{\theta}}$ and a recursive formula for higher moments, which is what is needed to develop further terms and establish better approximations. In particular, a formal Cornish-Fisher expansion implies the remainder in (26) is actually $O(\frac{1}{n^{2}\alpha^{2}})$.

Bias of $\hat{\theta}_{\textnormal{BMM}}$ is asymptotically negligible. The motivation for the Bayesian median of means is to improve on the mean squared error of the sample mean by trading off some variance for bias. For this approach to work, it is paramount that the squared bias introduced is not larger than the reduction in variance. The following corollary provides an asymptotic assurance.

Corollary 4.39**.**

Let $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}[\theta,\sigma^{2}]$, with $\mathbb{E}[|\hat{\theta}_{i}|^{3}]<\infty$, and consider $Y_{j}=\sum_{i=1}^{n}p^{(j)}_{i}\hat{\theta}_{i}$, where $p^{(j)}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}\operatorname{Dir}(\alpha,\ldots,\alpha)$, $j=1,\ldots,J$. Take and $J=n$ and $\alpha=1$. If $\operatorname{median}(Y|\hat{\bm{\theta}})$ is unique, then almost surely

[TABLE]

Proof 4.40**.**

If $\alpha=1$ and $J=O(n)$, then Proposition 4.37 implies the bias between $\operatorname{median}(Y|\hat{\bm{\theta}})$ and $\theta$ is of order $O\left(\frac{1}{n}\right)$, while Proposition 4.23 imply the bias between $\hat{\theta}_{\textnormal{BMM}}$ and $\operatorname{median}(Y|\hat{\bm{\theta}})$ is of order $O\left(\frac{1}{n}\right)$, so the overall bias of $\hat{\theta}_{\textnormal{BMM}}$ is $O(\frac{1}{n})$. From Proposition 4.23, the variance is of order $O\left(\frac{1}{n}\right)$.

The crucial implication of the above corollary is that, in terms of mean squared error, the squared bias incurred by the Bayesian median of means is of order $O(\frac{1}{n^{2}})$, and so it is negligible with respect to the variance, which is of order $O(\frac{1}{n})$. Put another way, for $n$ large enough, the expected decrease in variance of the Bayesian median of means is sure to make its mean squared error smaller than that of the sample mean. Determining how large $n$ should be, however, for the asymptotic regime to be sufficiently accurate depends on the underlying distribution of $\hat{\theta}_{i}$; in particular, note the MLE generally exhibits a similar behavior of trading-off some bias for variance, but it does so by requiring the distribution to be specified.

Deterministic approximation to $\hat{\theta}_{\textnormal{BMM}}$. Proposition 4.37 also suggests a new, robust mean estimator when $\mathbb{E}[|\hat{\theta}_{i}|^{3}]<\infty$: the approximate Bayesian median of means, defined as

[TABLE]

Note $\hat{\theta}_{\textnormal{aBMM}}$ is deterministic and simpler to compute than the original $\hat{\theta}_{\textnormal{BMM}}$, since it doesn’t require the sampling of auxiliary Dirichlet random variables. It only has one hyperparameter, namely $\alpha$, and Section 4.3 considers how to set it. Example 4.41 shows that $\hat{\theta}_{\textnormal{aBMM}}$ is accurate in replicating the performance of $\hat{\theta}_{\textnormal{BMM}}$ for a Skewnormal example; further comparisons are provided in Section 5.1. Better approximations are possible by using additional terms of the Cornish-Fisher expansion, but at the expense of more computation and assumptions on the distribution of $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$.

Example 4.41** (Approximate BMM).**

Let $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$ be a sample of $n=1000$ iid samples from a Skewnormal distribution, with location $\xi=0$, scale $\omega=20$ and shape $\beta=10$. The sample mean is $\overline{\theta}\approx 15.625$ and the true mean is $\theta\approx 15.878$ (with variance $147.87$ and skewness $0.956$).

Figure 6 illustrates the result of applying both the approximate and the exact Bayesian median of means procedures to such data, with $J=n=1000$ and varying values of $\alpha$. The violin plot in the figure, drawn from $100$ simulations, shows the distribution of $\hat{\theta}_{\textnormal{BMM}}$ is approximately a scaled Normal, in accordance with Proposition 4.16, and that $\hat{\theta}_{\textnormal{aBMM}}$ and $\hat{\theta}_{\textnormal{BMM}}$ agree with each other the bigger $\alpha$ is.

Figure 6 considers the same setting, but with fixed $\alpha=1$ and varying values of $n$. Note the sampled $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$ change with $n$, which is why there are variations in each case. Overall, it is clear that $\hat{\theta}_{\textnormal{aBMM}}$ provides a decent approximation to $\hat{\theta}_{\textnormal{BMM}}$, and they are generally similar to $\overline{\theta}$, especially as $n$ increases, in which case all estimators approach the value of $\theta$. For these simulations, $J$ is kept fixed at $1000$.

4.3 Choosing $\alpha$

The Bayesian median of means algorithm is fully specified, except for the choice of hyperparameter $\alpha$. This section discussesr how to set $\alpha$, from picking it independently of the data to more data-driven choices. For its computational simplicity and statistical properties, $\alpha=1$ is taken to be the recommended value.

4.3.1 Setting $\alpha=1$

Recall large values of $\alpha$ approximate $\hat{\theta}_{\textnormal{BMM}}$ to $\overline{\theta}$, so they induce less bias but more variance. Ideally, one would set $\alpha$ large enough so as to have minimum bias while still keeping the variance in control. From this regard, recall \hat{\theta}_{\textnormal{BMM}}=\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{median}}}(Y_{1},\ldots,Y_{J}), where

[TABLE]

The second term above is not affected by $\alpha$, so one could pick $\alpha$ to have the first term of the same asymptotic order of the second. This amounts to setting $\alpha=1-2/n$. For sufficiently large $n$, it makes sense to simply set $\alpha=1$.

Computationally, the choice $\alpha=1$ is also advantageous, since in this case sampling $(p_{1},\ldots,p_{n})\sim\operatorname{Dir}_{n}(1,\ldots,1)$ can be done by sampling $n-1$ Uniform random variables, $U_{1},\ldots,U_{n-1}$, and ordering them. Then, set $U_{0}=0$, $U_{n}=1$ and take $p_{i}=U_{(i)}-U_{(i-1)}$ for $i=1,\ldots,n$. For extremely large values of $n$, sorting is more expensive than sampling Gammas, the usual way of obtaining a Dirichlet draw; still, for the case $\alpha=1$, Exponentials can be sampled instead of Gammas, which are much faster.

Other canonical, data-independent choices for $\alpha$ are $\alpha=4$, so the first two terms of the Edgeworth expansion (27) match that of the regular bootstrap, and $\alpha=0.5$, which amounts to the half-sampled bootstrap scheme considered in [Friedman and Hall, 2007].

4.3.2 Picking $\alpha$ via prior information

One of the advantages of the Bayesian median of means is that prior information can be easily incorporated. If the user is confident the underlying distribution is close to Normal, a higher $\alpha$ should be set; if they expect a distribution with high or infinite variance, such as the importance sampling estimator in Section 1, then lower values of $\alpha$ are better. In particular, $\alpha$ can be attributed its own prior, say $\alpha\sim\text{Beta}(\beta_{1},\beta_{2})$, and have its parameters determined via an Empirical Bayes approach. More generally, the prior might depend on sample quantities such as standard deviation or skewness.

Prior information can also help set $\alpha$ through asymptotic considerations. Recall the bias of $\hat{\theta}_{\textnormal{BMM}}$ is of order $O(1/J)+O(1/n\alpha)$, so, for instance, $\alpha=1/\sqrt{n}$ and $J=n\alpha$ imply the bias squared is of the same order as the variance, which might be desirable if it is known beforehand that variance is a bigger concern than bias.

4.3.3 Picking $\alpha$ adaptively

It is also possible to pick $\alpha$ depending on the sampled values $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$ by using cross-validation. That is, for a given choice of $\alpha$, split $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$ into $k$ folds, and use $\hat{\theta}_{\textnormal{BMM}}$ with $k-1$ folds worth of data to estimate, say, the sample mean of the unseen fold as a proxy for $\theta$. Average the errors over the folds to obtain an error estimate, and pick the $\alpha$ that yields lowest error. Similarly, if $\alpha$ is set much smaller than $1$, most samples $\mathbf{p}^{(1)},\ldots,\mathbf{p}^{(J)}$ will contain coordinates very close to zero. Setting them to zero amounts to not using some of the $\hat{\theta}_{i}$ in creating $Y_{j}$, so these $\hat{\theta}_{i}$ can be thought of as out-of-bag samples, and the $\alpha$ that best predicts the $\hat{\theta}_{i}$ not used is selected. Note, however, that the number of folds now becomes another hyperparameter to be determined. Furthermore, if the sample is highly skewed or with large variance then cross-validation is expected to fail, since there might be severe mismatches between the folds.

5 Empirical Results

This section considers the empirical behavior of the Bayesian median of means in a variety of settings. The results will generally be compared against the sample mean, a standard nonparametric location estimator, using mean squared error loss. By default, the Bayesian median of means will use $\alpha=1$ and $J=n$. The full procedure is given in Algorithm 1 below. Note this is readily parallelizable. All the code to generate the figures and examples can be found at https://github.com/paulo-o/bmm.

The examples below range from low to high-variance distributions, coming from real and simulated data. Section 5.1 compares the behavior of the Bayesian median of means, $\hat{\theta}_{\textnormal{BMM}}$, against four other candidates: the sample mean, $\overline{\theta}$; the approximate Bayesian median of means, $\hat{\theta}_{\textnormal{aBMM}}$; the classical median of means, $\hat{\theta}_{\textnormal{MM}}$; and the sample median, \mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{median}}}(\hat{\bm{\theta}}). Section 5.2 considers the issue of developing confidence intervals for $\hat{\theta}_{\textnormal{BMM}}$ in an efficient manner, while Section 5.3 applies $\hat{\theta}_{\textnormal{BMM}}$ to procedures such as importance sampling, cross-validation and bagging. Note in some of these cases the estimators $\hat{\theta}_{i}$ are no longer independent, as was assumed throughout the paper, but the conditional results still hold, and it is interesting to see how the Bayesian median of means fares in this context.

5.1 Comparison with sample mean and median of means

Recall both $\hat{\theta}_{\textnormal{BMM}}$ and $\hat{\theta}_{\textnormal{MM}}$ can be thought of as interpolating between the sample mean and the sample median. The simulations below explore how $\hat{\theta}_{\textnormal{BMM}}$, $\hat{\theta}_{\textnormal{aBMM}}$ and $\hat{\theta}_{\textnormal{MM}}$ fare in settings that favor these extremes. Since the the Bayesian median of mean has its hyperparameter set at $\alpha=1$, a fair comparison would have the median of means with as few groups $g$ as possible. Since $g=1$ and $g=2$ lead to the sample mean, $g=3$ is used in the examples below.

Example 5.42** (Skewnormal).**

Suppose $\hat{\theta}_{i}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}\text{SkewN}(\xi,\sigma,\lambda)$, with location parameter $\xi=0$, scale $\sigma=1000$ and shape $\lambda=0$ (so this is a Normal distribution). The boxplot of $\overline{\theta}$, $\hat{\theta}_{\textnormal{aBMM}}$, $\hat{\theta}_{\textnormal{BMM}}$, $\hat{\theta}_{\textnormal{MM}}$ and \mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{median}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{median}}}(\hat{\theta}_{1},\ldots,\hat{\theta}_{n}), evaluated over $1000$ simulations, is shown on the left side of Figure 7. The numbers in the figure correspond to the mean squared error of each estimator. As expected, all estimators are unbiased, but the median has greater variance, contributing to worse performance. The median of means, being closer to the median than the Bayesian median of means, shows slightly degraded results. Both the Bayesian median of means and its approximation, as well as the sample mean, exhibit similar performance.

On the other hand, consider increasing the skewness of the distribution to $\lambda=40$. Now, the underlying distribution of the $\hat{\theta}_{i}$ is no longer symmetric, so the median becomes severely biased, as shown on the right side of Figure 7. Note, however, that $\hat{\theta}_{\textnormal{BMM}}$ and $\hat{\theta}_{\textnormal{aBMM}}$ incur minimal bias and still have a similar performance to $\overline{\theta}$.

Example 5.43** (Pareto).**

Let $\hat{\theta}_{i}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}\text{Pareto}(\xi,\sigma,\lambda)$ with location parameter $\xi=0$, scale $\sigma=1000$ and two possible shapes $\lambda=4$ (lower skewness) and $\lambda=2.5$ (higher skewness). In both cases, displayed in Figure 8, there is enough bias in using the median that it performs far worse than any other estimator. The sample mean works well when the skewness is small, but it is still comparable to both $\hat{\theta}_{\textnormal{BMM}}$ and $\hat{\theta}_{\textnormal{aBMM}}$. Once the skewness and variance increase, the sample mean performs worse than either, particularly because it gives weights to extreme sample points. The median of means exhibits slightly worse performance than both the exact and approximate Bayesian median of means.

Example 5.44** (Discrete).**

Assume now a distribution that is particularly favorable to the sample median relative to the sample mean: suppose $\hat{\theta}_{i}$ is distributed as

[TABLE]

Note $\mathbb{E}[\hat{\theta}_{i}]=0$ and the distribution is symmetric, so the median is unbiased. On the other hand, the variance is $\mathbb{V}[\hat{\theta}_{i}]=n^{4-p}\sigma^{2}$. For almost all samples drawn $\hat{\theta}_{i}=0$, which is the true value of the parameter. Rarely, however, $\hat{\theta}_{i}$ attains a large value that is capable of throwing the sample mean off, while the sample median remains immune.

Figure 9 shows the result of the different estimators when $n=1000$, $\sigma=30$ and $p=1$ (high probability of extremes) or $p=1.5$ (lower). Note in both cases the sample median attains zero mean squared error, while the sample mean behaves poorly. When $p$ is bigger, the extreme values seldom survive the median operation, so the median of means achieves good performance, and less so the Bayesian median of means, as it doesn’t give weight zero to any single $\hat{\theta}_{i}$. When the probability of extreme events increase, the median of means becomes relatively more susceptible to extreme measurements, and displays worse results. Also, because the underlying distribution is symmetric, \mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{skew}}}{\accentset{\textstyle\text{\smash{\raisebox{-5.59721pt}{\widehatsym}}}}{\operatorname{skew}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-3.91806pt}{\widehatsym}}}}{\operatorname{skew}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-2.7986pt}{\widehatsym}}}}{\operatorname{skew}}}(\hat{\bm{\theta}})\approx 0, so $\hat{\theta}_{\textnormal{aBMM}}$ is virtually the same as the sample mean. In this case, the approximation fails, as the neglected higher-order terms in the expansion become consequential.

The behavior of the estimators in this last example can be understood as follows. While the median of means is much better protected against extremes because of the median operation, it is less efficient in the sense that it underutilizes the data, as the final estimator has only used $n/g$ datapoints, where $g$ is the number of groups. The Bayesian median of means is somewhere in between the two. Indeed, write

[TABLE]

In each case, $\mathbb{E}[\overline{\theta}]=\mathbb{E}[Y_{1}]=\mathbb{E}[\overline{\theta}_{1}]$, but the variances are very different:

[TABLE]

It is clear $\mathbb{V}[Y_{1}]\geq\mathbb{V}[\overline{\theta}]$, while $\mathbb{V}[\overline{\theta}_{1}]\geq\mathbb{V}[Y_{1}]$ if $n\alpha\geq(\alpha+1)n/g$. With $\alpha=1$, this is always the case (since $g\geq 2$), and the inequality becomes more pronounced as $\alpha\to\infty$, that is, as $\hat{\theta}_{\textnormal{BMM}}\to\overline{\theta}$.

However, while the terms in $\hat{\theta}_{\textnormal{BMM}}$ have less variance than those of $\hat{\theta}_{\textnormal{MM}}$, the Bayesian median of means is less immune to extremes. Indeed, the blocks $\overline{\theta}_{i}$ are completely independent of each other, while the $Y_{j}$ are dependent, as they all rely on the same values of $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}$. Still, this dependence is weaker than that of the sample mean, where all terms are equal to $\overline{\theta}$. In particular, the sample mean has an asymptotic breakdown point of zero, so it suffices to have one outlier to severely alter it, while the median of means has an asymptotic breakdown point equal or bigger to $g/2n$, since at least $g/2$ outliers are necessary, but not sufficient, to damage it. The Bayesian median of means sits somewhere in-between as $\alpha$ varies, though the notion of breakdown point is harder to quantify in this case.

Consider now how the estimators behave as $n$ changes. Recall that as $n\to\infty$, both $\hat{\theta}_{\textnormal{BMM}}$ and $\overline{\theta}$ should become indistinguishable. In the examples below, $n=0,10,20,30,\ldots,2000$, and $100$ simulations are used for each $n$ to obtain the mean squared error. As before, the median of means estimators uses $g=3$ groups. All the figures below have the $y$-axis in $\log_{10}$ scale.

Example 5.45** (Exponential-$t$ distribution).**

Let $\hat{\theta}_{1},\ldots,\hat{\theta}_{n}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}\text{Expo}(\lambda)+t_{2.5}(0,\sigma)$, where $\sigma=\frac{1}{\lambda}=30$, so the $\hat{\theta}_{i}$ are not symmetric and have a distribution with relatively heavy tails. Figure 10 shows the decaying MSE for each estimator, except the median. Note the approximate and exact Bayesian median of means have virtually the same performance, while the sample mean exhibits consistently bigger MSE. The median of means displays the worst performance in this scenario because it relies too much on the median, which has MSE close to $17$ for large $n$ due to the bias.

Example 5.46** (Lognormal).**

Suppose $\hat{\theta}_{i}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}\text{Lognormal}(\theta,\sigma^{2})$, with $\mu=4$ and $\sigma=1$ (so the shape parameter is 1, the location is 0 and the scale is $e^{4}$). Recall the Lognormal distribution is in the exponential family, so the sample mean is the MLE. Figure 11 shows the median of means has consistently worse performance. The Bayesian median of means and the sample mean are very close to each other, although the Bayesian median of means is generally better, in particular for small $n$. The approximate Bayesian median of means is not shown because it is virtually identical to its exact counterpart. The median is not shown because its MSE revolves around $1250$ and does not diminish with $n$ (the bias squared in using the median is $(e^{4.5}-e^{4})^{2}\approx 1254$).

Example 5.47** (Maximal bias).**

For a given variance level, Proposition 4.25 ensures that the maximum distance between mean and median is the standard deviation. This can be achieved by sampling $\hat{\theta}_{i}$ as

[TABLE]

for sufficiently small $\varepsilon$, so $\theta\approx 0$ while $\operatorname{median}(\hat{\theta}_{1})=\sigma$. Note this represents the worst case for any median-based procedure, since the bias is maximal. Figure 12 shows the result of applying the sample mean, the Bayesian median of means and the median of means to this distribution. For small $n$ the sample mean dominates by a small amount, but that gap quickly disappears as $n$ increases.

5.2 Confidence intervals

Constructing confidence intervals for $\hat{\theta}_{\textnormal{BMM}}$ or $\hat{\theta}_{\textnormal{aBMM}}$ is not so straightforward since the samples $Y_{1},\ldots,Y_{J}$ are not independent, and confidence intervals for the median are usually based on asymptotic considerations. An alternative is to use the bootstrap. That is, sample $\hat{\theta}_{1}^{(b)},\ldots,\hat{\theta}_{n}^{(b)}$ with replacement from $\{\hat{\theta}_{i}\}_{i=1}^{n}$, $b=1,\ldots,B$ times. For each bootstrap sample $b$, find the Bayesian median of means estimator, $\hat{\theta}_{\textnormal{BMM}}^{(b)}$, and use the empirical distribution to generate confidence intervals for $\theta$ under the assumption that $P[\hat{\theta}_{\textnormal{BMM}}-\theta\leq t]\approx P[\hat{\theta}_{\textnormal{BMM}}^{(b)}-\hat{\theta}_{\textnormal{BMM}}\leq t]$. The examples below show that this procedure yields coverage near the prescribed levels.

Unfortunately, it might be hard or unfeasible to obtain $B$ bootstrap samples for the Bayesian median of means since it involves $O(Bn^{2})$ operations. A computational shortcut is to fix the Dirichlet draws $\mathbf{p}^{(1)},\ldots,\mathbf{p}^{(n)}$ for all $b=1,\ldots,B$ samples, reducing the complexity to $O(Bn)$, which is considerably faster and similar to bootstrapping the sample mean. One would not expect this to significantly change the estimates, since the $\hat{\theta}_{i}^{(b)}$ are being sampled independently with replacement, but it does add correlation across the samples. In the examples below, the effect of fixing the Dirichlet draws is negligible, while the computational speedup is considerable. Also, both percentile and BCa intervals were analyzed, and the difference was again minor, likely due to the stabilizing effect of the Dirichlet averages. See Algorithm 2 for the full description of the percentile bootstrap used.

Example 5.48** (CIs for Exponential, Pareto and Normal distributions).**

Figure 13 shows the confidence intervals at the $95\%$ level for $\hat{\theta}_{\textnormal{BMM}}$ when (i) $\hat{\theta}_{i}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}\text{Expo}(1/3)+5$; (ii) $\hat{\theta}_{i}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}\text{Pareto}(4,10)$; and (iii) $\hat{\theta}_{i}\operatorname{\stackrel{{\scriptstyle\text{iid}}}{{\sim}}}N(0,1)$. The $x$-axis refers to $1000$ different draws of $\{\hat{\theta}_{i}\}_{i=1}^{n}$, the black dots represent $\hat{\theta}_{\textnormal{BMM}}$ and the grey lines are the intervals. The dotted line in red is the actual value of $\theta$, and the red crosses are instances where the interval doesn’t cover. Note the coverage over these $1000$ draws is close to the nominal level of $95\%$.

5.3 Applications

Many statistical procedures are based on the idea of aggregating multiple unbiased estimators. This section considers the performance of the Bayesian median of means in importance sampling, bagging and cross-validation. The takeaway is that as long as the estimators being aggregated are low variance and low skewness, the Bayesian median of means performs well and comparable to the sample mean. However, once the estimators exhibit high variance then the sample mean will frequently underperform, and using the Bayesian median of means yields significant gains in mean squared error.

5.3.1 Importance sampling

A class of problems that depend on the aggregation of unbiased estimates with potentially huge variance is importance sampling, as previously illustrated in Example 1.1. As another application of importance sampling, consider estimating the number of Fibonacci permutations.

A Fibonacci permutation refers to a permutation of $m$ objects such that element $i$ is restricted to positions $i-1$, $i$ or $i+1$ ($1$ can be placed in position $1$ or $2$, and $m$ can placed in position $m-1$ or $m$). For example, if $m=7$, then

[TABLE]

are examples of Fibonacci permutations, while $[3,1,2,4,5,7,6]$, $[1,2,5,4,3,6,7]$ and $[2,1,4,3,6,7,5]$ are not. How to efficiently estimate the number of Fibonacci permutations for arbitrary $m$?

On the one hand, it is not hard to see that the number of such Fibonacci permutations is given by $\text{Fib}(m+1)$, where $\text{Fib}(m)$ denotes the $m$-th element in the Fibonacci sequence, $0,1,1,2,3,5,8,\ldots$. Indeed, let $P_{m}$ be the number of Fibonacci permutations of $m$ elements. Object $1$ can be placed in positions $1$ or $2$; if in position $1$, then there are $P_{m-1}$ possibilities left; if in position $2$, then position $1$ must necessarily have object $2$, and thus there are $P_{m-2}$ possibilities left. Hence, $P_{m}=P_{m-1}+P_{m-2}$ and $P_{1}=1$ give the Fibonacci recursion, so $P_{m}=\text{Fib}(m+1)$.

On the other hand, importance sampling can also be used for this task. The fact that the solution is known means it is possible to faithfully investigate whether using the Bayesian median of means yields any improvement to the method. Also, there are many other similar problems that do not admit closed-form solutions (see [Diaconis et al., 2001]), so this serves as a benchmark.

To use importance sampling, let $Z_{m}$ denote the number of Fibonacci permutations of $m$ objects, which is the target of the exercise, and take $p(x)=\frac{1}{Z_{m}}$ to be the Uniform distribution, which is hard to sample from. Consider an alternative distribution that is easier to simulate: (i) the first object’s position is to be sampled uniformly at random from $\left\{1,2\right\}$; (ii) if the first object went to position $1$, then the second object’s position is to sampled uniformly at random from $\left\{2,3\right\}$; if the first object went to position $2$, then the second object must necessarily go to position $1$; (iii-) keep picking among available alternatives for each object’s position until all objects are placed. This defines a sequence of conditional distributions for $X_{i}=(X_{i}^{(1)},X_{i}^{(2)},\ldots,X_{i}^{(m)})$:

[TABLE]

Denote this distribution by

[TABLE]

and note that if $\eta_{i,j}$ is the number of neighboring positions available to place $X_{i}^{(j)}$, (so $\eta_{i,1}=2$ for all $i$; $\eta_{i,2}=2$ if $X_{i}^{(1)}=1$ and $\eta_{i,2}=1$ if $X_{i}^{(1)}=2$, etc.), then

[TABLE]

that is, $q(X_{i})$ is the reciprocal of the product of the available neighbors for each position when constructed sequentially.

Now, since it is easy to sample from $q$, $Z_{m}$ can be estimated using the usual importance sampling estimator:

[TABLE]

where $\eta_{i}$ refers to the neighbors available in constructing $X_{i}$. This can be done using Algorithm 3 below.

Figure 14 shows histograms for $\hat{\theta}_{i}$, for $i=1,\ldots,1000$, when $m=20,50,80,100$. Note the problem becomes harder as $m$ increases, in particular due to the presence of extremely large but unlikely values of $\hat{\theta}_{i}$, corresponding to rare permutations (so $q(X_{i})$ is very small, and thus $1/q(X_{i})$ is very big). These large weights are essential for unbiasedness, but add substantial variance to the estimates.

Thus, importance sampling requires generating Fibonacci permutations $X_{i}$ from $q$ and calculating

[TABLE]

where $\eta_{i}^{(j)}$ refer to the number of possible positions object $j$ could take when the previous objects $1,\ldots,j-1$ have been placed at positions $X_{i}^{(1)},\ldots,X_{i}^{(j-1)}$. The usual importance sampling estimate is obtained by aggregating these estimates via the sample mean:

[TABLE]

As an alternative, one can consider aggregating $\hat{\theta}_{i}$ via the Bayesian median of means.

Figure 15 has the results of using $\hat{\theta}_{\textnormal{IS}}$, $\hat{\theta}_{\textnormal{BMM}}$ and $\hat{\theta}_{\textnormal{aBMM}}$ when $m=20,50,80,100$ and $n=1000$. It becomes quite clear that $\hat{\theta}_{\textnormal{BMM}}$, and even more so $\hat{\theta}_{\textnormal{aBMM}}$, improve on the usual importance sampling estimates by severely biasing the large weights towards the mean (as illustrated in Figure 3), resulting in larger gains in terms of MSE the larger the $m$.

Figure 16 shows the relative MSE, measured as $\mathbb{E}[(\hat{\theta}-\theta)^{2}/\theta^{2}]$, for $\overline{\theta}$, $\hat{\theta}_{\textnormal{BMM}}$ and $\hat{\theta}_{\textnormal{aBMM}}$. The top picture displays the case in which $n=100$ samples are used, and the lower picture has $n=1000$. Note that higher $m$ makes $\hat{\theta}_{\textnormal{BMM}}$ and $\hat{\theta}_{\textnormal{aBMM}}$ better since the variance component of MSE overwhelms the bias. For small $m$, all estimators perform comparably. Also, higher $n$ means all estimators are closer to each other, and the differences in relative MSE are lower.

5.3.2 Cross-validation

The Bayesian median of means can be used to improve many statistical methods that involve resampling the data. For instance, consider cross-validation, which is used to estimate the test error of a procedure, say linear regression. One divides the data into folds, and try to predict each fold using the other folds as training data. The estimated test error in each fold is aggregated via the sample mean, but consider using the Bayesian median of means instead. The gains should be especially evident when the estimates from different folds vary significantly.

Below, both methods are compared by estimating the test error in a variety of datasets. In each case, $20\%$ of the data is set aside as a validation set, and the remaining $80\%$ is used to generate the usual cross-validation (CV) estimates as well as the Bayesian median of means aggregation (BMM). Then, the full training data is used to train the statistical procedure, and the validation set provides an accurate estimate of test error (error). For each dataset, this is repeated for 100 different seeds, and the root mean squared error (RMSE) between CV and error, and BMM and error are computed.

The datasets were chosen from [James et al., 2013] and include the following:

1. (i)

Advertising: given advertising budgets for TV, radio and newspaper, predict sales of a product ($n=200$ data points, $p=4$ features). 2. (ii)

Auto: given several features about a car, such as horsepower, weight and acceleration, predict its overall displacement ($n=392$, $p=9$). 3. (iii)

College: using features from different colleges, such as whether they are private, number of accepted or enrolled students and number of graduate students, predict its graduation rate ($n=777$, $p=19$). 4. (iv)

Credit: knowing a person’s income, age, education and other financial data, predict their credit card balance ($n=400$, $p=12$). 5. (v)

Heart: given features from people with chest pain, such as their age, sex and fitness measurements, predict whether they have heart disease or not ($n=297$, $p=17$).

The first four datasets above involve regression tasks, while the last one is a classification problem. For the first four, the underlying procedures used were linear regression, random forest (with $20$ different trees) and $k$-nearest neighbors (with $5$ nearest neighbors). All features available were used, and dummies were created when necessary. The loss functions used was mean squared error for regression, and misclassification error for classification. The Bayesian median of means was used with $J=1000$ Dirichlet draws.

Generally, Tables 1-5 show that the Bayesian median of means adds a small amount of stability to the cross-validation estimates, resulting in lower mean squared error. While the improvement is consistent across different datasets, it is often small. This is likely the consequence of estimates from different folds begin similar, so the Bayesian median of means cannot leverage its gains in variance reduction. In fact, in some cases, the usual cross-validation exhibit smaller standard error, such as in Table 3, but the only case in which cross-validation exhibits smaller MSE is in Table 4 for linear regression and $k$-nearest neighbors. Some care must be taken in interpreting these results, as real data represents only one realization from the probability distribution characterizing the data-generating mechanism. Still, the results support the conclusion that the Bayesian median of means is a competitive way to perform aggregation in data resampling schemes.

5.3.3 Ensembling

The Bayesian median of means can also be used to improve many statistical algorithms that rely on ensemble learning, such as bagging. In this case, a particular algorithm, say a regression tree, is trained on different subsets of the data to produce distinct final estimates, with the purpose of achieving some variance reduction. This section considers using the Bayesian median of means to boost these gains.

The advantages in using the Bayesian median of means are greater the more varying the estimates are. Thus, highly non-convex methods or high-variance data provide a setting in which aggregation via the Bayesian median of means can outperform classical bagging. Note, however, that the data here are not independent anymore, since there is a large overlap between training folds. Still, all the conditional guarantees from Section 4.2 hold.

To investigate how these methods fare in real datasets, consider the performance of using both the sample mean and the Bayesian median of means in conjunction with bagging. Six randomly chosen datasets from the UCI data repository ([Dua and Graff, 2017]) were picked to span different domains of interest:

1. (i)

bike: with features such as weather and number of registered users at the moment, the target is the number of bikes rented at a given time ($n=17379$ data points, $p=15$ features); 2. (ii)

demand: for given a collection of business metrics from a large logistics company, such as orders of different types, predict the total daily orders ($n=60$, $p=11$); 3. (iii)

fires: predict the burned area of forest fires in Portugal using meteorological and spatial data ($n=517$, $p=10$); 4. (iv)

GPUs: predict the running time of for the multiplication of two $2048\times 2048$ matrices using a GPU OpenCL SGEMM kernel with features regarding the computational environment and past runs ($n=241600$, $p=16$); 5. (v)

news: given features regarding an online news webpage such as number of words, links, date of publication and subject area, predict how popular the article will be ($n=39644$, $p=59$); 6. (vi)

superconductors: given features about superconductors, predict the critical temperature ($n=21263$, $p=80$).

Minimal modifications were made to the datasets, such as converting categorical predictors to dummy variables.

For each dataset, $15\%$ of the data were reserved for testing, and $20$ regression trees were grown, each fitted with different subsets of the training data and using all the features available. The $20$ predictions from each tree were aggregated using the sample mean and the Bayesian median of means.

As measures of test error for each dataset, consider the proportional reduction in MSE and MAD when using $\hat{\theta}_{\textnormal{BMM}}$:

[TABLE]

where $y_{i}$ is the realized outcome in the test data for a given dataset, $\hat{y}_{\text{SM},i}$ is the corresponding prediction using the sample mean for aggregation, and $\hat{y}_{\text{SM},i}$ is the prediction using the Bayesian median of means. While these measures try to address how well this modified bagging performs, the exercise must be viewed with caution since the underlying data distribution is not known.

Table 6 shows an estimate for the proportional reduction in MSE and MAd, for each dataset; larger entries favor the Bayesian median of means.

In most cases, the Bayesian median of means provides an improvement over traditional bagging both in terms of MSE and MAD, although the improvements are modest. Whenever there is decrease in performance, it is generally small compared to the gains. The largest improvements happens on datasets with a small number of samples, so the 20 estimates to be aggregated are much less stable. From Section 4.2.3, the smaller the variance or skewness of the estimates the more similar the procedure becomes to regular bagging. Overall, the gains in using the Bayesian median of means in datasets from a wide variety of areas displays its potential as an ensembling technique.

6 Conclusion

This paper introduced the Bayesian median of means, a non-parametric aggregation procedure that leads to estimates with relatively small variance at the expense of some bias. Asymptotically, the added bias is negligible with respect to the reduction in variance, much like maximum likelihood estimates in the parametric setting. Furthermore, the Bayesian median of means is asymptotically unbiased, and essentially reduces to the sample mean, a widely used location estimator, when the variance or skewness of the underlying sample is small. Computationally, it is easy to implement and can be parallelized in a straightforward way.

The Bayesian median of means is, however, a randomized procedure, so a deterministic approximation was developed, and dubbed the approximate Bayesian median of means. It resembles a shrinkage-type estimator, and has similar performance to the Bayesian median of means while being much faster to run.

Both methods were empirically tested on a variety of datasets and simulations. Overall, the Bayesian median of means was shown to be competitive with the sample mean in settings of low variance, and performed significantly better as the variance increased. When tested in real datasets, it displayed small but consistent gains over the sample mean, and its performance on applications such as importance sampling, cross-validation and bagging showed it can be useful in adding robustness to data resesampling schemes.

There are several directions for future work. First, more refined concentration bounds are likely possible, and with it a better understanding of the theoretical properties of the algorithm. Second, the approximate Bayesian median of means calls for an a theoretical investigation of its own, including issues of minimaxity and admissibility. Third, more extensive empirical analyses can illuminate the extent to which the procedure can overcome the sample mean, or fail altogether.

Lastly, this paper intends to encourage research into designing better general-purpose aggregation procedures, and stimulate further work in an area that has far-reaching applications in statistics and beyond.