Empirical priors and posterior concentration rates for a monotone density

TL;DR

This paper introduces a new empirical Bayes prior for monotone density estimation that achieves near-optimal posterior concentration rates and demonstrates practical advantages over traditional methods like Dirichlet process mixtures.

Contribution

The paper develops a novel empirical Bayes prior tailored for monotone densities, ensuring desirable concentration properties and near-minimax convergence rates.

Findings

Empirical Bayes prior achieves near-minimax posterior concentration rates.

Numerical results show improved performance over Dirichlet process mixtures.

Proposed method satisfies empirical prior concentration conditions.

Abstract

In a Bayesian context, prior specification for inference on monotone densities is conceptually straightforward, but proving posterior convergence theorems is complicated by the fact that desirable prior concentration properties often are not satisfied. In this paper, I first develop a new prior designed specifically to satisfy an empirical version of the prior concentration property, and then I give sufficient conditions on the prior inputs such that the corresponding empirical Bayes posterior concentrates around the true monotone density at nearly the optimal minimax rate. Numerical illustrations also reveal the practical benefits of the proposed empirical Bayes approach compared to Dirichlet process mixtures.

Figures

Tables

Table 1. Table 1: Exponential example—coverage probability and mean length of the 95% posterior credible regions for the empirical Bayes (EB) and Dirichlet process mixture model (DPM), for several values of x 𝑥 x and n 𝑛 n .

		Coverage Prob.		Mean Length
$n$	$x$	EB	DPM	EB	DPM
100	0.5	0.943	0.987	0.371	0.459
	1.0	0.984	0.965	0.270	0.338
	2.0	0.986	0.938	0.150	0.185
	3.0	0.990	0.891	0.086	0.098
200	0.5	0.934	0.981	0.285	0.371
	1.0	0.962	0.955	0.219	0.285
	2.0	0.983	0.914	0.120	0.154
	3.0	0.970	0.842	0.068	0.082

Table 2. Table 2: Half-normal example—same results as in Table  2 .

		Coverage Prob.		Mean Length
$n$	$x$	EB	DPM	EB	DPM
100	0.5	0.896	0.974	0.313	0.370
	1.0	0.969	0.961	0.291	0.361
	2.0	0.971	0.897	0.178	0.205
	3.0	0.986	0.636	0.054	0.043
200	0.5	0.865	0.991	0.242	0.294
	1.0	0.964	0.954	0.235	0.304
	2.0	0.956	0.885	0.137	0.166
	3.0	0.951	0.717	0.043	0.045

Equations

$$f_{\theta }(x)={\int }_0^{\infty }k(x\mid \mu )\theta (d\mu ),$$

$$\theta (d\mu )=\sum _{s=1}^S{\omega }_s{\delta }_{{\mu }_s}(d\mu ).$$

$$L_n(\theta )=\prod _{i=1}^n\sum _{s=1}^S{\omega }_{s}k(X_i\mid {\mu }_s)$$

$${\mathcal{L}}_n=\{\theta \in {\Theta }_n:L_n(\theta )\ge e^{-dn{\epsilon }_n^2}{L}_n(\hat{\theta })\},d>0,$$

$${\Theta }_n=\{\theta =(\omega ,\mu )\in \Delta (S_n)\times [t,T{]}^{S_n}\},0<t<T<\infty ,$$

$${\Pi }^n(d\theta )\propto L_n(\theta ){\Pi }_n(d\theta ),$$

$$S_n\propto {\epsilon }_n^{-1}=n{\epsilon }_n^2(\log n{)}^{-1},c_n\propto n{\epsilon }_n^{-2},\text{and}{\delta }_n\log (T_n/t_n)\lesssim \log n,$$

$$\mathsf{E}_{f^{\star}}\bigl{[}\Pi^{n}(\{f:d(f^{\star},f)>M\varepsilon_{n}\})\bigr{]}\to 0,\quad n\to\infty.$$

$$\text{there exists b,r>0 such that f⋆(x)≤e−bxr for all large x}.$$

$$\max\Bigl{\{}\int\Bigl{(}\log\frac{f^{\star}}{f_{\theta^{\dagger}}}\Bigr{)}f^{\star}\,dx,\,\int\Bigl{(}\log\frac{f^{\star}}{f_{\theta^{\dagger}}}\Bigr{)}^{2}f^{\star}\,dx\Bigr{\}}\leq\varepsilon_{n}^{2}.$$

$$\underset{n\to \infty }{\mathrm{liminf}}{e}^{Cn{\epsilon }_n^2}{\Pi }_n({\mathcal{L}}_n)>0,\text{with Pf⋆-probability 1}.$$

