A Bayesian nonparametric approach to log-concave density estimation

TL;DR

This paper introduces a Bayesian nonparametric method for estimating log-concave densities that achieves near-minimax convergence rates, with practical computational approaches demonstrated through simulations.

Contribution

It proposes a novel Bayesian approach using an exponentiated Dirichlet process mixture prior for log-concave density estimation, with theoretical convergence guarantees.

Findings

Posterior converges at near-minimax rate in Hellinger distance.

Two computationally feasible approximations are developed.

Empirical Bayes approach is demonstrated through simulations.

Abstract

The estimation of a log-concave density on $\mathbb{R}$ is a canonical problem in the area of shape-constrained nonparametric inference. We present a Bayesian nonparametric approach to this problem based on an exponentiated Dirichlet process mixture prior and show that the posterior distribution converges to the log-concave truth at the (near-) minimax rate in Hellinger distance. Our proof proceeds by establishing a general contraction result based on the log-concave maximum likelihood estimator that prevents the need for further metric entropy calculations. We also present two computationally more feasible approximations and a more practical empirical Bayes approach, which are illustrated numerically via simulations.

Figures

Tables

Table 1. Table 1: Frequencies out of 100 experiments when the empirical Bayes credible set contained the true function values at points x = 0.5 , 1 , … , 3 𝑥 0.5 1 … 3 x=0.5,1,...,3 . From top to bottom the sample size increases from n = 50 𝑛 50 n=50 until n = 500 𝑛 500 n=500 .

n \ x	0.5	1	1.5	2	2.5	3
50	0.51	0.81	0.88	0.85	0.87	0.93
200	0.68	0.97	0.94	0.87	0.94	0.97
500	0.83	0.96	0.88	0.87	0.86	0.91

Equations

$$w_{-}^{\prime }(x)=\underset{s\nearrow x}{\lim }{w}^{\prime }(s)\text{and}{w}_{+}^{\prime }(x)=\underset{s\searrow x}{\lim }{w}^{\prime }(s).$$

$${\mathcal{F}}_{\alpha ,\beta }:=\{f\in \mathcal{F}:f(x)\le e^{\beta -\alpha |x|}\forall x\in \mathbb{R}\}.$$

$$\Pi_{n}\Big{(}f\in\mathcal{F}:\int_{\mathbb{R}}f_{0}\Big{(}\log\frac{f_{0}}{f}\Big{)}\leq\varepsilon_{n}^{2},\quad\int_{\mathbb{R}}f_{0}\Big{(}\log\frac{f_{0}}{f}\Big{)}^{2}\leq\varepsilon_{n}^{2}\Big{)}\geq\exp(-Cn\varepsilon_{n}^{2}).$$

$${\Pi }_n(f\in \mathcal{F}:h(f,f_0)\ge M{\epsilon }_n|X_1,...,X_n)\to 0$$

$$f_w(x)=\frac{e^{w(x)}}{{\int }_{\mathbb{R}}{e}^{w(y)}dy},$$

$$f(x)={\int }_x^{\infty }\frac{1}{u}dP(u),$$

$$w(x)={\gamma }_1{\int }_0^{\infty }\frac{u\wedge x}{u}dP(u)-{\gamma }_{2}x,$$

$$W(x)={\gamma }_1{\int }_0^{b_n-a_n}\frac{u\wedge (x-a_n)}{u}dP(u)-{\gamma }_2(x-a_n),$$

$$\Pi (f:h(f,f_0)\ge M(\log n)n^{-2/5}\mid X_1,...,X_n)\to 0$$

$${\Pi }_n(f:h(f,f_0)\ge M{\epsilon }_n\mid X_1,...,X_n)\to 0$$

$${\epsilon }_n=\max \left((\log n)n^{-2/5},(b_n-a_n)n^{-4/5}\right).$$

$$\pi (a,b)\ge C{e}^{-c_1|a{|}^q-c_2|b-a{|}^r}\text{for all }a<b$$

$${\Pi }_n(f:h(f,f_0)\ge M(\log n)n^{-2/5}\mid X_1,...,X_n)\to 0$$

$$\begin{array}{r}{\theta }_i\stackrel{iid}{\sim }\bar{H}\text{mathds}{1}_{[0,b_n-a_n]},\text{for i=1,...,N},\\ p=(p_1,...p_N)\sim Dir({\alpha }_1,....,{\alpha }_N),\\ {\gamma }_i\stackrel{iid}{\sim }{p}_{{\gamma }_i},i=1,2,\\ f_{\theta ,p,{\gamma }_1,{\gamma }_2}(x)=\frac{\exp \{{\gamma }_1{\sum }_{i=1}^N\frac{{\theta }_i\wedge (x-a_n)}{{\theta }_i}{p}_i-{\gamma }_2(x-a_n)\}\text{mathds}{1}_{[a_n,b_n]}(x)}{{\int }_{a_n}^{b_n}\exp \{{\gamma }_1{\sum }_{i=1}^N\frac{{\theta }_i\wedge (u-a_n)}{{\theta }_i}{p}_i-{\gamma }_2(u-a_n)\}du},\end{array}$$

$$p_i\sim V_i\prod _{j=1}^{i-1}(1-V_j),\text{where Vi∼Beta(1,H(R+)), i=1,...,N−1.}$$

$${\Pi }_n^{\prime }(f:h(f,f_0)\ge M{\epsilon }_n\mid X_1,...,X_n)\to 0$$

$${\bar{f}}_n(x)=\exp \left({\bar{\gamma }}_1\sum _{i=1}^{\bar{N}}\frac{z_i\wedge (x-a_n)}{z_i}{\bar{p}}_i-{\bar{\gamma }}_2(x-a_n)+{\bar{\gamma }}_3\right)\text{mathds}{1}_{[a_n,b_n]}(x).$$

$$\overline{\mathcal{F}}^{\xi,\eta}=\bigg{\{}f\in\mathcal{F}:\mu_{f}:=\int xf(x)dx=\xi,\quad\sigma_{f}^{2}:=\int(x-\mu_{f})^{2}f(x)dx=\eta\bigg{\}}$$

$${\tilde{\mathcal{F}}}^{\xi ,\eta }=\{f\in \mathcal{F}:|{\mu }_f|\le \xi ,|{\sigma }_f^2-1|\le \eta \}.$$

$$\sup_{g_{0}\in\overline{\mathcal{F}}^{0,1}}P_{g_{0}}^{n}\Big{(}h(\hat{g}_{n},g_{0})\geq Ln^{-2/5}\Big{)}\leq C\exp\left(-cn^{1/(4+2\varepsilon)}\right)+D\exp\left(-dL^{2}n^{1/5}\right),$$

$$P_{f_0}^n{\varphi }_n\to 0,\text{ and }\underset{f\in \mathcal{F}:h(f,f_0)\ge M{\epsilon }_n}{\sup }{P}_f^n(1-{\varphi }_n)\le e^{-(C+4)n{\epsilon }_n^2}$$

$$P_{f_0}^n{\varphi }_n=P_{g_0}^n(h({\hat{g}}_n,g_0)\ge M_0{\epsilon }_n)\le P_{g_0}^n(h({\hat{g}}_n,g_0)\ge L_0{n}^{-2/5})\le C{e}^{-c{n}^{1/5}}\to 0$$

$$P_f^n(1-{\varphi }_n)$$

$$\le P_f^n(h(f_0,f)-h(f,{\hat{f}}_n)<M_0{\epsilon }_n)$$

$$\le P_f^n((M-M_0){\epsilon }_n<h(f,{\hat{f}}_n)).$$

$$\underset{f\in \mathcal{F}:h(f,f_0)\ge M{\epsilon }_n}{\sup }{P}_f^n(1-{\varphi }_n)$$

$$=\underset{g\in {\overline{\mathcal{F}}}^{0,1}}{\sup }{P}_g^n(h({\hat{g}}_n,g)>L_n{n}^{-2/5})$$

$$\le C(\rho )\exp \left(-c(\rho )n^{\frac{1}{4+2\epsilon (\rho )}}\right)+D\exp \left(-d{L}_n^2{n}^{1/5}\right)$$

$$\underset{f\in \mathcal{F}:h(f,f_0)\ge M{\epsilon }_n}{\sup }{P}_f^n(1-{\varphi }_n)$$

$$\begin{split}h^{2}(f_{0},f_{0,n})&\leq 2\int_{\mathbb{R}}f_{0}\Bigg{(}1-\frac{1}{\sqrt{\int f_{0}\psi_{n}}}\Bigg{)}^{2}dx+\frac{2}{\int f_{0}\psi_{n}}\int_{\mathbb{R}}f_{0}(1-\sqrt{\psi_{n}})^{2}dx\\ &\leq 2\frac{(\int f_{0}\psi_{n}-1)^{2}}{\int f_{0}\psi_{n}}+\frac{2}{\int f_{0}\psi_{n}}\int_{\mathbb{R}\backslash[-t_{n},t_{n}]}f_{0}dx=o(n^{-1}).\end{split}$$

$$\displaystyle\int_{\mathbb{R}}f_{0,n}\Big{(}\log\frac{f_{0,n}}{f_{W}}\Big{)}^{k}$$

$$\times {\left(1+\log {\Vert \frac{f_{0,n}}{{\bar{f}}_n}\Vert }_{L^{\infty }([a_n,b_n])}+\log {\Vert \frac{{\bar{f}}_n}{f_W}\Vert }_{L^{\infty }([a_n,b_n])}\right)}^k$$

$${\bar{f}}_n(x)=\exp \left({\bar{\gamma }}_1\sum _{i=1}^{\bar{N}}\frac{z_i\wedge (x-a_n)}{z_i}{\bar{p}}_i-{\bar{\gamma }}_2(x-a_n)+{\bar{\gamma }}_3\right)\text{mathds}{1}_{[a_n,b_n]}(x),$$

Full text

AMS 2000 subject classifications: 62G07, 62G20.

Keywords: Density estimation, log-concavity, Dirichlet mixture, posterior distribution, convergence rate, nonparametric hypothesis testing.

1 Introduction

Nonparametric shape constraints offer practitioners considerable modelling flexibility by providing infinite-dimensional families that cover a wide range of parameters whilst also including numerous common parametric families. Log-concave densities on $\mathbb{R}$, that is densities whose logarithm is a concave function taking values in $[-\infty,\infty)$, constitute a particularly important shape-constrained class. This class includes many well-known parametric densities that are frequently used in statistical modelling, including the Gaussian, uniform, Laplace, Gumbel, logistic, gamma distributions with shape parameter at least one, Beta$(\alpha,\beta)$ distributions with $\alpha,\beta\geq 1$ and Weibull distributions with parameter at least one.

One of the original statistical motivations for considering log-concave density estimation was the problem of estimating a unimodal density with unknown mode. While this is a natural constraint in many applications, the nonparametric MLE over this class does not exist [3]. Since the class of log-concave densities equals the class of strongly unimodal densities [20], Walther [38] argues that this class provides a natural alternative to the full set of all unimodal densities. The class of log-concave densities also preserves many of the attractive properties of Gaussian distributions, such as closure under convolution, marginalization, conditioning and taking products. One can therefore view log-concave densities as a natural infinite-dimensional surrogate for Gaussians that retain many of their important features yet allow substantially more freedom, such as heavier tails. For these reasons, estimation of log-concave densities has received significant attention in recent years, particularly concerning the performance of the log-concave MLE [9, 31, 4, 5, 10, 7, 23, 22].

Outside density estimation, log-concavity as a modelling assumption has found applications in many statistical problems, such as mixture models [38, 1], tail index estimation [26], clustering [5], regression [10] and independent component analysis [30]. For general reviews of inference with log-concave distributions and estimation under shape constraints, see [39] and [16] respectively.

The Bayesian approach provides a natural way to encode shape constraints via the prior distribution, for instance under monotonicity [32, 21, 29, 28] or convexity contraints [33, 18, 17]. We present here a Bayesian nonparametric method for log-concave density estimation on $\mathbb{R}$ based on an exponentiated Dirichlet process mixture prior, which we show converges to a log-concave truth in Hellinger distance at the (near-)minimax rate. To the best of our knowledge, this is the first Bayesian nonparametric approach to this problem. We also study two computationally motivated approximations to the full Dirichlet process mixture based on standard Dirichlet process approximations, namely the Dirichlet multinomial distribution and truncating the stick-breaking representation (see Chapter 4.3.3 of [14]). We further propose both an empirical and hierarchical Bayes approach that have clear practical advantages, while behaving similarly to the above in simulations. All of these priors are easily implementable using a random walk Metropolis-Hastings within Gibbs sampling algorithm, which we illustrate in Section 3.

An advantage of the Bayesian method is that point estimates and credible sets can be approximately computed as soon as one is able to sample from the posterior distribution. In particular, the posterior yields easy access to statements on Bayesian uncertainty quantification as we show numerically in Section 3. Our numerical results suggest that pointwise credible sets have reasonable coverage at moderate sample sizes.