$$\int_{\Theta_{n}}\bigl{[}\mathsf{E}_{f^{\star}}\{\pi_{n}(\theta)^{p}\}\bigr{]}^{1/p}\,d\theta\leq e^{Kn\varepsilon_{n}^{2}}.$$

$$L_n(\theta )$$

$$=\sum _{(n_1,\dots ,n_S)}{\omega }_1^{n_1}\cdots {\omega }_S^{n_S}\sum _{(s_1,\dots ,s_n)}\prod _{s=1}^S{\mu }_s^{-n_s}1({\mu }_s\ge {\hat{X}}_s),$$

$$L_n(\theta )\ge \sum _{(n_1,\dots ,n_S)}{\omega }_1^{n_1}\cdots {\omega }_S^{n_S}\sum _{(s_1,\dots ,s_n)}\prod _{s=1}^S{\mu }_s^{-n_s}1({\hat{\mu }}_s\ge {\hat{X}}_s).$$

$$\mathsf{E}\{L_n(\theta )\}\ge \sum _{(n_1,\dots ,n_S)}\mathsf{E}({\omega }_1^{n_1}\cdots {\omega }_S^{n_S})\sum _{(s_1,\dots ,s_n)}\prod _{s=1}^S\mathsf{E}({\mu }_s^{-n_s})1({\hat{\mu }}_s\ge {\hat{X}}_s),$$

$$\mathsf{E}({\omega }_1^{n_1}\cdots {\omega }_S^{n_S})\ge \frac{\Gamma (c+S)c^n}{\Gamma (c+S+n)}{\hat{\omega }}_1^{n_1}\cdots {\hat{\omega }}_S^{n_S}.$$

$$\mathsf{E}({\mu }_s^{-n_s})={\int }_{{\hat{\mu }}_s}^{\infty }\frac{\delta {\hat{\mu }}_s^{\delta }}{{\mu }_s^{\delta +n_s+1}}d{\mu }_s=\frac{\delta }{\delta +n_s}{\hat{\mu }}_s^{-n_s}\ge \frac{1}{1+n{\delta }^{-1}}{\hat{\mu }}_s^{-n_s}.$$

$$\mathsf{E}\{L_n(\theta )\}\ge \frac{\Gamma (c+S)c^n}{\Gamma (c+S+n)}{e}^{-S\log n}{L}_n(\hat{\theta }),\text{all large n}.$$

$$\frac{\Gamma (c+S)c^n}{\Gamma (c+S+n)}\ge e^{-Kn{\epsilon }_n^2}.$$

$$\mathsf{E}\{L_n(\theta )\}/L_n(\hat{\theta })>e^{-Dn{\epsilon }_n^2}$$

$$\delta T^{\delta }{\mu }_s^{-(\delta +1)}1({\mu }_s\ge t),$$

$$\frac{(c+S)^{c+S+1/2}}{c^{c}}\,\Bigl{(}\frac{T}{t}\Bigr{)}^{\delta S}.$$

$$d(f^{\star },f)\le d(f^{\star },f_T^{\star })+d(f_T^{\star },f).$$

$$d(f^{\star },f_T^{\star })=\frac{1}{2}{\int }_0^{\infty }|f^{\star }(x)-f_T^{\star }(x)|dx=1-F^{\star }(T).$$

$$1-F^{\star }(T_n)\lesssim n^{-1}{\epsilon }_n,$$

$$d(f^{\star },f)>M{\epsilon }_n$$

$$⟹d(f_T^{\star },f)>M{\epsilon }_n-d(f^{\star },f_T^{\star })=M{\epsilon }_n-\{1-F^{\star }(T_n)\}$$

$$⟹d(f_T^{\star },f)>(M/2){\epsilon }_n,\text{say}.$$

$$A_n=\{f:d(f^{\star },f)>M{\epsilon }_n\}\text{and}{B}_n=\{f:d(f_{T_n}^{\star },f)>(M/2){\epsilon }_n\}.$$

$${\Pi }^n(A_n)\le {\Pi }^n(B_n).$$

$${\Pi }^n(B_n)={\Pi }^n(B_n)1({\mathcal{X}}_n)+{\Pi }^n(B_n)1({\mathcal{X}}_n^c)\le {\Pi }^n(B_n)1({\mathcal{X}}_n)+1({\mathcal{X}}_n^c),$$

Full text

Keywords and phrases: Density estimation; empirical Bayes; Grenander estimator; mixture model; shape constraint.

1 Introduction

Let $X_{1},\ldots,X_{n}$ be iid samples from a density function $f^{\star}$, supported on the positive half-line, assumed to be monotone non-increasing. Nonparametric inference on a monotone density has received considerable attention in the literature, dating back to Grenander, (1956), with a wide range of applications (e.g., Robertson et al., 1988; Groeneboom and Jongbloed, 2014). Theoretical properties of estimators have been studied in Prakasa Rao, (1969), Groeneboom, (1985), and Balabdaoui and Wellner, (2007), among others, with the behavior of the Grenander estimator at the origin being a now-classical example of inconsistency of the maximum likelihood estimator (Woodroofe and Sun, 1993) and failure of bootstrap (e.g., Kosorok, 2008; Sen et al., 2010).

From a Bayesian point of view, constructing a prior and corresponding posterior distribution for the monotone density is at least conceptually straightforward thanks to the mixture representation of Williamson, (1956); see Section 2. This makes it possible to construct priors for monotone densities using the standard tools, such as finite mixture models, Dirichlet processes, etc. However, theoretical analysis of the corresponding posterior distribution is complicated by the fact that, unless the support of $f^{\star}$ is known, the usual Kullback–Leibler property (e.g., Ghosal et al., 1999; Wu and Ghosal, 2008) used to prove posterior convergence results may not be satisfied; in fact, it could be the the Kullback–Leibler divergence of $f$ from $f^{\star}$ could be infinite for all $f$ in a set of prior probability 1. Therefore, the general theorems in, e.g., Ghosal et al., (2000) and Walker et al., (2007) cannot be applied. Salomond, (2014) worked around this to show, among other things, that the Bayesian posterior distribution based on various mixture priors has concentration rate within a logarithmic factor of the minimax optimal rate, $n^{-1/3}$, with respect to Hellinger or $L_{1}$ distance.

A recent trend in the Bayesian literature is asymptotic concentration results for empirical Bayes posteriors; see, e.g., Szabó et al., (2013), Rousseau and Szabo, (2017), and Donnet et al., (2018). These papers propose to extend the classical techniques and results to handle the case where the prior involves data in some way, e.g., through a plug-in estimator of a hyperparameter. However, given that the usual support conditions fail in the problem considered here, even with a fixed prior, it seems unlikely that these new techniques would apply to empirical Bayes monotone density estimation. Martin and Walker, (2017), building on Martin and Walker, (2014) and Martin et al., (2017), recently proposed a new empirical Bayes approach, one that constructs the empirical prior specifically so that the desirable posterior concentration rate properties are achieved. In particular, the empirical prior is designed to satisfy the prior support conditions—a variation on the Kullback–Leibler property—so this approach seems ideally suited for cases, like monotone density estimation, where satisfying the prior support condition is problematic. Here, in Section 2, I will construct a simple and intuitively appealing empirical prior and establish, in Sections 3–4, that the corresponding empirical Bayes posterior concentration rate is nearly minimax optimal. Beyond these desirable theoretical properties, in Section 5, I will show that the proposed empirical Bayes approach has a number of practical benefits compared to the Dirichlet process mixture model, including improved computational efficiency and finite-sample performance.

2 An empirical prior

The starting point here is the representation in Williamson, (1956) of a monotone density as a scale mixture of uniforms, i.e., for any monotone density density $f$, there exists a mixing distribution $\theta$, supported on a subset of $[0,\infty)$, such that $f=f_{\theta}$, where

[TABLE]

and the kernel $k(x\mid\mu)=\mu^{-1}1(x\leq\mu)$ is the $\mathsf{Unif}(0,\mu)$ density. From here, a prior for $f$ can be defined by introducing a prior for $\theta$ and using the mapping $\theta\mapsto f_{\theta}$.

As is typical, I will model $\theta$ as a (finite) discrete distribution, i.e.,

[TABLE]

This makes $f_{\theta}$ a finite mixture of uniforms. For the moment, fix the number of support points $S$. Then the mixing distribution can be expressed as a finite-dimensional parameter $\theta=(\omega,\mu)$, where $\omega=(\omega_{1},\ldots,\omega_{S})$ is the vector of mixture weights and $\mu=(\mu_{1},\ldots,\mu_{S})$ is the corresponding vector of mixture locations. The theory will require that $S=S_{n}$ be increasing with $n$ at a suitable rate; see Section 3. For the remainder of this section, I will focus on specifying a prior for $\theta=(\omega,\mu)$, given $S$. The prior here will be empirical in the sense that it depends on data in a particular way.

The general construction of an empirical prior in Martin and Walker, (2017) selects an appropriate data-driven center, e.g., the prior mode. Their motivation is to replace the usual Kullback–Leibler/prior concentration property with an “empirical” version. Towards this, write the likelihood function for $\theta=(\omega,\mu)$, with $S$ fixed, as

[TABLE]