The Bayesian approach also permits inference about multiple quantities, such as functionals, in a unified way using the posterior distribution. A particular functional of interest is the mode of a log-concave density. While the pointwise limiting distribution of the log-concave MLE is known [2], it depends in a complicated way on the unknown density making it difficult to use to construct a confidence interval for the mode. An alternative approach to constructing a confidence interval based on comparing the log-concave MLE with the mode constrained MLE has recently been proposed [6]. For the Bayesian, the marginal posterior of the mode provides a natural approach to both estimation and uncertainty quantification. Indeed, it is easy to construct Bayesian credible intervals as we demonstrate numerically in Section 3. Whether such an approach is theoretically justified from a frequentist perspective is a subtle question related to the semiparametric Bernstein-von Mises phenomenon (Chapter 12 of [14]) that is, however, beyond the scope of this article. We also note that other constraints, such as a known mode [8], can similarly be enforced through suitable prior calibration.

Given the good performance of the log-concave MLE, one might expect that Bayesian procedures, being driven by the likelihood, behave similarly well. This is indeed the case, as we show below. Our proof relies on the classic testing approach of Ghosal et al. [12] with interesting modifications in the log-concave setting. The existence and optimality of the MLE in Hellinger distance is closely linked to a uniform control of bracketing entropy [35]. In our setting, one can exploit the affine equivariance of the log-concave MLE (Remark 2.4 of [10]) to circumvent the need to control the metric entropy of the whole space by reducing the problem to studying a subset satisfying restrictions on the first two moments of the underlying density. This is a substantial reduction, since obtaining sharp entropy bounds in even this reduced case is highly technical, see Theorem 4 of Kim and Samworth [23]. One can then use the MLE to construct suitable plug-in tests with exponentially decaying type-II errors as in Giné and Nickl [15] that take full advantage of the extra structure of the problem compared to the standard Le Cam-Birgé testing theory for the Hellinger distance [24]. Indeed, a naive attempt to control the entropy directly, as is standard in the Bayesian nonparametrics literature (e.g. [12]), results in an overly small set on which the prior must place most of its mass. This leads to unnecessary restrictions on the prior, which in particular are not satisfied by the priors we consider in Section 2, see Remark 1. Beyond this, there remain significant technical hurdles to proving that the prior places sufficient mass in a Kullback-Leibler neighbourhood of the truth, in particular related to the approximation of log-concave densities using piecewise log-linear functions with suitably spaced knots.

The paper is structured as follows. In Section 2 we introduce our priors and present our main results, both on general contraction for log-concave densities and for the specific priors considered here. In Section 3 we present a simulation study, including a more practical empirical Bayes implementation, with some discussion in Section 4. In Section 5 we present the proofs of the main results with technical results placed in Section 6. Some additional simulations are found in Section 7.

Notation: For two probability densities $p$ and $q$ with respect to Lebesgue measure $\lambda$ on $\mathbb{R}$, we write $h^{2}(p,q)=\int(\sqrt{p}-\sqrt{q})^{2}$ for the squared Hellinger distance, $K(p,q)=\int p\log\tfrac{p}{q}$ for the Kullback-Leibler divergence and $V=\int p(\log\tfrac{p}{q})^{2}$. We denote by $P_{f_{0}}^{n}$ the product probability measure corresponding to the joint distribution of i.i.d random variables $X_{1},\dots,X_{n}$ with density $f_{0}$ and write $P_{f_{0}}=P_{f_{0}}^{1}$. For a function $w$, we denote by $w_{-}^{\prime}$ and $w_{+}^{\prime}$ its left and right derivatives respectively, that is

[TABLE]

Let $\mathbb{R}^{+}=[0,\infty)$ and for two real numbers $a,b$, let $a\wedge b$ and $a\vee b$ denote the minimum and maximum of $a$ and $b$ respectively. Finally, the symbols $\lesssim$ and $\gtrsim$ stand for an inequality up to a constant multiple, where the constant is universal or (at least) unimportant for our purposes.

2 Main Results

Consider i.i.d. density estimation, where we observe $X_{1},...,X_{n}\sim f_{0}$ with $f_{0}=e^{w_{0}}$ an unknown log-concave density to be estimated. Let $\mathcal{F}$ denote the class of upper semi-continuous log-concave probability densities on $\mathbb{R}$. For $\alpha>0$ and $\beta\in\mathbb{R}$, denote

[TABLE]

By Lemma 1 of Cule and Samworth [4], for any log-concave density $f_{0}$ there exist constants $\alpha_{f_{0}}>0$ and $\beta_{f_{0}}\in\mathbb{R}$ such that $f_{0}(x)\leq e^{\beta_{f_{0}}-\alpha_{f_{0}}|x|}$ for all $x\in\mathbb{R}$. Consequently, any upper semi-continuous log-concave density $f_{0}$ belongs to $\mathcal{F}_{\alpha,\beta}$ for $0<\alpha\leq\alpha_{f_{0}}$ and $\beta\geq\beta_{f_{0}}$.

We establish a general posterior contraction theorem for priors on log-concave densities using the general testing approach introduced in [12], which requires the construction of suitable tests with exponentially decaying type-II errors. We construct plug-in tests based on the concentration properties of the log-concave MLE, similar to the linear estimators considered in [15, 27]. The MLE has been shown to converge to the truth at the minimax rate in Hellinger distance in Kim and Samworth [23] and the following theorem relies heavily on their result.

Theorem 1**.**

Let $\mathcal{F}$ denote the set of upper semi-continuous, log-concave probability densities on $\mathbb{R}$ and let $\Pi_{n}$ be a sequence of priors supported on $\mathcal{F}$. Consider a sequence $\varepsilon_{n}\rightarrow 0$ such that $n^{-2/5}\lesssim\varepsilon_{n}\lesssim n^{-3/8-\rho}$ for some $\rho>0$ and suppose there exists a constant $C>0$ such that

[TABLE]

Then for sufficiently large $M$,

[TABLE]

in $P_{f_{0}}^{n}$-probability as $n\to\infty$.

The upper bound $\varepsilon_{n}\lesssim n^{-3/8-\rho}$ is an artefact of the proof arising from the exponential inequality for the log-concave MLE that we use to construct our tests, see Lemma 1. Since our interest lies in obtaining the optimal rate $n^{-2/5}$, possibly up to logarithmic factors, it plays no further role in our results. It is typical in Bayesian nonparametrics to require metric entropy conditions, which come from piecing together tests for Hellinger balls into tests for the complements of balls, see for instance Theorem 7.1 of [12]. The lack of such a condition in Theorem 1 is tied to the optimality and specific structure of the log-concave MLE. Using the affine equivariance of the MLE (Remark 2.4 of [10]), one can reduce the testing problem to considering alternatives in the class $\mathcal{F}$ restricted to have zero mean and unit variance. Unlike the whole space $\mathcal{F}$, the bracketing Hellinger entropy of this latter set can be suitably controlled, thereby avoiding the need for additional entropy bounds.

Remark 1**.**

Obtaining sharp entropy bounds for log-concave function classes is a highly technical task and such bounds are only available for certain restricted subsets. Even in the case of mean and variance restrictions (Theorem 4 of [23]) and compactly supported and bounded densities (Proposition 14 of [22]), the proofs are lengthy and require substantial effort. To use such bounds for the classic entropy-based approach to prove posterior contraction would therefore require the prior to place most of its mass on the above types of restricted sets. For instance, the prior might be required to place all but exponentially small probability on $\mathcal{F}_{\alpha,\beta}$ for some given $\alpha>0$, $\beta\in\mathbb{R}$. Such a prior construction is undesirable in practice and in fact none of our proposed priors satisfy such a restriction.

We now introduce a prior on log-concave densities based on an exponentiated Dirichlet process mixture. For any measurable function $w:\mathbb{R}\rightarrow\mathbb{R}$, define the density

[TABLE]

which is well-defined if $\int_{\mathbb{R}}e^{w(y)}dy<\infty$. Recall that any monotone non-increasing probability density on $\mathbb{R}^{+}$ has a mixture representation [40]

[TABLE]

where $P$ is a probability measure on $\mathbb{R}^{+}$. Khazaei and Rousseau [21] and Salomond [29] used the above representation to obtain a Bayesian nonparametric prior for monotone non-increasing densities. Unfortunately, such a convenient mixture representation is unavailable for log-concave densities and so the prior construction is somewhat more involved. Integrating the right-hand side of the last display, we obtain a function $w:\mathbb{R}^{+}\rightarrow\mathbb{R}$ as follows:

[TABLE]

where $\gamma_{1}>0$, $\gamma_{2}\in\mathbb{R}$ and $P$ is a probability measure on $[0,\infty)$. Since its (left and right) derivative is monotone decreasing, $w$ is concave. While not every concave function can be represented in this way, any log-concave density on $[0,\infty)$ can be approximated arbitrary well in Hellinger distance by a function of the form $e^{w}/(\int e^{w})$, where $w$ is as above with $P$ a discrete probability measure, see Proposition 1. Translating the above thus gives a natural representation for a prior construction for log-concave densities on $\mathbb{R}$.

Consider therefore the following possibly $n$-dependent prior on the log-density $w:[a_{n},b_{n}]\rightarrow\mathbb{R}$, where possibly $a_{n}\rightarrow-\infty$ and $b_{n}\rightarrow\infty$:

[TABLE]

where

• •

$P\sim DP(H\mathds{1}_{[0,b_{n}-a_{n}]})$, the Dirichlet process with base measure $H\mathds{1}_{[0,b_{n}-a_{n}]}=H(\mathbb{R}^{+})\bar{H}\mathds{1}_{[0,b_{n}-a_{n}]}$, where $0<H(\mathbb{R}^{+})<\infty$, $\bar{H}$ is a probability measure on $\mathbb{R}^{+}$ and every subset $U\subset[0,b_{n}-a_{n}]$ satisfies $H(U)\gtrsim\lambda(U)/(b_{n}-a_{n})^{\eta}$ for some $\eta\geq 0$,

• •

$\gamma_{i}\sim p_{\gamma_{i}}$, $i=1,2,$ where $p_{\gamma_{1}}$, $p_{\gamma_{2}}$ are probability densities on $[0,\infty)$ and $\mathbb{R}$ respectively, satisfying $p_{\gamma_{i}}(|x|)\gtrsim e^{-c_{i}x^{1/4}}$, $c_{i}>0$, for all $x\in[0,\infty)$ and $x\in\mathbb{R}$ respectively,

• •

$\gamma_{1}$, $\gamma_{2}$, and $P$ are independent.

We denote by $\Pi_{n}$ the full prior induced by $f_{W}$, where $W$ is drawn as above. Some typical draws from the prior are plotted in Figure 1.

Remark 2**.**

If $(b_{n}-a_{n})$ grows polynomially in $n$, then $H$ must have polynomial tails. On the other hand, if $(b_{n}-a_{n})$ grows more slowly than any polynomial, one can relax this condition. For instance, if $H$ has a density $h$ with respect to the Lebesgue measure, then it is sufficient that $\min_{t\in[0,b_{n}-a_{n}]}h(t)\gtrsim n^{-\lambda}$ for some $\lambda>0$. In particular, if $(b_{n}-a_{n})\lesssim\log n$, then $h$ may have exponential tails.

We comment on several aspects of our prior. Firstly, since Dirichlet process draws are atomic with probability one, the prior draws (3) will be piecewise linear and concave. Moreover, we could add any concave function to (3), such as an $-\gamma_{3}x^{2}$-type term, and still have a suitable concave prior. This permits greater modelling flexibility but complicates computation. In any case, the prior described above gives optimal contraction rates and can be computed in practice, so we restrict our attention to it. Another point to note is that if $(b_{n}-a_{n})\rightarrow\infty$ and $H$ is supported on the whole of $\mathbb{R}^{+}$, then the Dirichlet process base measure has total mass $H(\mathbb{R}^{+})\bar{H}([a_{n},b_{n}])\leq H(\mathbb{R}^{+})$ for fixed $n$. This has the interpretation of assigning the prior more weight as $n\rightarrow\infty$, up to the full prior weight $H(\mathbb{R}^{+})$. An alternative would be to re-weight the base measure to have full mass $H(\mathbb{R}^{+})$ to give it equal weight for all $n$. This plays no role asymptotically and so we restrict to the first case for technical convenience.

A potentially more serious issue is that for fixed $n$, the support of the prior draws may not contain the support of the true density $f_{0}$, in which case observations outside $[a_{n},b_{n}]$ cause the likelihood to be identically zero. While this is not a problem for $n$ large enough if $-a_{n},b_{n}\rightarrow\infty$ fast enough (see Theorem 3), it can be an issue for finite $n$. In practice, if one has an idea of the support of $f_{0}$, it is enough to select $[a_{n},b_{n}]$ large enough to contain $\text{supp}(f_{0})$. A more pragmatic solution is to use an empirical Bayes approach and make the prior data-dependent by setting $a_{n}:=X_{(1)}$, $b_{n}:=X_{(n)}$ the first and last order statistics. This ensures that the likelihood is never zero and the posterior is always well-defined. Indeed, the MLE is supported on $[X_{(1)},X_{(n)}]$ and so this can be thought of as plugging-in an estimate of the approximate support based on the likelihood. Moreover, since this approach yields the smallest support $[a_{n},b_{n}]$ with non-zero likelihood, it also brings computational advantages. In particular, it can prevent the need to simulate the posterior distribution on potentially very large regions of $\mathbb{R}$ where the posterior draws are essentially indistinguishable from zero. The empirical Bayes method behaves very similarly to the prior (3) in simulations and we would advocate this approach in practice.

We first present a contraction result when the true density $f_{0}$ has known compact support.

Theorem 2**.**

Let $f_{0}\in\mathcal{F}_{\alpha,\beta}$ for some $\alpha>0$, $\beta\in\mathbb{R}$ and suppose further that $f_{0}$ is compactly supported. Let $a_{n}\equiv a$ and $b_{n}\equiv b$ for all $n$ and denote by $\Pi_{n}=\Pi$ the prior described above. If $\text{supp}(f_{0})\subset[a,b]$, then

[TABLE]

in $P_{f_{0}}^{n}$-probability for some $M=M(\alpha,\beta)>0$.

If $\text{supp}(f_{0})$ is not contained in a compact set or is unknown, it suffices to let $-a_{n},b_{n}\rightarrow\infty$ fast enough. A slightly stronger lower bound on the tail of $p_{\gamma_{1}}$ is consequently required, depending on the size of $(b_{n}-a_{n})$.

Theorem 3**.**

Let $f_{0}\in\mathcal{F}_{\alpha,\beta}$ for some $\alpha>0$, $\beta\in\mathbb{R}$ and let $\Pi_{n}$ denote the prior described above with $-a_{n},b_{n}\gg\log n$. Assume further that $(b_{n}-a_{n})\lesssim n^{\mu/5}$ and that the prior density $p_{\gamma_{1}}$ for $\gamma_{1}$ satisfies the stronger lower bound $p_{\gamma_{1}}(x)\gtrsim e^{-c_{1}x^{1/(4+\mu)}}$ for some $0\leq\mu\leq 2$. Then

[TABLE]

in $P_{f_{0}}^{n}$-probability for some $M=M(\alpha,\beta)>0$ and

[TABLE]

Theorem 2 follows immediately from Theorem 3 and so its proof is omitted. If $(b_{n}-a_{n})=O((\log n)n^{2/5})$, then we obtain the minimax rate for log-concave density estimation in Theorem 3, up to a logarithmic factor. Since the Hellinger distance dominates the total variation distance, the above also implies posterior convergence in total variation at the same rate $\varepsilon_{n}$ given in (4). We also note that all the above statements are proved uniformly over $f_{0}\in\mathcal{F}_{\alpha,\beta}$.

The posterior mean, also considered in the simulation study, is not necessarily log-concave. Nevertheless one can construct log-concave density estimators by separately computing the posterior mean for each parameter $\theta,p,\gamma_{1},\gamma_{2}$ and then constructing the corresponding log-concave density according to (2) and (3). Another approach is to take the smallest Hellinger ball accumulating, say, $50\%$ of the posterior mass and sample an arbitrary log-concave density from that ball. It is straightforward to verify that both of these estimators achieve the minimax concentration rate (up to the same logarithmic factor).

It is also of interest to obtain a fully Bayesian procedure that does not require the user to define the support of the prior draws. We therefore consider a hierarchical prior where one places a prior on the end points $a$ and $b$, now not necessarily depending on $n$. This method has the advantage of employing a prior that does not depend on the data, but is slightly more involved computationally than the simple empirical Bayes approach. Assign to $(a,b)$ a prior supported on the open half-space $\{(a,b):a<b\}$ that has a Lebesgue density $\pi(a,b)$ satisfying

[TABLE]

and some $c_{1},c_{2},C,q,r>0$. Such a distribution can be easily constructed by first drawing $a\sim\pi_{1}$, where the Lebesgue density $\pi_{1}(a)\geq Ce^{-c|a|^{q}}$, and then independently drawing $(b-a)|a\sim\pi_{2}$, where $\pi_{2}$ is a Lebesgue density on $(0,\infty)$ satisfying $\pi_{2}(b-a)\geq Ce^{-c|b-a|^{r}}$. Conditionally on $(a_{n},b_{n})=(a,b)$, the prior is then exactly as above. This hierarchical construction leads to a fully Bayesian procedure that again contracts at the (near-)minimax rate.

Theorem 4**.**

Let $f_{0}\in\mathcal{F}_{\alpha,\beta}$ for some $\alpha>0$, $\beta\in\mathbb{R}$ and let $\Pi_{n}$ denote the prior described above with hyperprior $(a,b)\sim\pi(a,b)$ satisfying (5). Then

[TABLE]

in $P_{f_{0}}^{n}$-probability for some $M=M(\alpha,\beta)>0$.

Dirichlet process mixture priors are popular in density estimation due to the conjugacy of the posterior distribution, thereby providing methods that are highly efficient computationally. However, due to the exponentiation (2), this conjugacy property no longer holds, resulting in a less attractive prior choice that brings computational challenges. In practice, it is common to use approximations of the Dirichlet process to speed up computations, see for instance Chapter 4.3.3 of [14].

We firstly consider the Dirichlet multinomial distribution as a replacement for the Dirichlet process in our prior. By the proof of Theorem 3, the underlying true log-concave density can be well approximated by a piecewise log-linear density with at most $N=Cn^{1/5}\log n$ knots, for some large enough constant $C>0$. In view of this, it is reasonable to take $N$ atoms in the distribution. The corresponding prior on log-concave densities then takes the form

[TABLE]

where $\alpha_{i},$ $i=1,...,N,$ are chosen such that $\alpha_{i}=\alpha/N$ for some arbitrary $0<\alpha\leq H(\mathbb{R}^{+})$.

An alternative choice for the mixing prior is to truncate the stick-breaking representation of the Dirichlet process at a fixed level. Similarly to the Dirichlet multinomial distribution, we truncate the stick-breaking process at level $N=Cn^{1/5}\log n$, resulting in the same hierarchical prior as in (6) with the only difference being that the distribution of $p$ in the $N$-simplex is given by

[TABLE]

Both of these computationally more efficient approximations have the same theoretical guarantees as the full exponentiated Dirichlet process prior $\Pi_{n}$ or its hierarchical Bayes equivalent.

Corollary 1**.**

Let $f_{0}\in\mathcal{F}_{\alpha,\beta}$ for some $\alpha>0$, $\beta\in\mathbb{R}$ and let $\Pi_{n}^{\prime}$ denote either the prior (6) or (7). If $-a_{n},b_{n}\gg\log n$, $(b_{n}-a_{n})\lesssim n^{\mu/5}$ and the prior density $p_{\gamma_{1}}$ for $\gamma_{1}$ satisfies the stronger lower bound $p_{\gamma_{1}}(x)\gtrsim e^{-c_{1}x^{1/(4+\mu)}}$ for some $0\leq\mu\leq 2$, then

[TABLE]

in $P_{f_{0}}^{n}$-probability for some $M=M(\alpha,\beta)>0$ and $\varepsilon_{n}$ given by (4).

If we additionally assign a hyperprior $\pi(a,b)$ satisfying (5) to $(a,b)$ and no longer require the strong lower bound for $p_{\gamma_{1}}$, then the above holds for the resulting posterior with $\varepsilon_{n}=(\log n)n^{-2/5}$.

The proofs of Theorem 3 and Corollary 1 establish the small-ball probability (1) by approximating a log-concave density in $\mathcal{F}_{\alpha,\beta}$ with a suitable piecewise log-linear density. This approximation requires several key properties, which make its construction non-standard and technically involved, and it may be of independent interest. The proof of Proposition 1 is deferred to Section 6.1.

Proposition 1**.**

Let $f_{0}\in\mathcal{F}_{\alpha,\beta}$ and $([a_{n},b_{n}])_{n}$ be a sequence of compact intervals such that $[-\tfrac{8}{5\alpha}\log n,\tfrac{8}{5\alpha}\log n]\subset[a_{n},b_{n}]$ and $(b_{n}-a_{n})=o(n^{4/5})$. For any $n\geq n_{0}$, where $n_{0}$ is an integer depending only on $\alpha$ and $\beta$, there exists a log-concave density $\bar{f}_{n}$ that is piecewise log-linear with $\bar{N}\leq C(\alpha,\beta)n^{1/5}\log n$ knots $z_{1},...,z_{\bar{N}}\in[0,b_{n}-a_{n}]$ satisfying the following properties:

• (i)

$h^{2}(f_{0},\bar{f}_{n})\leq C(\alpha,\beta)[(\log n)^{2}n^{-4/5}+(b_{n}-a_{n})^{2}n^{-8/5}]$,

• (ii)

$\{x\in\mathbb{R}:\bar{f}_{n}(x)>0\}=[a_{n},b_{n}]$,

• (iii)

the knots are $cn^{-6/5}\log n$-separated for some universal constant $c>0$,

• (iv)

$f_{0}(x)\leq C(\alpha,\beta)\bar{f}_{n}(x)$* for all $x\in[a_{n},b_{n}]$,*

• (v)

there exist $\bar{\gamma}_{1}\in[0,2(b_{n}-a_{n})n^{4/5}]$, $|\bar{\gamma}_{2}|\leq n^{4/5}$, $\bar{\gamma}_{3}\in\mathbb{R}$ and $(\bar{p}_{1},..,\bar{p}_{\bar{N}})$ satisfying $p_{i}\geq 0$ and $\sum_{i=1}^{\bar{N}}p_{i}=1$, such that

[TABLE]

It is relatively straightforward to establish an approximation of $f_{0}$ satisfying $(i)$. However, approximating $f_{0}$ by $\bar{f}_{n}$ in a Kullback-Leibler type sense, as in (1), necessitates control of the support of $\bar{f}_{n}$ via $(ii)$ and uniform control of the ratio $f_{0}/\bar{f}_{n}$ via $(iv)$. The most difficult property to establish is the polynomial separation of the points in $(iii)$. This is needed to ensure that the Dirichlet process prior simultaneously puts sufficient mass in a neighbourhood of each of the knots $z_{i}$, $i=1,\dots,\bar{N}$. Setting $[a_{n},b_{n}]=[-\tfrac{8}{5\alpha}\log n,\tfrac{8}{5\alpha}\log n]$ yields the following corollary.

Corollary 2**.**

Let $f_{0}\in\mathcal{F}_{\alpha,\beta}$. For any $n\geq n_{0}$, where $n_{0}$ is an integer depending only on $\alpha$ and $\beta$, there exists a log-concave density $\bar{f}_{n}$ supported on $[-\tfrac{8}{5\alpha}\log n,\tfrac{8}{5\alpha}\log n]$ that is piecewise log-linear with $O(n^{1/5}\log n)$ knots and satisfies $h^{2}(f_{0},\bar{f}_{n})\leq C(\alpha,\beta)(\log n)^{2}n^{-4/5}$. Moreover, we may take the knots to be $cn^{-6/5}\log n$-separated for some universal constant $c>0$.

3 Simulation study

We present a simulation study to assess the performance of the proposed log-concave priors for density estimation. In particular, we investigate the prior based on the truncated stick breaking representation (7), firstly with deterministically chosen support $[a_{n},b_{n}]$, secondly its empirical Bayes counterpart with support $[X_{(1)},X_{(n)}]$, where $X_{(1)}$ and $X_{(n)}$ denote the smallest and largest observations, respectively, and thirdly the hierarchical Bayes version where the parameters $a$ and $b-a$ are endowed with independent Cauchy and half-Cauchy distributions, respectively. In all cases we plot the posterior mean and $95\%$ pointwise credible sets, and compare them with the log-concave maximum likelihood estimator (computed using the R function “mlelcd”).