If $\varepsilon_{n}$ is the target convergence rate and $\hat{\theta}$ is a maximizer of the likelihood $L_{n}$ over a suitable set $\Theta_{n}$, a sieve, then Martin and Walker, (2017) defined

[TABLE]

which is effectively a “neighborhood” of $\hat{\theta}$ in $\Theta_{n}$, an empirical or data-dependent version of the Kullback–Leibler neighborhood in classical Bayesian nonparametric studies (e.g., Schwartz, 1965). Like in the familiar Bayesian settings, the goal is for the prior to charge $\mathcal{L}_{n}$ with a sufficient amount of mass; see Condition LP in Section 4. But the fact that $\mathcal{L}_{n}$ is data-dependent means that the prior must also be so, thus, an empirical prior.

Here I take the prior mode equal to $\hat{\theta}=(\hat{\omega},\hat{\mu})$, a sieve maximum likelihood estimator. The particular sieve is of the form

[TABLE]

where $t=t_{n}$ and $T=T_{n}$ might also depend on $n$; see Section 3. Then the empirical prior for $\theta=(\omega,\mu)$ I propose here is as follows:

• •

$\omega$ and $\mu$ are independent;

• •

$\omega$ has an $S$-dimensional Dirichlet distribution, $\mathsf{Dir}_{S}(\hat{\alpha})$, on the simplex $\Delta(S)$, where $\hat{\alpha}_{s}=1+c\hat{\omega}_{s}$, $s=1,\ldots,S$, and $c=c_{n}$ is non-stochastic;

• •

$\mu_{1},\ldots,\mu_{S}$ are independent with $\mu_{s}\sim\mathsf{Par}(\hat{\mu}_{s},\delta)$, a Pareto distribution with scale parameter/lower bound $\hat{\mu}_{s}$ and non-stochastic shape parameter $\delta=\delta_{n}$.

It is easy to see that $\hat{\theta}=(\hat{\omega},\hat{\mu})$ is the mode. The Pareto prior for the components of $\mu$ is convenient because it is conjugate to the uniform mixture kernel. It is also important for the proofs in Section 4 that the prior for $\mu_{s}$ be supported on $[\hat{\mu}_{s},\infty)$, which is easily arranged with a Pareto distribution. The to-be-determined constants $(c,\delta)$ control the spread of the prior for $(\omega,\mu)$ around its mode $(\hat{\omega},\hat{\mu})$.

To summarize, $f$ is modeled as $f_{\theta}$ and an empirical prior on $f$ is induced by specifying an empirical prior for $\theta$ and using the mapping $\theta\mapsto f_{\theta}$. In what follows, $\Pi_{n}$ will denote the empirical prior for $\theta$ on the sieve (3) as described above. With a slight abuse of notation, I will also use $\Pi_{n}$ to denote the corresponding empirical prior for $f=f_{\theta}$; the meaning should be clear from the context. Given this prior, the corresponding posterior distribution $\Pi^{n}$ for $\theta$ is defined as

[TABLE]

where $L_{n}(\theta)$ is the likelihood function in (1). Again, with a slight abuse of notation, I will also write $\Pi^{n}$ for the empirical Bayes posterior for the monotone density $f$.

Remark 1*.*

The prior support condition eluded to above could be immediately achieved by taking the prior to be degenerate at the $\hat{\theta}$ that corresponds to Grenander’s estimator $\hat{f}=f_{\hat{\theta}}$, the nonparametric maximum likelihood estimator. Of course, the posterior based on this trivial empirical prior is also degenerate at $\hat{f}$ and, therefore, inherits the concentration rate of Grenander’s estimator. However, achieving the target rate is only a first objective. By using a non-degenerate prior, the posterior will have spread, leaving open the possibility for uncertainty quantification; see Sections 5–6.

Remark 2*.*

There is a clinical version of the model that is perhaps more natural for applications. In particular, when $n$ is large, the sieve ought to contain the $\theta$ corresponding to Grenander’s estimator, so practical applications could dispense with the sieves altogether—which eliminates the need to specify $S$ and $[t,T]$, and to maximize the likelihood over the sieve—and center the prior directly on Grenander’s estimator; see Section 5. However, establishing the concentration rate for this clinical version requires control on the mixture support size in Grenander’s estimator but, to my knowledge, no such results are available in the literature. Just like in Ghosal and van der Vaart, (2001, Sec. 4), a reasonable conjecture is that the sieve estimator above is the same as Grenander’s, in which case, the clinical version is also covered by Theorems 1–2.

3 Posterior concentration rate