Consider first the posterior distribution arising from the prior with deterministic support $[a_{n},b_{n}]$. We have drawn random samples of size $n=50,200,500$ and $2500$ from a gamma distribution with shape and rate parameters 2 and 1, respectively. We took the number of linear pieces in the exponent of the prior to be $m=n^{1/5}\log n$, set $[a_{n},b_{n}]=[-2.3\log n,2.3\log n]$, endowed the break-point parameters $\theta=(\theta_{1},...,\theta_{m})$ with independent uniform priors on $[0,b_{n}-a_{n}]$, assigned the weight parameters $p$ a stick-breaking distribution truncated at level $m$, and endowed $\gamma_{1}$ and $\gamma_{2}$ a half Cauchy and a Cauchy distribution, respectively, with location parameter 0 and scale parameter 1. Since the posterior distribution does not have a closed-form expression, we drew approximate samples from the posterior using a random walk Metropolis-Hastings within Gibbs sampling algorithm for 10000 iterations out of which the first 5000 are discarded as burn-in period. In Figure 2, we have plotted the true distribution (solid red), posterior mean (solid blue), $95\%$ pointwise credible band (dashed blue) and the maximum likelihood estimator (solid green). The data is represented by a histogram on the figures.

We see that the posterior mean gives an adequate estimator for the true log-concave density with similar, if not superior, performance compared to the more jagged maximum likelihood estimator, and the $95\%$ pointwise credible bands mostly contain the true function except for points close to zero. We further investigate the frequentist coverage properties of the pointwise Bayesian credible sets. We repeat the above experiment for the empirical Bayes procedure 100 times (each with 2000 iterations out of which half were discarded as burn in) and report the frequencies where the density $f(x)$ at given points $x\in\{0.5,1,1.5,2,2.5,3\}$ is inside of the corresponding credible interval. We consider sample sizes $n=50,200$ and $500$ and report the empirical coverage probabilities in Table 1. One can see that in this particular example we get quite reliable uncertainty quantification, especially for larger sample sizes. It should be noted, however, that the frequentist coverage of Bayesian credible sets is a delicate subject in nonparametric statistics, see for instance Szabó et al. [34], and is beyond the scope of this article.

We next investigate the behaviour of the empirical and hierarchical Bayes versions of the proposed prior. We again simulate $n=50,200,500$ and $2500$ independent draws from a Gamma(2,1) distribution and set the compact support of the prior densities to be $[a_{n},b_{n}]=[X_{(1)},X_{(n)}]$, that is the smallest and largest observations, for the empirical Bayes procedure and $[a,b]$, with $a\sim Cauchy(0,1)$ and $b-a\sim Cauchy_{+}(0,1)$ independent, for the hierarchical Bayes method. As before we set $m=n^{1/5}\log n$ and endowed the parameters $\theta,p,\gamma_{1}$ and $\gamma_{2}$ with the same priors as above. We ran the algorithm again for 10000 iterations, taking the first half of the chain as a burn-in period. We plot the outcomes in Figures 3 and 4 for the empirical and hierarchical Bayes procedures, respectively. One can see that for $n\geq 500$ observations, the posterior mean (solid blue) closely resembles the underlying gamma density (solid red), while the fit is already reasonable for $n=200$. The pointwise $95\%$-credible bands contain the true density, even near zero, which was problematic in case of the prior with support selected deterministically. Comparing Figures 2, 3 and 4, we see that the empirical and hierarchical Bayes approaches of selecting the support $[a_{n},b_{n}]$ in a data-driven way outperform a deterministic selection. We also note that the algorithm for the empirical Bayes method was considerably faster than the others due to the smaller support, which reduces the computation time of the normalizing constants $\int e^{w(y)}dy$ of the densities.

We then investigate the performance of the posterior distribution corresponding to the empirical and hierarchical Bayes methods for recovering different log-concave densities and again compare them with the MLE. We have considered a standard normal distribution, a gamma distribution with shape parameter 2 and rate parameter 1, a beta distribution with shape parameters 2 and 3, and a Laplace distribution with location parameter 0 and dispersion parameter 1. In all four examples we have taken sample size $n=1500$. The posterior mean (solid blue), the $95\%$ pointwise credible bands (dashed blue) and the MLE (green) are plotted in Figures 5 and 6 for the empirical and hierarchical Bayes procedures, respectively. All four subpictures for both data-driven methods show satisfactory results, both for recovery using the posterior mean and for uncertainty quantification using the pointwise credible bands. We note that the displayed plots convey typical behaviour and are representative of multiple simulations. We hence draw the conclusion that the proposed method seems to work well in practice for various choices of common log-concave densities.

Lastly, we investigate the performance of the proposed Bayesian procedures for estimating the mode of the underlying log-concave density. We consider the standard normal distribution and take i.i.d. random samples of size ranging from 50 to 20000. We run the Gibbs sampler for 20000 iterations and take the first half of the iterations as burn-in period. For each posterior draw we compute the mode and use the resulting histogram to approximate the one-dimensional marginal posterior. The histograms from the empirical Bayes procedure are displayed in Figure 7. One can see that the posterior concentrates around the true mode (i.e. 0) as the sample size increases.

The marginal posterior concentrates substantially slower than $n^{-1/2}$-rate. This is as expected, since the best possible minimax rate for estimating the mode $m_{0}$ of a unimodal or log-concave density $f_{0}$ satisfying $f_{0}^{\prime\prime}(m_{0})<0$ is $n^{-1/5}$, see [19, 2]. Indeed, the mode of the log-concave MLE attains this rate [2]. Interestingly, the marginal posterior does not seem to be Gaussian, which may be linked to the irregular asymptotic distribution of the mode of the log-concave MLE. This rather complicated distribution equals the mode of the second derivative of the lower invelope of a certain Gaussian process, see [2] for full details. A better understanding of the limiting shape of the marginal posterior would be interesting, but is beyond the scope of this article.

In the supplementary material we provide additional simulations for the marginal posterior for the mode from the empirical Bayes posterior for different underlying log-concave densities, namely the beta and gamma distributions. We also numerically investigate the applicability of our log-concave Bayesian prior for estimating mixtures of log-concave densities.

4 Discussion

We have proposed a novel Bayesian procedure for log-concave density estimation. The prior is defined on compactly supported densities, where the support can be chosen either deterministically, empirically or using a fully Bayesian hierarchical procedure. We have shown theoretically that both the deterministic and fully Bayesian choices of the support give (near-)optimal posterior contraction rates, and have demonstrated the good small sample performance of the posterior for all three methods in a simulation study. We have also plotted the $95\%$ pointwise credible bands, which in our simulation study provide reliable frequentist uncertainty quantification. However, this might depend heavily on the choice of the underlying true density and it is unclear at present whether our methods generally provide trustworthy frequentist uncertainty quantification. The rigourous study of this question is beyond the scope of the present paper.

In our simulation study, we further investigated the behaviour of the marginal posterior for the mode functional. A natural next question is whether one can obtain semiparametric Bernstein-von Mises type results for the mode. In view of the irregular behaviour of the log-concave MLE, this is an interesting problem as it is unclear whether the limiting distribution of the posterior is indeed Gaussian.

A possible application of our proposed approach is clustering based on mixture models. Assuming that clusters have log-concave densities instead of (say) Gaussians broadens their modelling flexibility. We have executed a small simulation study to explore this direction. For simplicity we have considered a mixture of only two log-concave densities and modified our prior accordingly. In the considered examples (see the Supplementary material) our procedure performs reasonably well. However, we should note that the computational time is much worse than using simple Gaussian kernels. Extending this to mixtures with more than two (possibly unknown number of) components seems to be possible, but requires optimization of the Gibbs sampler and perhaps introducing other approximation steps, which are beyond the scope of the present paper.

Another natural question if whether one can extend these results to multivariate density estimation, especially in view of the difficulty of computing the log-concave MLE in higher dimensions. Since our present prior construction is based on using a mixture representation to model a decreasing function, which corresponds to the derivative of the concave exponent, this will require new ideas. A possible approach is presented in Hannah and Dunson [18], who place a prior over all functions that are the maximum of a set of hyperplanes. This yields a prior on the set of convex functions that could potentially be adapted to the multivariate log-concave setting.

5 Proofs

Define the following classes of log-concave densities with mean and variance restrictions:

[TABLE]

and

[TABLE]

Let $\hat{f}_{n}$ denote the log-concave MLE based on i.i.d. random variables $X_{1},\dots,X_{n}$ arising from a density $f_{0}\in\mathcal{F}$.

The proof of Theorem 1 relies on a concentration inequality for the log-concave MLE based on data from moment-restricted densities. This is the content of the following lemma, whose proof is essentially contained in Kim and Samworth [23] for the more difficult case of general $d\geq 1$. However, we require a sharper probability bound than they provide and so make some minor modifications to their argument. The proof is deferred to Section 6.2.

Lemma 1**.**

For every $\varepsilon>0$, there exist positive constants $L_{0},C,c,n_{0}$, depending only on $\varepsilon$, and positive universal constants $D,d>0$, such that for all $L\geq L_{0}$ and $n\geq n_{0}$,

[TABLE]

where $\hat{g}_{n}$ denotes the log-concave maximum likelihood estimator based on an i.i.d. sample $Z_{1},...,Z_{n}$ from $g_{0}$.

Proof of Theorem 1.

As in the proof of Theorem 2.1 of [12], using the lower bound on the small ball probability from (1), it suffices to construct tests $\phi_{n}=\phi_{n}(X_{1},...,X_{n};f_{0})$ such that

[TABLE]

for $n$ large enough, where the constant $C>0$ matches that in (1).

For $M_{0}$ a constant to be chosen below, set $\phi_{n}=\mathds{1}\{h(\hat{f}_{n},f_{0})\geq M_{0}\varepsilon_{n}\}$, where $\hat{f}_{n}$ is the log-concave MLE based on i.i.d. observations $X_{1},...,X_{n}$ from a density $f_{0}\in\mathcal{F}$. Let $\mu_{f_{0}}=\mathbb{E}X_{i}$, $\sigma_{f_{0}}^{2}=\textnormal{Var}(X_{i})$ and define $Z_{i}=(X_{i}-\mu_{f_{0}})/\sigma_{f_{0}}$, so that $\mathbb{E}Z_{i}=0$ and $\textnormal{Var}(Z_{i})=1$. Further set $g_{0}(z)=\sigma_{f_{0}}f_{0}(\sigma_{f_{0}}z+\mu_{f_{0}})$ and $\hat{g}_{n}(z)=\sigma_{f_{0}}\hat{f}_{n}(\sigma_{f_{0}}z+\mu_{f_{0}})$, so $g_{0}\in\overline{\mathcal{F}}^{0,1}$. By affine equivariance (Remark 2.4 of [10]), $\hat{g}_{n}$ is the log-concave maximum likelihood estimator of $g_{0}$ based on $Z_{1},...,Z_{n}$.

Using the invariance of the Hellinger distance under affine transformations and Lemma 1 with $\varepsilon=1/2$, the type-I error satisfies

[TABLE]

as $n\rightarrow\infty$ for $M_{0}$ large enough since $\varepsilon_{n}\gtrsim n^{-2/5}$. For $f\in\mathcal{F}$ such that $h(f,f_{0})\geq M\varepsilon_{n}$,

[TABLE]

Since $\varepsilon_{n}\lesssim n^{-3/8-\rho}$ implies $\varepsilon_{n}\lesssim n^{-3/8-\rho^{\prime}}$ for any $0<\rho^{\prime}\leq\rho$, we may take $\rho>0$ arbitrarily small. Applying Lemma 1 with $\varepsilon(\rho)>0$ to be chosen below and $L_{n}=(M-M_{0})\varepsilon_{n}n^{2/5}$, which satisfies $L_{n}\geq L_{0}$ for $M>0$ large enough since $\varepsilon_{n}\gtrsim n^{-2/5}$, yields

[TABLE]

for all $n\geq n_{0}(\rho)$. Since $\varepsilon_{n}\lesssim n^{-3/8-\rho}$ by assumption, it follows that $n\varepsilon_{n}^{2}\lesssim n^{1/4-2\rho}=o(n^{1/(4+2\varepsilon(\rho))})$ for $\varepsilon(\rho)>0$ small enough. Therefore,

[TABLE]

Since we can make the constant in the exponent arbitrarily large by taking $M>0$ large enough, this completes the proof. ∎

Proof of Theorem 3.

Let $f_{0}\in\mathcal{F}_{\alpha,\beta}$ for some $\alpha>0$ and $\beta\in\mathbb{R}$. We may restrict to a suitable compactly supported density approximating $f_{0}$ using the first paragraph of the proof of Theorem 2 of Ghosal and van der Vaart [13]. For completeness we reproduce their argument in this paragraph. Let $\psi_{n}(x)=\mathds{1}_{[-t_{n},t_{n}]}(x)$ for $t_{n}=a^{\prime}\log n$ for some $a^{\prime}>\alpha^{-1}$. Define new observations $\bar{X}_{1},...,\bar{X}_{\bar{n}}$ from the original observations $X_{1},...,X_{n}$ by rejecting each $X_{i}$ independently with probability $1-\psi_{n}(X_{i})$. Since $P_{f_{0}}[-t_{n},t_{n}]^{c}\leq 2e^{\beta}\alpha^{-1}e^{-\alpha t_{n}}=o(n^{-1})$, the probability that at least one of the $X_{i}$’s is rejected is $o(1)$ and so the posterior based on the original and modified observations are the same with $P_{f_{0}}^{n}$-probability tending to one. Since posterior contraction is defined via convergence in $P_{f_{0}}^{n}$-probability, this implies that the posterior contraction rates are the same. The new observations come from a density $f_{0,n}$ that is proportional to $f_{0}\psi_{n}$, which is log-concave and upper semi-continuous. Since $|1-\int f_{0}\psi_{n}|\leq P_{f_{0}}[-t_{n},t_{n}]^{c}=o(n^{-1})$,

[TABLE]

It therefore suffices to establish contraction for the posterior based on the new observations about the density $f_{0,n}=f_{0}\psi_{n}/\int f_{0}\psi_{n}$.

Under the assumed conditions on $(b_{n}-a_{n})$, $\varepsilon_{n}$ given by (4) satisfies $n^{-2/5}\lesssim\varepsilon_{n}\lesssim n^{-3/8-\rho}$ for some $\rho>0$ small enough. We thus apply Theorem 1 so that we need only show the small-ball probability (1). Note that $f_{0,n}(x)\leq e^{\beta-\alpha|x|}(1+o(n^{-1}))$, so that $f_{0,n}\in\mathcal{F}_{\alpha,2\beta}$ for $n$ large enough. Since $-a_{n},b_{n}\gg\log n$ and $b_{n}-a_{n}=o(n^{4/5})$, we may construct an approximation $\bar{f}_{n}$ of $f_{0,n}$ based on the interval $[a_{n},b_{n}]$ for $n$ large enough using Proposition 1. By Lemma 8 of [13],

[TABLE]

for $k=1,2$. By Proposition 1(i) and (iv), the first term in the first bracket and the second term in the second bracket are $O((\log n)^{2}n^{-4/5}+(b_{n}-a_{n})^{2}n^{-8/5})$ and $O(1)$ respectively.

By Proposition 1(v), $\bar{f}_{n}$ has representation

[TABLE]

where $(z_{i})_{i=1}^{\bar{N}}\subset[0,b_{n}-a_{n}]$ are the knots written in increasing order, $\bar{N}=\bar{N}_{n}=O(n^{1/5}\log n)$ and $\sum_{i=1}^{\bar{N}}\bar{p}_{i}=1$. Let $\bar{w}_{n}(x)=(\log\bar{f}_{n}(x)-\bar{\gamma}_{3})\mathds{1}_{[a_{n},b_{n}]}(x)-\infty\mathds{1}_{\mathbb{R}\backslash[a_{n},b_{n}]}(x)$ so that $\bar{f}_{n}=f_{\bar{w}_{n}}$ using the transformation (2). We may thus without loss of generality take $\bar{\gamma}_{3}=0$ since it is contained in the normalization (2).

Suppose that $f_{w}$ is a (log-concave) density with support equal to $[a_{n},b_{n}]$ and such that $\|\bar{w}_{n}-w\|_{L^{\infty}([a_{n},b_{n}])}\leq c\varepsilon_{n}$. Since $\int e^{w}=e^{O(\varepsilon_{n})}\int e^{\bar{w}_{n}}$, it follows that for $x\in[a_{n},b_{n}]$, $\bar{f}_{n}(x)/f_{w}(x)\leq e^{O(\varepsilon_{n})}=e^{o(1)}$. Since $h^{2}(\bar{f}_{n},f_{w})\lesssim\varepsilon_{n}^{2}$ by Lemma 3.1 of [36], we can conclude that

[TABLE]

for some $c>0$. It therefore suffices to lower bound the prior probability of the left-hand set.

Fix $\delta>0$ to be chosen sufficiently large below. Since the $z_{i}$ are $n^{-6/5}$-separated by Proposition 1(iii), we can find a collection of disjoint intervals $(U_{i})_{i=1}^{\bar{N}}$ in $[a_{n},b_{n}]$ with Lebesgue measure $\lambda(U_{i})=\varepsilon_{n}^{\delta}$ and such that $z_{i}\in U_{i}$ for $i=1,...,\bar{N}$. Further denote $U_{0}:=\mathbb{R}\backslash\cup_{i=1}^{M_{n}}U_{i}$. Let $W$ be a prior draw of the form (3) with parameters $\gamma_{1},\gamma_{2}$ and $P$. Writing $p_{i}=P(U_{i})$, $\bar{p}_{0}=0$ and using the triangle inequality, for any $x\in[a_{n},b_{n}]$,

[TABLE]

where we have used in the second to last line that the maximal distance between the (piecewise) lines $(y\wedge a)/a$ and $(y\wedge b)/b$ occurs at $y=a\wedge b$ and in the last line that $\theta>z_{i}/2$ for all $\theta\in U_{i}$ for a sufficiently large choice of the parameter $\delta>0$. By Proposition 1(v), we have $\bar{\gamma}_{1}\leq 2n^{4/5}(b_{n}-a_{n})\lesssim n^{(4+\mu)/5}$. Furthermore, by the separation of the knots, $z_{i}\geq z_{1}\geq cn^{-6/5}$, $i=1,...,\bar{N}$, and so by the assumptions on the $(U_{i})$, the second term is bounded by $2c^{-1}\bar{\gamma}_{1}n^{6/5}\varepsilon_{n}^{\delta}\leq\tilde{c}\varepsilon_{n}$ for some $\delta,\tilde{c}>0$ large enough.

The remaining three terms are independent under the prior and so can be dealt with separately. By the assumptions on the base measure of the Dirichlet process, we have that $\sum_{i=0}^{\bar{N}}H(U_{i})\leq H(\mathbb{R}^{+})$ and $H(U_{i})=H(\mathbb{R}^{+})\bar{H}(U_{i})\gtrsim\lambda(U_{i})/(b_{n}-a_{n})^{\eta}\geq\varepsilon_{n}^{\delta}/(b_{n}-a_{n})^{\eta}\geq\varepsilon_{n}^{\delta^{\prime}}$ for $i=1,...,\bar{N}$ and some $\delta^{\prime}>\delta$. For $i=0$, note that $\lambda(U_{0})\geq(b_{n}-a_{n})-\bar{N}\varepsilon_{n}^{\delta}\gtrsim 1$. Using the lower the bounds for the $\lambda(U_{i})$, which come from the polynomial separation of the knots in Proposition 1(iii), we can apply Lemma 10 of [13] to get

[TABLE]

From the tail assumption on the density of $\gamma_{1}$ and the upper bound on $\bar{\gamma}_{1}$, we have

[TABLE]

By Proposition 1(v), $|\bar{\gamma}_{2}|\lesssim n^{4/5}$, which, combined with the tail bound on the density of $\gamma_{2}$, yields

[TABLE]

since $\varepsilon_{n}/(b_{n}-a_{n})\rightarrow 0$ no faster than polynomially in $n$. Combining the above, we have that $\Pi_{n}(\|\bar{w}_{n}-W\|_{L^{\infty}([a_{n},b_{n}])}\leq(3+\tilde{c})\varepsilon_{n})\geq e^{-(2+c^{\prime})n\varepsilon_{n}^{2}}.$ ∎

Proof of Theorem 4.

Since the proof follows that of Theorem 3, we only specify the details where the present proof differs. Using the same arguments, we restrict to studying posterior contraction based on observations arising from the log-concave density $f_{0,n}=f_{0}\psi_{n}/\int f_{0}\psi_{n}$ for $\psi_{n}(x)=\mathds{1}_{[-t_{n},t_{n}]}(x)$ with $t_{n}=a^{\prime}\log n$ for some $a^{\prime}>\alpha^{-1}$. We again apply Theorem 1 so that we only need to show the small-ball probability (1) for $\varepsilon_{n}=(\log n)n^{-2/5}$.

Writing $\Pi_{a,b}$ for the prior conditional on $(a,b)$ and setting $\Delta_{n}=\{(a,b):-3t_{n}\leq a\leq-2t_{n},2t_{n}\leq b\leq 3t_{n}\}$, the small ball probability in (1) is lower bounded by

[TABLE]

Using the lower bound assumption (5) on $\pi(a,b)$, the last integral is lower bounded by $Ce^{-c_{1}(3t_{n})^{q}-c_{2}(6t_{n})^{r}}\int_{\Delta_{n}}da\,db\geq Ct_{n}^{2}e^{-c_{3}(\log n)^{q\vee r}}\gtrsim e^{-C^{\prime}n\varepsilon_{n}^{2}}$. It thus suffices to lower bound the infimum in the last display.

Since $[-\tfrac{8}{5\alpha}\log n,\tfrac{8}{5\alpha}\log n]\subset[a,b]$ and $(b-a)=O(\log n)=o(n^{4/5})$ for all $(a,b)\in\Delta_{n}$, we may apply Proposition 1 to construct an approximation $\bar{f}_{n}$ of $f_{0,n}$ based on the interval $[a,b]$ for any $(a,b)\in\Delta_{n}$. One can then proceed exactly as in the proof of Theorem 3 to lower bound the prior small-ball probability by $e^{-Cn\varepsilon_{n}^{2}}$ for fixed $(a,b)\in\Delta_{n}$. Since all constants in that argument depend only on $\alpha$, $\beta$ and the prior hyperparameters, the lower bound is uniform over all $(a,b)$ with $(b-a)=O(n^{2/5})$ (to ensure $\varepsilon_{n}$ in (4) takes the value $(\log n)n^{-2/5}$) and $-a,b\geq\tfrac{8}{5\alpha}\log n$. In particular, the lower bound is uniform over $\Delta_{n}$. ∎

Proof of Corollary 1.

We use the notation employed in the proof of Theorem 3. By (10), it suffices to lower bound the prior probability of an $L^{\infty}$-small ball about $\bar{w}_{n}$, where $\bar{f}_{n}=f_{\bar{w}_{n}}$ is the approximation (9). Since $\bar{N}\leq N$ (at least for $n$ large enough), we can add additional breakpoints to the piecewise linear function $\bar{w}_{n}$ with weights $\bar{p}_{i}=0$, $i=\bar{N}+1,...,N$, without changing $\bar{w}_{n}$. Without loss of generality, pick any such additional breakpoints to be no smaller than $cn^{-6/5}$. Using similar computations to (11), for any $x\in[a_{n},b_{n}]$,

[TABLE]

The first and fourth terms are bounded from above by $\varepsilon_{n}$ with prior probability at least $e^{-n\varepsilon_{n}^{2}}$ by (13) and (14), respectively, for both priors. For the second term note, similarly to the proof of Theorem 3, that $z_{1}\geq cn^{-6/5}$. Taking $\theta_{i}\in[z_{i},z_{i}+cn^{-6/5}\varepsilon_{n}/\bar{\gamma}_{1}]$, the second term is bounded by $\bar{\gamma}_{1}\sum_{i=1}^{N}p_{i}\varepsilon_{n}/\bar{\gamma}_{1}=\varepsilon_{n}$. The probability of this set under the base measure is H([z_{i},z_{i}+cn^{-6/5}\varepsilon_{n}/\bar{\gamma}_{1}])\gtrsim cn^{-6/5}\varepsilon_{n}/\big{(}\bar{\gamma}_{1}(b_{n}-a_{n})^{\eta}\big{)} by the assumptions on $\bar{H}$. The joint probability that $\theta_{i}\in[z_{i},z_{i}+cn^{-6/5}\varepsilon_{n}/\bar{\gamma}_{1}]$ for every $i=1,...,N$ is therefore bounded from below by a multiple of $(cn^{-6/5}\varepsilon_{n}/\bar{\gamma}_{1})^{N}/(b_{n}-a_{n})^{\eta N}\gtrsim e^{-c_{1}N\log n}\geq e^{-c_{2}n\varepsilon_{n}^{2}}$, for sufficiently large constants $c_{1},c_{2}>0$.

It remains to show that the third term is bounded from above by $\varepsilon_{n}$ with probability at least $e^{-cn\varepsilon_{n}^{2}}$ for some $c>0$. In the case where $(p_{1},...,p_{N})$ is endowed with a Dirichlet distribution, this statement follows from (12). In the case of the truncated stick-breaking prior, writing $(\bar{p}_{1},\dots,\bar{p}_{N})$ in decreasing order, we note that there exist $0\leq\bar{v}_{1},...,\bar{v}_{N-1}\leq 1$ such that $\bar{p}_{i}=\prod_{j=1}^{i-1}(1-\bar{v}_{j})\bar{v_{i}}$, $i=1,...,N-1$, and $\bar{p}_{N}=\prod_{j=1}^{N-1}(1-\bar{v}_{j})$. Define for $i=1,2,...,N-1$ the intervals

[TABLE]

For $v_{i}\in I_{i}\subset[0,1]$, $i=1,..,N-1$, we have

[TABLE]