The previous section described an empirical prior that, when combined with the likelihood via Bayes’s formula, leads to a posterior distribution $\Pi^{n}$ in (4) that can be used for inference on the monotone density $f$. But why is this a reasonable approach? To answer this question, I will provide conditions prior inputs—$c$, $\delta$, $S$, $t$, and $T$—such that the empirical Bayes posterior distribution $\Pi^{n}$ for $f$ concentrates around the true $f^{\star}$ at nearly the optimal minimax rate. Proofs of the two theorems are given in Section 4 and finite-sample performance of the posterior is investigated in Section 5.

Let $d$ denote the Hellinger or $L_{1}$ distance. Then the optimal rate with respect to $d$ is $n^{-1/3}$; see van der Vaart and Wellner, (1996, Theorem 2.7.5) and Ghosal et al., (2000, Example 3.2). Under certain conditions, this rate can be achieved, within a logarithmic factor, by the nonparametric maximum likelihood estimator (e.g., Birgé, 1989; Balabdaoui and Wellner, 2007) and by certain nonparametric Bayesian methods (Salomond, 2014). The following theorem establishes the near-optimal concentration rate for the proposed empirical Bayes approach, in the case where $f^{\star}$ has a bounded support.

Theorem 1**.**

Let the true density $f^{\star}$ be monotone non-increasing with support $[0,T^{\star}]$, with $f^{\star}(0)<\infty$, and let $\varepsilon_{n}=(\log n)^{1/3}n^{-1/3}$ be the target rate. If the prior inputs $(c,\delta,S,t,T)=(c_{n},\delta_{n},S_{n},t_{n},T_{n})$ satisfy

[TABLE]

where $T_{n}$ and $\delta_{n}$ are non-decreasing and $t_{n}$ is non-increasing, then there exists a constant $M>0$ such that the posterior distribution $\Pi^{n}$ satisfies

[TABLE]

The conclusion here is similar to that in Theorem 1 of Salomond, (2014). Indeed, the concentration rate above is the same as that obtained by a suitable Dirichlet process mixture model, which is minimax optimal up to the logarithmic factor.

The latter condition in (5) deserves some explanation. Since the upper bound $T^{\star}$ is finite, any sufficiently large $T_{n}$ would suffice for estimating $f^{\star}(x)$ for $x$ near $T^{\star}$. Similarly, since a draw $f$ from the proposed prior satisfies $f(0)\leq t_{n}^{-1}$, any $t_{n}$ less than $f^{\star}(0)^{-1}>0$ would suffice for estimating $f^{\star}(x)$ for $x$ near 0. Therefore, $\log(T_{n}/t_{n})$ could be very slowly increasing, or even bounded, which means $\delta_{n}$ can grow as fast as order $\log n$. Intuitively, slowly increasing $\log(T_{n}/t_{n})$ indicates some certainty about the support of $f^{\star}$, in which case, a larger $\delta_{n}$ and, hence, a smaller Pareto variance, is reasonable. On the other hand, if the support is uncertain, one could take $T_{n}/t_{n}$ to be polynomial in $n$, in which case $\delta_{n}$ must be bounded and, hence, the Pareto variance must be bounded away from zero.

The only serious assumption on $f^{\star}$ in Theorem 1 is that the support is bounded. It turns out that this bounded-support condition can be replaced by a condition on the tails of $f^{\star}$. Condition C4 in Salomond, (2014) states:

[TABLE]

The next result is analogous to Theorem 2 in Salomond, (2014).

Theorem 2**.**

Let the true density $f^{\star}$ be monotone non-increasing, with $f^{\star}(0)<\infty$, whose support is $[0,\infty)$. Assume that $f^{\star}$ satisfies (7) for a given $r$, and set the target rate equal to $\varepsilon_{n}=(\log n)^{1/3+1/r}n^{-1/3}$. Let the prior inputs be as in (5), but with $T_{n}\gtrsim(\log n)^{1/r}$. Then the conclusion (6) of Theorem 1 holds with the modified rate $\varepsilon_{n}$.

Note that the rate in the unbounded support case is slightly slower, by only a logarithmic factor, than in the bounded support case. Actually, Theorem 1 can be viewed as a special case of this Theorem 2: if the support is bounded, then (7) holds for “$r=\infty$” and so the rate in Theorem 2 agrees with that in Theorem 1. The only subtlety is that the choice of $T=T_{n}$ depends on $r$, a feature of $f^{\star}$, which is typically unknown in real examples. If one is willing to assume a positive lower bound $r_{0}$ on $r$, then $T=T_{n}$ can be chosen with $r=r_{0}$. It is not known if the rate in Theorem 2 is optimal so, even though working with a lower bound on $r_{0}$—corresponding to a larger $T_{n}$—will slow down the rate slightly, there is no practical difference compared to the rate with the true $r$. Salomond, (2014) does not say whether his Theorem 2 requires knowledge of the tail exponent $r$, but the tails of the Dirichlet process base measure generally affect the posterior concentration rates (e.g., Ghosal and van der Vaart, 2001, Theorem 5.1). So if a target rate depending on $r$ is to be achieved, then this requires some $r$-dependent condition on the base measure and, therefore, to check this condition, $r$ must be known.