Hence for $p_{i}:=(1-v_{1})(1-v_{2})...(1-v_{i-1})v_{i}$, $i=1,...,N-1$, $(p_{1},..,p_{N})$ is in the $N$-dimensional simplex and $\bar{\gamma}_{1}\sum_{i=1}^{N}|\bar{p}_{i}-p_{i}|\leq\bar{\gamma}_{1}N\varepsilon_{n}n^{-4/5}N^{-1}\lesssim\varepsilon_{n}$. Finally, we note that for $v_{i}\sim\text{Beta}(a,b)$, we have $P(v_{i}\in I_{i})\gtrsim(\varepsilon_{n}n^{-4/5}N^{-2})^{a\vee b}$ and we can therefore conclude

[TABLE]

for some large enough constants $c_{1},c_{2}>0$, thereby completing the proof.

For the hierarchical case where we assign a prior to $(a,b)$, the proof follows as in that of Theorem 4 using (15) and the lower bound for the small-ball probability just derived. ∎

6 Technical results

6.1 Proof of Proposition 1

In this section, we construct the piecewise log-linear approximation for an upper semi-continuous log-concave density given in Proposition 1. In particular, we require that the number of knots in the approximating function does not grow too quickly and that the knots are polynomially separated, thereby rendering the construction somewhat involved. The proof relies on firstly approximating any continuous concave function on a given compact interval using a piecewise linear function. One then splits $\text{supp}(f_{0})$ into sets, depending on the size of both $\log f_{0}$ and $|(\log f_{0})^{\prime}|$, and obtains suitable piecewise linear approximations defined locally on each of these sets. Piecing together these local functions gives the desired global approximation.

We now construct a piecewise linear approximation of a continuous concave function $w$ on a compact interval $[a,b]$. For any partition $a=x_{0}<x_{1}<\dots<x_{m}=b$ of $[a,b]$, let $\widetilde{w}_{m}$ denote the piecewise linear approximation of $w$ given by

[TABLE]

where $\theta_{i}:=\int_{x_{i-1}}^{x_{i}}w(s)ds$ and $x_{i}^{*}:=\frac{x_{i}+x_{i-1}}{2}$. On $[a,x_{1}^{*}]$ and $(x_{m}^{*},b]$, the function is defined by linearly extending the piecewise linear function defined above, that is

[TABLE]

The function $\widetilde{w}_{m}$ takes value $\widetilde{w}_{m}(x_{i}^{*})=\frac{1}{x_{i}-x_{i-1}}\int_{x_{i-1}}^{x_{i}}w(s)ds$ at the midpoint $x_{i}^{*}=\frac{x_{i}+x_{i-1}}{2}$ of the interval $[x_{i-1},x_{i}]$ and interpolates linearly in between.

Lemma 2**.**

Let $w:[a,b]\to\mathbb{R}$ be a continuous concave function, where $-\infty<a<b<\infty$. For any partition $a=x_{0}<x_{1}<\dots<x_{m}=b$ of $[a,b]$, let $\widetilde{w}_{m}$ denote the piecewise linear approximation of $w$ defined in (16) and (17). Then $\widetilde{w}_{m}$ is a concave function.

Proof.

By rescaling, we may without loss of generality assume that $[a,b]=[0,1]$. Note that $\widetilde{w}_{m}$ is concave if and only if

[TABLE]

for $i=2,\dots,m-1$. Indeed, since $\widetilde{w}_{m}$ is piecewise linear, it is concave if and only if at every point where the derivative is discontinuous (i.e. a knot), the left derivative is greater than or equal to the right derivative. The above statement follows since the derivative of $\widetilde{w}_{m}$ is discontinuous (at most) at the points $x_{i}^{*}$, $i=2,\dots,m-1$, where the desired inequality is:

[TABLE]

which is equivalent to (18).

To see that (18) holds if $w$ is concave, we argue by contradiction and suppose that there exists $i$ such that $\widetilde{w}_{m}(x_{i}^{*})<\frac{x_{i}^{*}-x_{i-1}^{*}}{x_{i+1}^{*}-x_{i-1}^{*}}\widetilde{w}_{m}(x_{i+1}^{*})+\frac{x_{i+1}^{*}-x_{i}^{*}}{x_{i+1}^{*}-x_{i-1}^{*}}\widetilde{w}_{m}(x_{i-1}^{*})$. Consider the linear function $l$,

[TABLE]

In particular, we have that $\widetilde{w}_{m}(x_{i-1}^{*})-l(x_{i-1}^{*})=\widetilde{w}_{m}(x_{i+1}^{*})-l(x_{i+1}^{*})=0$ and $\widetilde{w}_{m}(x_{i}^{*})<l(x_{i}^{*})$. We further denote $g:=w-l$ and observe that

[TABLE]

It follows that $\tilde{g}_{m}(x)=\widetilde{w}_{m}(x)-l(x)$ for all $x\in[0,1]$ and hence by the mean value theorem,

[TABLE]

for some $\xi_{i}\in[x_{i-1},x_{i}]$. One can similarly prove the existence of two points, $\xi_{i-1}\in[x_{i-2},x_{i-1}]$ and $\xi_{i+1}\in[x_{i},x_{i+1}]$, such that $\tilde{g}_{m}(x_{i-1}^{*})=g(\xi_{i-1})$ and $\tilde{g}_{m}(x_{i+1}^{*})=g(\xi_{i+1})$. Using the above results, we deduce the existence of three points $\xi_{i-1}<\xi_{i}<\xi_{i+1}$ such that $g(\xi_{i-1})=0=g(\xi_{i+1})$ and $g(\xi_{i})<0$, which is a contradiction since $g$ is concave by the concavity of $w$ and $l$. ∎

Lemma 3**.**

Let $w:[a,b]\rightarrow\mathbb{R}$ be a continuous concave function with $w_{+}^{\prime}(a)-w_{-}^{\prime}(b)\leq M$ and where $-\infty<a<b<\infty$. Then there exists a partition $a=x_{0}<x_{1}<\dots<x_{m}=b$ of $[a,b]$ with $\min_{i=1,\dots,m}(x_{i}-x_{i-1})\geq(b-a)(2m)^{-2}$ and such that

[TABLE]

where $\widetilde{w}_{m}$ is the piecewise linear approximation of $w$ defined in (16) and (17) and $C>0$ is a universal constant (i.e. not depending on $a,b,m$).

Proof.

By translation we may without loss of generality take $a=0$. Recall that since $w$ is a continuous concave function, it has left and right derivatives at every point $x\in[0,b]$. Define $\Delta w^{\prime}(x)=w_{-}^{\prime}(x)-w_{+}^{\prime}(x).$ For every $r\geq 1$, let $\mathcal{P}_{1}:=\{x_{i,1}:=\frac{ib}{r},i=0,\dots,r\}$ be the uniform partition of $[0,b]$ and let $\tilde{x}_{i,2}$, $i=1,\dots,r_{2}$, be the points such that $\Delta w^{\prime}(\tilde{x}_{i,2})\geq M/r$, setting $r_{2}=0$ if no such point exists. By concavity of $w$,

[TABLE]

so that $r_{2}\leq r$.

Consider a new partition $\mathcal{P}_{2}:=\{x_{0,2}<\dots<x_{r_{2}^{\prime},2}\}$ of $[0,b]$, consisting of the points $\{x_{i,1}\}\cup\{\tilde{x}_{i,2}\}\cup\{\tilde{x}_{i,2}-br^{-2}\}\cup\{\tilde{x}_{i,2}+br^{-2}\}$ written in increasing order. Note that $r_{2}^{\prime}\leq r+3r_{2}$. Colour in red all the points of the form $\tilde{x}_{i,2}$ and $\tilde{x}_{i,2}-br^{-2}$, so that each red point is the left endpoint of an interval of length at most $br^{-2}$. This colouring will be used to keep track of points that have a close neighbour on the right.

We next refine the partition $\mathcal{P}_{2}$ by adding the point $y$ between $x_{i,2}$ and $x_{i+1,2}$

[TABLE]

if

[TABLE]

Denote by $r_{3}$ the total number of points $y$ added in this manner to the sequence. We further add the points $y-br^{-2}$, $y+br^{-2}$ and colour in red all points of the form $y$ and $y-br^{-2}$, similarly to the previous case. Repeating this procedure results in a new partition that separates intervals where the derivative decreases by at most $2M/r$. Denote by $\mathcal{P}_{3}:=\{0=x_{0,3}<x_{1,3}<\dots<x_{r_{3}^{\prime},3}=b\}$ this new partition. We now show that $r_{3}^{\prime}\leq 7r$.

Let $y$ by any point added in the way just described. Suppose by contradiction that $w_{+}^{\prime}(x_{i,2})-w_{-}^{\prime}(y)<M/r$. By definition, we know that for all $x>y$, $w_{+}^{\prime}(x_{i,2})-w_{-}^{\prime}(x)>2M/r$. Subtracting the two inequalities gives $w_{-}^{\prime}(y)-w_{-}^{\prime}(x)>M/r$. However, since the right derivative of a concave function is right continuous, taking the limit $x\to y^{+}$ (and restricting to the points $x$ where $w$ is differentiable) yields $\Delta w^{\prime}(y)\geq M/r$. This is a contradiction however, because if this were the case, $y$ would already belong to $\mathcal{P}_{2}$. Since $w_{+}^{\prime}(x_{i,2})-w_{-}^{\prime}(y)\geq M/r$, using a similar argument to (19) gives $r_{3}\leq r$ so that $r_{3}^{\prime}\leq 7r$.

Finally, if the function $w$ is not differentiable at the point $x_{i,3}^{*}=\frac{x_{i,3}+x_{i-1,3}}{2}$, we split $[x_{i-1,3},x_{i,3}]$ into two parts in such a way that $w$ is differentiable at the midpoints of both new intervals and each interval has size at least $(x_{i,3}-x_{i-1,3})/3$. We add the points separating the new intervals to the previous partition, thereby obtaining $\mathcal{P}_{4}:=\{0=x_{0,4}<x_{1,4}<\dots<x_{\nu,4}=b\}$. The cardinality of $\mathcal{P}_{4}$, satisfies $\nu+1\leq 14r+1$. We now create a new partition $\mathcal{P}$ with polynomially separated points using the following algorithm.

Set $\mathcal{P}=\mathcal{P}_{4}$, keeping track of all the points coloured red. Set $\tilde{x}=x_{0,4}$. 2. 2.

If $b-\tilde{x}\leq br^{-2}$, remove all points in $\mathcal{P}$ strictly between $\tilde{x}$ and $b$ skip to Step 4. 3. 3.

Set $y=\inf\{t\in\mathcal{P}:t>\tilde{x}+br^{-2}\}$. Remove all elements of $\mathcal{P}$ between $\tilde{x}$ and $y$. If at least one element was removed, add to $\mathcal{P}$ the point $s=\tilde{x}+br^{-2}+\varepsilon$ for some $0<\varepsilon<br^{-2}\wedge(y-\tilde{x}-br^{-2})$ such that $w$ is differentiable at $(s+y)/2$ and $(s+\tilde{x})/2$. Colour $\tilde{x}$ red to mark that $s-\tilde{x}<2br^{-2}$. Set $\tilde{x}:=s$. If no point was removed from $\mathcal{P}$, set $\tilde{x}:=y$. Go to Step 2. 4. 4.

If $\tilde{x}=b$ then stop. Otherwise set $y=\max\{t\in\mathcal{P}:t<\tilde{x}\}$ and remove $\tilde{x}$ from $\mathcal{P}$. If $b-y>2br^{-2}$, add the point $s:=b-br^{-2}-\varepsilon$ to $\mathcal{P}$ and colour it red, where $0<\varepsilon<(b-y-2br^{-2})\wedge br^{-2}$ is such that $w$ is differentiable at $(y+s)/2$. If $b-y\leq 2br^{-2}$, add the point $s:=(y+b)/2$ and colour both $y$ and $s$ red.

Relabel the final partition $\mathcal{P}:=\{0=x_{0}<x_{1}<\dots<x_{m}=b\}$ and note that

[TABLE]

By construction $\min_{i=0,\dots,m-1}(x_{i+1}-x_{i})\geq\tfrac{1}{2}br^{-2}\geq\tfrac{1}{2}bm^{-2}$ and

[TABLE]

since if $x_{i}$ is coloured red, $x_{i+1}-x_{i}\leq 2br^{-2}\leq 2C_{0}^{2}bm^{-2}$.

We now show that $\|w-\widetilde{w}_{m}\|_{\infty}=O(bm^{-2})$. If $x_{i-1}$ is red, then by the mean value theorem, there exists $\xi_{i}\in J_{i}:=[x_{i-1},x_{i}]$ such that $\widetilde{w}_{m}(x_{i}^{*})=w(\xi_{i})$. Using the Lipschitz continuity of $w$ and $\widetilde{w}_{m}$,

[TABLE]

If $x_{i-1}$ is not red, Taylor expanding $w$ at the points $x_{i}^{*}:=\frac{x_{i}+x_{i-1}}{2}$ (at which $w$ is differentiable by the construction of $\mathcal{P})$ gives

[TABLE]

Due to the construction of $\mathcal{P}$,

[TABLE]

where $w^{\prime}(\xi_{i})$ stands here for some value in the interval $[w^{\prime}_{+}(\xi_{i}),w^{\prime}_{-}(\xi_{i})]$ for some point $\xi_{i}\in J_{i}$. We then deduce that $|R_{i}(x)|\leq\frac{2M}{r}\frac{C_{0}b}{2m}\leq C_{0}^{2}bM/m^{2}.$

Since $\widetilde{w}_{m}$ is piecewise linear, we can write

[TABLE]

where $\widetilde{w}_{m}^{\prime}$ denotes the left or right derivative of $\widetilde{w}_{m}$ at $x_{i}^{*}$, depending on whether $x<x_{i}^{*}$ or $x>x_{i}^{*}$. We now show that $|\widetilde{w}_{m}^{\prime}(x_{i}^{*})-w^{\prime}(x_{i}^{*})|\leq\max\{w_{+}^{\prime}(x_{i-1})-w^{\prime}(x_{i}^{*}),w^{\prime}(x_{i}^{*})-w_{-}^{\prime}(x_{i+1})\}$ for $i=1,\dots,m-1$. Consider the case of right derivatives (the same argument also works for left derivatives). Using the definition of $\widetilde{w}_{m}$ and that $\theta_{i}=\int_{J_{i}}w(s)ds$,

[TABLE]

and

[TABLE]

By the continuity and concavity of $w$, $(v-u)w_{-}^{\prime}(x_{i+1})\leq w(v)-w(u)\leq(v-u)w_{+}^{\prime}(x_{i-1})$ for any $x_{i-1}\leq u\leq v\leq x_{i+1}$. Combining all of the above yields

[TABLE]

We remark that $\max\{w_{+}^{\prime}(x_{i-1})-w^{\prime}(x_{i}^{*}),w^{\prime}(x_{i}^{*})-w_{-}^{\prime}(x_{i+1})\}\leq\frac{5M}{r}$. Indeed, since $x_{i}-x_{i-1}>C_{0}^{2}b/m^{2}$, the point $x_{i}$ is not equal to $x_{j,4}^{\prime}$ for any $j$ and hence both $\Delta w^{\prime}(x_{i})<M/r$ and $w_{+}^{\prime}(x_{i})-w_{-}^{\prime}(x_{i+1})\leq 2M/r$ hold. Together with $w_{+}^{\prime}(x_{i-1})-w^{\prime}(x_{i}^{*})\leq 2M/r$ and $w^{\prime}(x_{i}^{*})-w^{\prime}_{-}(x_{i})\leq 2M/r$, this verifies the preceding statement. Then

[TABLE]

We hence deduce that $|w^{\prime}(x_{i}^{*})-\widetilde{w}_{m}^{\prime}(x_{i}^{*})|\leq 5MC_{0}m^{-1}.$ Finally, using (20) and the fact that $\int_{J_{i}}(x-x_{i}^{*})dx=0$,

[TABLE]

Collecting together all the pieces, we have that for any $x\in[x_{i-1},x_{i}]$,

[TABLE]

∎

Lemma 4**.**

Any piecewise linear concave function $w:[a,b]\rightarrow\mathbb{R}$ with $N$ knots $\{z_{1},...,z_{N}\}$ can be written in the form

[TABLE]

with parameters $0\leq\gamma_{1}\leq(w_{+}^{\prime}(a)-w_{-}^{\prime}(b))(b-a)$, $|\gamma_{2}|\leq|w_{-}^{\prime}(b)|$, $\gamma_{3}\in\mathbb{R}$, $\sum_{i=1}^{N}p_{i}=1$ and $p_{i}\geq 0$ for $i=1,\dots,N$.

Proof.

The left derivative of $w$ is a step function $g:\,(a,b]\mapsto\mathbb{R}$ with $g(a+\varepsilon)=w_{+}^{\prime}(a)$ for sufficiently small $\varepsilon>0$ and $g(b)=w_{-}^{\prime}(b)$. By shifting this function vertically by $-w_{-}^{\prime}(b)$, we arrive at a non-negative, bounded, monotone decreasing step function, which can therefore be written as a monotone decreasing probability density times a normalizing constant $\gamma_{1}$. It is easy to see that $\gamma_{1}\leq(b-a)(w_{+}^{\prime}(a)-w_{-}^{\prime}(b))$. The step function can therefore be represented as

[TABLE]

with $P_{N}$ an atomic probability measure with $N$ atoms on $[0,b-a]$ and $\gamma_{1}$ the normalizing constant. Integrating the step function $g$ yields

[TABLE]

which is equal to $w$ for an appropriately chosen constant $C>0$. ∎

Proof of Proposition 1.

Let $\psi_{n}(x)=\mathds{1}_{[-s_{n},s_{n}]}(x)$ for $s_{n}=\tfrac{4}{5\alpha}\log n$. The log-concave density function $f_{1}=f_{1,n}=f_{0}\psi_{n}/\int f_{0}\psi_{n}$ supported on $[-s_{n},s_{n}]$ satisfies $|1-\int f_{0}\psi_{n}|\leq P_{f_{0}}[-s_{n},s_{n}]^{c}\leq 2e^{\beta}\alpha^{-1}n^{-4/5}$. Arguing as in (8), one has $h^{2}(f_{0},f_{1,n})\leq 12e^{\beta}\alpha^{-1}n^{-4/5}$ for $n\geq(4e^{\beta}/\alpha)^{5/4}$.

We write $f_{1}=e^{w_{1}}$ and construct the approximating function $\bar{f}_{n}$ according to the value of $w_{1}$ and its left and right derivatives $w_{1,-}^{\prime}$ and $w_{1,+}^{\prime}$. Let

[TABLE]

where $D>0$ is some fixed constant, $|w_{1,\pm}^{\prime}(x)|=\max(|w_{1,+}^{\prime}(x)|,|w_{1,-}^{\prime}(x)|)$ and $j_{n}=\lceil\log_{2}(n^{4/5}/D)\rceil-1$. In fact the set where the left and right derivatives of the concave function $w_{1}$ do not agree has measure zero. Note that the above sets are all disjoint except $A_{2,j_{n}}^{n}$ and $A_{3}^{n}$: since $j_{n}$ is the smallest integer such that $2^{-j_{n}-1}n^{4/5}\leq D$, these last two sets may overlap. In particular, we can express $[a_{n},b_{n}]$ as the almost disjoint union of the above sets. Write $B_{n}=(\cup_{j=0}^{j_{n}}A_{2,j}^{n})\cup A_{3}^{n}\subset[a_{n},b_{n}]$ and note that by the concavity of $w_{1}$, this is an interval. Since $\|f_{1}\|_{\infty}\leq 2e^{\beta}$ for $n\geq(4e^{\beta}/\alpha)^{5/4}$, the set $A_{1}^{n}$ consists of at most two intervals, each of width $O(n^{-4/5}\log n)$. Using again the boundedness of $f_{1}$, the definition of $A_{0}^{n}$ and that $|\text{supp}(f_{1})|\lesssim\log n$,

[TABLE]

so that in particular, $B_{n}\neq\emptyset$ for $n\geq n_{0}(\alpha,\beta)$ large enough.

We now construct a partition $\mathcal{P}_{n}$ of $B_{n}$ based on which we take the piecewise linear approximation (16)-(17) of the function $w_{1}$. Note that $A_{2,j}^{n}$ consists of at most two disjoint intervals, $A_{2,j,+}^{n}$ and $A_{2,j,-}^{n}$, which by the boundedness of $f_{1}$ are each of length $O(2^{j+1}n^{-4/5}\log n)$. Let $\mathcal{P}_{n,A_{2,j,+}^{n}}$ denote the partition of the interval $A_{2,j,+}^{n}$ given by Lemma 3 with partition size $m=\lceil 2^{-j/2}n^{3/5}|A_{2,j,+}^{n}|^{1/2}/\sqrt{\log n}\rceil=O(n^{1/5})$ and let $\mathcal{P}_{n,A_{2,j,-}^{n}}$ be the analogous partition constructed on $A_{2,j,-}^{n}$. Similarly, $A_{3}^{n}\cap(A_{2,j_{n}}^{n})^{c}$ consists of a single interval of length $O(\log n)$. Let $\mathcal{P}_{n,A_{3}^{n}\cap(A_{2,j_{n}}^{n})^{c}}$ denote the corresponding partition of $A_{3}^{n}\cap(A_{2,j_{n}}^{n})^{n}$ given by Lemma 3 with partition size $m=\lceil n^{1/5}(D|A_{3}^{n}\cap(A_{2,j_{n}}^{n})^{c}|/\log n)^{1/2}\rceil=O(n^{1/5}).$ Define the overall partition

[TABLE]

of $B_{n}$, which has $O(j_{n}n^{1/5})=O(n^{1/5}\log n)$ points. The associated piecewise linear function $\widetilde{w}_{n}$ defined in (16)-(17) based on $\mathcal{P}_{n}$ is concave by Lemma 2 and by construction corresponds to the partition given in Lemma 3 for each of the sets comprising $B_{n}$. It therefore satisfies the conclusions of Lemma 3 on each such set (with the appropriate $m$), so that in particular,

• •

$\|w_{1}-\widetilde{w}_{n}\|_{L^{\infty}(A_{2,j}^{n})}\leq Cn^{-2/5}\log n$ for some universal constant $C>0$ independent of $j$,

• •

the partition points in $A_{2,j}^{n}$ are distance at least $c2^{j}n^{-6/5}\log n\geq cn^{-6/5}\log n$ apart for some universal constant $c>0$ independent of $j$,

• •

on $A_{3}^{n}\cap(A_{2,j_{n}}^{n})^{c}$, we have the same $L^{\infty}$-bound with the partition points being $cn^{-2/5}\log n$-separated.

Moreover, since these intervals meet only at their boundaries, and the boundary points of the intervals are contained in the partition presented in Lemma 3, the interval boundaries will be contained in $\mathcal{P}_{n}$. Consequently, the separation property continues to hold even across the different subpartitions. In conclusion, we have shown that $\widetilde{w}_{n}$ is concave and piecewise linear with $O(n^{1/5}\log n)$ knots, which are $cn^{-6/5}\log n$-separated, and satisfies

[TABLE]

We now extend the approximating function to $[a_{n},b_{n}]\supset B_{n}$. Write $\mathcal{P}_{n}=(x_{i})_{i=0}^{M}$, where $\min(B_{n})=x_{0}<x_{1}<...<x_{M}=\max(B_{n})$ and $M=O(n^{1/5}\log n)$. Define $\bar{w}_{n}:[a_{n},b_{n}]\rightarrow\mathbb{R}$ as

[TABLE]

This is simply the function $\widetilde{w}_{n}$ extended linearly from the boundary points of $B_{n}$ with slope $w_{1,-}^{\prime}(x_{0})\wedge n^{4/5}\vee(-n^{4/5})$ and $w_{1,+}^{\prime}(x_{M})\wedge n^{4/5}\vee(-n^{4/5})$ on $[a_{n},x_{0}]$ and $[x_{M},b_{n}]$ respectively. We now verify that $\bar{w}_{n}$ is concave, for which it is enough to show that $\bar{w}_{n,+}^{\prime}(x_{0})\leq\bar{w}_{n,-}^{\prime}(x_{0})$ and $\bar{w}_{n,+}^{\prime}(x_{M})\leq\bar{w}_{n,-}^{\prime}(x_{M})$. For the first inequality, using (21), the concavity of $w_{1}$ and the boundary construction of $\widetilde{w}_{n}$ given by (16), $\bar{w}_{n,+}^{\prime}(x_{0})=\widetilde{w}_{n,+}^{\prime}(x_{0})=\widetilde{w}_{n,+}^{\prime}(x_{1}^{*})\leq w_{1,+}^{\prime}(x_{0})\leq w_{1,-}^{\prime}(x_{0})$. Since $x_{0}\in B_{n}$, it also holds that $|\bar{w}_{n,+}^{\prime}(x_{0})|\leq n^{4/5}$. The second inequality can be proved analogously.

Since $\log(1+z)=O(z)$ as $|z|\rightarrow 0$, it follows that $\log\int f_{0}\psi_{n}=O(n^{-4/5})$. Using this and (23),

[TABLE]