4 Proofs

4.1 General strategy

Martin and Walker, (2017) proposed a general strategy for constructing empirical priors such that the corresponding posterior distribution has the desired concentration properties. Their Theorem 1 lists three general conditions that include assumptions about the prior concentration, one local and one global, as well as an assumption about the approximation properties of the sieve. I will summarize these conditions in the context of iid data as being considered here. Let $\varepsilon_{n}$ be the target rate.

Condition S*.*

There exists a $\theta^{\dagger}=\theta_{n}^{\dagger}$ in the sieve $\Theta_{n}$ such that

[TABLE]

Condition LP*.*

For a given $d>0$ define $\mathcal{L}_{n}$ as in (2). Then there exists a constant $C>0$ such that the empirical prior $\Pi_{n}$ satisfies

[TABLE]

Condition GP*.*

Let $\pi_{n}$ be the density function for $\theta$ under the empirical prior. For a constant $p>1$, there exists $K>0$ such that

[TABLE]

4.2 Proof of Theorem 1

I will begin by checking Condition LP. For fixed $S$, if $\theta=(\omega,\mu)$ denotes the mixture weights and locations, respectively, then the likelihood function $L_{n}(\theta)$ in (1) for the discrete mixture model can be expressed as

[TABLE]

where $\hat{X}_{s}=\max_{i:s_{i}=s}X_{i}$ is the largest of the $n_{s}$ $X$ values in category $s$, relative to the partition determined by $(s_{1},\ldots,s_{n})$ with frequency table $(n_{1},\ldots,n_{S})$; the inner- and outer-most sums are over all such partitions and frequency tables, respectively. Since the prior for $\mu_{s}$ is supported on $[\hat{\mu}_{s},\infty)$, I can lower-bound the likelihood by

[TABLE]

The prior has $\omega$ and $\mu$ independent, and $\mu_{1},\ldots,\mu_{S}$ independent, so

[TABLE]

where expectation is with respect to the prior for $\theta=(\omega,\mu)$. The proof of Proposition 2 in Martin and Walker, (2017) gives a bound for the first expectation, i.e.,

[TABLE]

For the Pareto prior, $\mathsf{Par}(\hat{\mu}_{s},\delta)$, on $\mu_{s}$, we have

[TABLE]

Since $\delta=\delta_{n}$ is non-vanishing, the first term in the lower bound is at least $O(n^{-1})$. Plugging this bound back into the above expectation gives

[TABLE]

As in Proposition 2 of Martin and Walker, (2017), if $S=S_{n}$ is of the order $n\varepsilon_{n}^{2}(\log n)^{-1}$ and $c=n\varepsilon_{n}^{-2}$, then there exists a constant $K>0$ such that

[TABLE]

Similarly, for the second term, since $S\log n$ is of the order $n\varepsilon_{n}^{2}$, we can conclude that, for a suitable constant $D>0$,

[TABLE]

which, according to the argument in the proof of Proposition 3 in Martin and Walker, (2017), implies Condition LP.

Next, I check Condition GP. As a first step, we have that the density for the Dirichlet prior on $\omega$ is uniformly upper bounded by $(c+S)^{c+S+1/2}c^{-c}$, which does not depend on data. For the prior on $\mu$, the density function is upper bounded by

[TABLE]

which is also free of data. Then, for any $p>1$, the integral (8) is bounded by

[TABLE]

With $S=S_{n}$ of the order $n\varepsilon_{n}^{2}(\log n)^{-1}$ and $\delta\log(T/t)$ of order $\log n$, it follows that the second term in (9) is bounded by $e^{An\varepsilon_{n}^{2}}$. Similarly, for $c=n\varepsilon_{n}^{2}$ and $S=S_{n}$ of the order $n\varepsilon_{n}^{2}(\log n)^{-1}$, Martin and Walker, (2017) showed that the first term in (9) is also $e^{Bn\varepsilon_{n}^{2}}$ so, altogether, the relevant integral is bounded by $e^{Cn\varepsilon_{n}^{2}}$, hence Condition GP.

Finally, note that, for the case of $f^{\star}$ with bounded support $[0,T^{\star}]$ and $f^{\star}(0)<\infty$, Condition S on the sieve follows from Lemma 11 in Salomond, (2014). My sieve has a lower bound $t>0$ but, if it is small enough, then it does not affect Salomond’s calculations since there is no benefit to having a mixture location smaller than $f^{\star}(0)^{-1}>0$. Having checked the three conditions in Section 4.1, the conclusion of Theorem 1 follows from the general results in Martin and Walker, (2017).

4.3 Proof of Theorem 2

For a given $T$, possibly depending on $n$, write $f_{T}^{\star}$ for the normalized version of $f^{\star}$ to the interval $[0,T]$, i.e., $f_{T}^{\star}(x)=f^{\star}(x)/F^{\star}(T)$, where $F^{\star}$ is the distribution function corresponding to $f^{\star}$. Without loss of generality, let $d$ denote the $L_{1}$ distance. Then the triangle inequality implies that, for any density $f$,

[TABLE]

Moreover, a simple calculation shows that

[TABLE]

Under the condition (7) on the density $f^{\star}$, it is easy to check that the tail probability $1-F^{\star}(T)\lesssim\Gamma(r^{-1},T^{r})$, where $\Gamma(s,t)$ is the upper incomplete gamma function, i.e., $\Gamma(s,t)=\int_{t}^{\infty}y^{s-1}e^{-y}\,dy$. From the well-known asymptotic behavior of this gamma function, it follows that if $T_{n}\gtrsim(\log n)^{1/r}$, then

[TABLE]

where $\varepsilon_{n}=(\log n)^{1/3+1/r}n^{-1/3}$ as in the statement of the theorem; see, also, page 1390 in Salomond, (2014). Therefore, for any $f$,

[TABLE]

This effectively converts the problem into one with bounded support. To see this, define the two sets of densities

[TABLE]

Then the argument above implies that $A_{n}\subseteq B_{n}$ which, in turn, implies

[TABLE]

Define the event $\mathscr{X}_{n}=\{(x_{1},\ldots,x_{n})\in[0,\infty)^{n}:x_{(n)}\leq T_{n}\}$, where $x_{(n)}=\max_{i}x_{i}$. Based on the bound in (10), it is easy to check that $\mathsf{P}_{f^{\star}}(\mathscr{X}_{n}^{c})=o(1)$. Next, write

[TABLE]

so that

[TABLE]

The expectation on the right-hand side can be rewritten as

[TABLE]

Since $\mathsf{P}_{f^{\star}}(\mathscr{X}_{n})\to 1$, it remains to deal with the conditional expectation. The key observation is that the conditional distribution of $(X_{1},\ldots,X_{n})$, given $\mathscr{X}_{n}$, is iid $f_{T_{n}}^{\star}$ and, therefore,

[TABLE]

hence, the claim that this is effectively a bounded support problem. Moreover, all of the work in checking Conditions LP, GP, and S in the proof of Theorem 1 above applies here with $f^{\star}$ replaced by $f_{T_{n}}^{\star}$ and the modified rate. Since the general results in Martin and Walker, (2017) do not require that the “true parameter” be fixed, we can conclude that the right-hand side of (11) is $o(1)$ and, therefore,

[TABLE]

as was to be shown. This completes the proof of Theorem 2.

5 Numerical results

5.1 Computation

Here I will focus on the clinical version of the empirical prior suggested in Remark 2. That is, this version of the prior is centered on Grenander’s estimator, which is obtained from the output provided by the function grenander in the R package “fdrtool” (Klaus and Strimmer, 2015); more precisely, the clinical empirical prior is centered on the $\hat{\theta}=(\hat{\omega},\hat{\mu})$ for which $f_{\hat{\theta}}$ is Grenander’s estimator. For moderate to large $n$, there ought to be no difference between this and the empirical prior in Section 2, but the former has two clear advantages: first, there is no need to specify $S$ or the interval $[t,T]$; second, maximizing the likelihood over the sieve $\Theta_{n}$ is not so easy, but there is already efficient software available for computing Grenander’s estimator.

Since $S$ and $[t,T]$ are taken care of by the clinical formulation, it only remains to specify $c=c_{n}$ and $\delta=\delta_{n}$. In what follows, I will use

[TABLE]

a choice guided by the conditions in Theorems 1–2. The simulation experiments below suggest that this choice works well in various cases, though more work would be needed to determine if these values are “good” in any general sense; see Section 6.