By concavity, the slope of the linear extension on $[a_{n},x_{0}]$ satisfies $w_{1,-}^{\prime}(x_{0})=(\log f_{0})_{-}^{\prime}(x_{0})\leq(\log f_{0})_{+}^{\prime}(x)$ for all $x<x_{0}$ such that $f_{0}(x)>0$. Combining the above yields $\bar{w}_{n}(x)\geq\log f_{0}(x)-O(n^{-2/5}\log n)$ for all $x\in[a_{n},x_{0}]$. The same computation also gives the result for $x\in[x_{M},b_{n}]$, so that for some $C>0$,

[TABLE]

Define the log-concave density

[TABLE]

This function is piecewise log-linear, has $O(n^{1/5}\log n)$ knots and satisfies $(ii)$ and $(iii)$ by construction. We have

[TABLE]

The first integral is $O(n^{-4/5}\log n)$ by (22). Using (23), $\int_{B_{n}}e^{\bar{w}_{n}}=e^{o(1)}\int_{B_{n}}f_{1}$. Write the second integral as $\int_{B_{n}^{c}}\bar{f}_{n}=\int_{A_{0}^{n}}\bar{f}_{n}+\int_{A_{1}^{n}}\bar{f}_{n}$. By the definition of $A_{0}^{n}$ and (25), $\int_{A_{0}^{n}}e^{\bar{w}_{n}}\leq\int_{A_{0}^{n}}e^{-(4/5)\log n+n^{-2/5}\log n}\leq(x_{0}-a_{n}+b_{n}-x_{M})e^{o(1)}n^{-4/5}=O((b_{n}-a_{n})n^{-4/5})$. For the integral over $A_{1}^{n}$ we simply observe that by (25), $\bar{w}_{n}\leq\beta+n^{-2/5}\log n$, and recall that the measure of $A_{1}^{n}$ is at most $2n^{-4/5}(\beta+\frac{4}{5}\log n)$. Since $b_{n}-a_{n}\geq\frac{16}{5\alpha}\log n$, then $\int_{A_{1}^{n}}\bar{f}_{n}$ is also $O((b_{n}-a_{n})n^{-4/5})$. This implies that the second integral in (26) is $O((b_{n}-a_{n})n^{-4/5})$.

Using (22), (23), Lemma 3.1 of [36] and the above, the third term of (26) is bounded by a multiple of

[TABLE]

which establishes $(i)$.

Consider $(iv)$. Note that this is trivial if $f_{0}(x)=0$, so assume $f_{0}(x)\neq 0$. If $x\in B_{n}$, then by (23),

[TABLE]

If $x\in[a_{n},b_{n}]\backslash B_{n}$, then the result follows from (25).

Consider lastly $(v)$. Since $\bar{w}_{n}$ defined in (24) is piecewise linear with $|w_{+}^{\prime}(a_{n})|\vee|w_{-}^{\prime}(b_{n})|\leq n^{4/5}$, in view of Lemma 4 it takes the form

[TABLE]

with $M=O(n^{1/5}\log n)$, $\gamma_{1}\leq|w_{+}^{\prime}(a_{n})-w_{-}^{\prime}(b_{n})|(b_{n}-a_{n})\leq 2n^{4/5}(b_{n}-a_{n})$ and $|\gamma_{2}|\leq|w_{-}^{\prime}(b_{n})|\leq n^{4/5}$. This completes the proof. ∎

6.2 Proof of Lemma 1

Proof of Lemma 1.

It is shown in the proof of Theorem 5 of Kim and Samworth [23], p. 2772, that for $\eta=\eta(\varepsilon)\in(0,1)$ to be defined below and $L\geq L_{0}(\eta)=L_{0}(\varepsilon)$,

[TABLE]

so that it remains only to control $\sup_{g_{0}\in\overline{\mathcal{F}}^{0,1}}P_{g_{0}}^{n}(\hat{g}_{n}\not\in\widetilde{\mathcal{F}}^{1,\eta})$. Lemma 6 of Kim and Samworth [23] shows that this quantity is $O(n^{-1})$ as $n\to\infty$; we essentially follow their proof, suitably sharpening the probability bounds in the case $d=1$. Then

[TABLE]

By Lemma 13 of [11], there exist universal constants $\alpha_{0},\beta_{0}>0$ such that for all $x\in\mathbb{R}$,

[TABLE]

It is shown in [23], p. 2774, that the first term in (28) is bounded by $2\alpha_{0}^{-1}e^{\beta_{0}-\alpha_{0}\sqrt{n}}$. For the second term, by Remark 2.3 of Dümbgen et al. [10], one has $\sigma_{\hat{g}_{n}}^{2}\leq\widetilde{\sigma}_{n}^{2}:=n^{-1}\sum_{i=1}^{n}(Z_{i}-\bar{Z}_{n})^{2}\leq n^{-1}\sum_{i=1}^{n}Z_{i}^{2}$, where $\widetilde{\sigma}_{n}^{2}$ denotes the sample variance and $\bar{Z}_{n}=n^{-1}\sum_{i=1}^{n}Z_{i}$ the sample mean. Letting $C_{0}=C_{0}(\alpha_{0},\beta_{0},2)$ denote the constant in Lemma 5 below, we can apply that lemma to bound the second term in (28) by

[TABLE]

for $n\geq\max(16C_{0}^{2}/\eta,e^{2})$.

Consider now the third term in (28). For $\varepsilon>0$, let $\overline{\cal{P}}_{\varepsilon}^{1/10,1/2}$ denote the class of probability distributions on $\mathbb{R}$ such that $\mu_{P}=\int xdP(x)$ and $\sigma_{P}^{2}=\int(x-\mu_{P})^{2}dP(x)$ satisfy $|\mu_{P}|\leq 1/10$ and $1/2\leq\sigma_{P}^{2}\leq 3/2$ and

[TABLE]

where $\Gamma(t)=\int_{0}^{\infty}x^{t-1}e^{-x}dx$. This is exactly the same as the class $\mathcal{P}^{1/10,1/2}$ considered in Lemma 6 of [23], except that we have replaced the $4^{\text{th}}$-moment condition with a $(2+\varepsilon)$-moment condition. Following the rest of the proof of Lemma 6 of [23] (noting that $\sup_{g_{0}\in\overline{\mathcal{F}}^{0,1}}\int|x|^{2+\varepsilon}g_{0}(x)dx\leq\tau_{\varepsilon}/2$ and that uniform integrability of $\{Y_{n_{k}}^{2}:k\in\mathbb{N}\}$ in that proof follows from the $(2+\varepsilon)$-moment condition), one can similarly conclude that for some \eta=\eta\big{(}\alpha_{0},\beta_{0},\varepsilon\big{)}\in(0,1),

[TABLE]

where $\mathbb{P}_{n}=n^{-1}\sum_{i=1}^{n}\delta_{Z_{i}}$ is the empirical measure and $\delta_{x}$ is the Dirac measure at $x$. This last probability can be bounded by

[TABLE]

Using similar arguments to those used previously, the first two terms can be bounded by $Ce^{-c\sqrt{n}}$ and $Ce^{-cn^{1/4}}$, respectively, where the constants depend only on $\alpha_{0}$ and $\beta_{0}$. For the last term, using Lemma 5 with $C_{0}=C_{0}(\alpha_{0},\beta_{0},2+\varepsilon)$ the constant in that lemma,

[TABLE]

for $n\geq\max(4(2+\varepsilon)^{4+2\varepsilon}C_{0}^{2}/\tau_{\varepsilon}^{2},e^{2+\varepsilon})$. In conclusion, we have shown that for any $\varepsilon>0$ there exists $\eta=\eta(\varepsilon)\in(0,1)$ such that

[TABLE]

for all $n\geq n_{0}(\varepsilon)$, where $C(\varepsilon)$, $c(\varepsilon)$ and $n_{0}(\varepsilon)$ are positive constants. Together with (27), this establishes the result. ∎

Lemma 5**.**

Let $Z_{1},...,Z_{n}$ be i.i.d. random variables from a density $f\in\mathcal{F}_{\alpha,\beta}$, where $\alpha>0$ and $\beta\in\mathbb{R}$. Then for any $r\geq 1$, $t\geq r^{r}$ and $n\geq e^{r}$,

[TABLE]

Proof.

For notational convenience, write $\lambda=1/r\in(0,1]$. Let $x_{\lambda}=(1/\lambda)^{1/\lambda}$ and define

[TABLE]

where $\tau_{\lambda}=\Psi_{\lambda}(x_{\lambda})/x_{\lambda}=(\lambda e)^{1/\lambda}$. This defines a Young function, that is a convex, increasing function $\Psi_{\lambda}:\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$ with $\Psi_{\lambda}(0)=0$. Denote the corresponding Orlicz norm $\|X\|_{\Psi_{\lambda}}:=\inf\{a>0:\mathbb{E}[\Psi_{\lambda}(|X|/a)]\leq 1\}$. Note that the density function $g_{\lambda}$ of $|Z_{1}|^{1/\lambda}$ satisfies g_{\lambda}(x)=\lambda x^{-(1-\lambda)}\big{(}f(x^{\lambda})+f(-x^{\lambda})\big{)}\mathds{1}_{\{x\geq 0\}}\leq\lambda x^{-(1-\lambda)}e^{\beta-\alpha x^{\lambda}}\mathds{1}_{\{x\geq 0\}}. Then for fixed $a>\alpha^{-1/\lambda}$,

[TABLE]

If $K_{0}=K_{0}(a,\alpha,\beta,\lambda)>1$, it follows by convexity that $\||Z_{1}|^{1/\lambda}\|_{\Psi_{\lambda}}\leq aK_{0}$.

By Theorem 6.21 of Ledoux and Talagrand [25], there exists a constant $K_{\lambda}$ such that

[TABLE]

The first-term on the right-hand side can be bounded by the $\|\cdot\|_{2}$-norm of the same quantity, which equals the square root of $n\mathbb{E}(|Z_{i}|^{\frac{1}{\lambda}}-\mathbb{E}|Z_{i}|^{\frac{1}{\lambda}})^{2}\leq n\mathbb{E}|Z_{i}|^{\frac{2}{\lambda}}$. For any $\delta\in(0,1]$,

[TABLE]

Since this is finite for any $\delta\in(0,1]$, we can bound the $\|\cdot\|_{1}$-norm above by $C(\alpha,\beta,\lambda)\sqrt{n}$.

Note that for any random variable $X$, $\|X-\mathbb{E}X\|_{\Psi_{\lambda}}\leq 2\|X\|_{\Psi_{\lambda}}$. Indeed, setting $a=\|X\|_{\Psi_{\lambda}}$, since $\Psi_{\lambda}$ is convex and increasing,

[TABLE]

Using this and Lemma 2.2.2 of van der Vaart and Wellner [37],

[TABLE]

so that we have shown

[TABLE]

By Markov’s inequality, for any random variable $X$, $\mathbb{P}(|X|\geq x\|X\|_{\Psi_{\lambda}})=\mathbb{P}(\Psi_{\lambda}(|X|/\|X\|_{\Psi_{\lambda}})\geq\Psi_{\lambda}(x))\leq 1/\Psi_{\lambda}(x)$. Applying this to the above sum completes the proof. ∎

7 Additional simulations

We firstly provide some additional figures for the empirical Bayes posterior of the mode. We consider the $Beta(2,5)$ and $Gamma(2,1)$ distributions and take i.i.d. samples of size ranging from $n=50$ to $n=20000$ in both cases. As in the Gaussian example, we run the Gibbs sampler for $20000$ iterations of which half were considered as burn in and discarded. The resulting marginal posteriors for the mode based on the empirical Bayes posterior are displayed in Figures 8 and 9.

Next we provide a brief simulation study to demonstrate the (possible) applicability of our procedure for cluster analysis. Here we consider only the simple case where the underlying density is a mixture of two log-concave densities. We have modified our Bayesian procedure to accommodate mixtures of log-concave densities in a straightforward way, again using a Gibbs sampler. We consider various combinations of log-concave densities including two Gaussians, two Laplace distributions, a Laplace and a Gamma distribution and two Beta distributions. We consider two sample sizes, $n=100$ and $n=500$, see Figures 10 and 11, respectively. In all pictures we have plotted the posterior mean (solid blue), $95\%$ pointwise credible sets (dashed blue) and the underlying density (solid red). One can see that in all cases our procedure provides reasonable estimators and seemingly reliable uncertainty quantification.

Acknowledgements: The authors would like to thank Richard Samworth and Arlene Kim for helpful discussions. The authors would further like to thank the AE and two referees for their helpful suggestions which lead to an improved version of the manuscript.

All three authors received funding from the European Research Council under ERC Grant Agreement 320637 for this research. Ester Mariucci was further supported by the Federal Ministry for Education and Research through the Sponsorship provided by the Alexander von Humboldt Foundation, by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 314838170, GRK 2297 MathCoRe, and by Deutsche Forschungsgemeinschaft (DFG) through grant CRC 1294 ’Data Assimilation’. Botond Szabó also received funding from the Netherlands Organization for Scientific Research (NWO) under Project number: 639.031.654.