I will compare the empirical Bayes results with those from a Bayesian Dirichlet process mixture model. I take the Dirichlet process precision parameter as 1 and the base measure as an inverse gamma with data-driven choice of shape and scale parameters, $a$ and $b$. That is, I choose $a$ and $b$ so that the base measure mean and variance match the mean and variance of the data; this is to ensure that the base measure center and spread are consistent with the data. Other base measures could also be considered, such as gamma. A Pareto base measure would be a natural choice, given its conjugacy to the uniform kernel, but this does not satisfy the tail conditions in Salomond, (2014) needed to achieve the nearly-$n^{-1/3}$ concentration rate. Conjugacy would allow for very convenient posterior sampling via the slice sampler in Walker, (2007) and Kalli et al., (2011). But for other base measures, such as gamma and inverse gamma, rejection sampling is required in every MCMC iteration, so the computation is slower than with a Pareto base measure. Note that the proposed empirical Bayes approach can employ the computationally efficient Pareto prior since the strategic data-driven centering makes the tails irrelevant to the concentration rate properties. R codes for both the Bayes and empirical Bayes methods are available at https://www4.stat.ncsu.edu/~rmartin/.

5.2 Simulated data

Consider the case where $f^{\star}(x)=e^{-x}$, $x\geq 0$, is an exponential density. Figure 1 plots posterior samples from the respective posterior distributions, the posterior mean density, the 95% posterior credible band, the Grenander estimator, and the truth. Both posteriors generally follow the Grenander estimator, and both have roughly the same spread but, at least in this case, the empirical Bayes posterior is less rough and better captures the shape of the truth. The credible bands for both posteriors cover the truth at most $x$ values, but the Dirichlet process mixture misses in a few places because of its roughness. Another difference, not apparent from the plots, is that the empirical Bayes computations take about 50% less time than the Dirichlet process mixture.

Next I present a more detailed investigation into the concentration properties of the two posteriors. Specifically, I will consider the coverage properties of the 95% credible bands at select values of $x$. That is, for a fixed $x$, I get posterior samples of $f(x)$ and return the 0.025 and 0.975 quantiles as the 95% posterior credible interval. I repeat this process for 500 data sets and report the (Monte Carlo approximation of the) coverage probability and expected lengths for the two posteriors. Table 1 reports (Monte Carlo approximations of) the coverage probability and mean lengths for the 95% credible intervals based on 500 samples from an exponential distribution, where $f^{\star}(x)=e^{-x}$, $x\geq 0$. Table 2 reports the same but for a half-normal distribution, where $f^{\star}(x)=2\,\mathsf{N}(x\mid 0,1)$, $x\geq 0$.

The key observation across all the different setting is that the 95% posterior credible intervals from the Dirichlet process mixture model tend to under-cover, i.e., the coverage probability is less than the nominal 0.95, sometimes much less, while those from the empirical Bayes approach tend to over-cover. However, the higher coverage of the empirical Bayes intervals is not a result of being wider on average; in fact, the empirical Bayes credible intervals tend to be shorter than the Bayesian competitor’s. Coverage properties of nonparametric methods is a delicate matter and a detailed investigation is beyond the scope of this paper, but these results strongly suggest that the proposed empirical Bayes procedure—despite its “double-use” of the data—does not follow the data too closely and, as suggested in Remark 1, may provide valid uncertainty quantification.

5.3 Real data

Example 1*.*

Silverman, (1986, Table 2.1) presents data on the lengths of psychiatric treatment undergone by $n=86$ patients used as controls in a study of suicide risks. Figure 2 shows the data histogram, the posterior mean density, and the 95% credible bands for the two estimators, based on the same settings described in Section 5.1. As expected, both posterior means fit the data well, but the empirical Bayes estimate is more smooth. That the empirical Bayes credible band is also a bit narrower than that of the Dirichlet process mixture should not be a concern based on the simulation results above.

Example 2*.*

The Norwegian fire claims data is a common example in the actuarial science literature (e.g., Brazauskas and Kleefeld, 2016). I consider $n=820$ fire loss claims exceeding $500$ thousand Norwegian krones during the year 1988. Figure 3 shows the data, the posterior mean densities, and the 95% credible bands for the two methods. In this case, the empirical Bayes posterior mean is a bit more rough than that of the Dirichlet process mixture, but arguably fits the data histogram better. And the narrow credible band is expected since this data set is about ten times as large as that in Example 1.

6 Conclusion

This paper presents a unique approach to the specification of an empirical or data-dependent prior for nonparametric Bayesian-like inference on a density function. The chief novelty is the centering of the prior on a suitable estimator, and this is particularly advantageous in the present context where the usual Kullback–Leibler property may fail. The challenge is to choose the prior tails so that the corresponding posterior does not track the data too closely, and Theorems 1–2 provide sufficient conditions on the prior inputs to achieve the target posterior concentration rate. Whether the proposed formulation can achieve valid uncertainty quantification, e.g., in the sense of Szabó et al., (2015), remains an open question, but the numerical illustrations in Section 5 are promising.

Acknowledgment

This work is partially supported by the National Science Foundation, DMS–1737933. The author also thanks two anonymous the reviewers for their helpful feedback.