Maximum likelihood estimation of a multidimensional log-concave density

Introduction

Modern nonparametric density estimation began with the introduction of a kernel density estimator in the pioneering work of Fix and Hodges (1951), later republished as Fix and Hodges (1989). For independent and identically distributed real-valued observations, the appealing asymptotic theory of the mean integrated squared error was provided by Rosenblatt (1956) and Parzen (1962). This theory leads to an asymptotically optimal choice of the smoothing parameter, or bandwidth. Unfortunately, however, it depends on the unknown density $f$ through the integral of the square of the second derivative of $f$. Considerable effort has therefore been focused on finding methods of automatic bandwidth selection (cf. Wand and Jones, 1995, Chapter 3, and the references therein). Although this has resulted in algorithms, e.g. Chiu (1992), that achieve the optimal rate of convergence of the relative error, namely $O_{p}(n^{-1/2})$, where $n$ is the sample size, good finite sample performance is by no means guaranteed.

In this paper, we propose a fully automatic nonparametric estimator of $f$, with no tuning parameters to be chosen, under the condition that $f$ is log-concave – that is, $\log f$ is a concave function. The class of log-concave densities has many attractive properties and has been well-studied, particularly in the economics, sampling and reliability theory literature. See Section 2 for further discussion of examples, applications and properties of log-concave densities.

In Section 3, we show that if $X_{1},\ldots,X_{n}$ are independent and identically distributed random vectors with a log-concave density, then with probability one there exists a unique log-concave density $\hat{f}_{n}$ that maximises the likelihood function,

Before continuing, it is worth noting that without any shape constraints on the densities under consideration, the likelihood function is unbounded. To see this, we could define a sequence $(f_{n})$ of densities that represent successively close approximations to a mixture of $n$ ‘spikes’ (one on each $X_{i}$), such as $f_{n}(x)=n^{-1}\sum_{i=1}^{n}\phi_{d,n^{-1}I}(x-X_{i})$, where $\phi_{d,\Sigma}$ denotes the $N_{d}(0,\Sigma)$ density. This sequence satisfies $L(f_{n})\rightarrow\infty$ as $n\rightarrow\infty$ (cf. Figure 2). In fact, a modification of this argument may be used to show that the likelihood function remains unbounded even if we restrict attention to unimodal densities.

Figure 2 gives a diagram illustrating the structure of the maximum likelihood estimator on the logarithmic scale. This structure is most easily visualised for two-dimensional data, where one can imagine associating a ‘tent pole’ with each observation, extending vertically out of the plane. For certain tent pole heights, the graph of the logarithm of the maximum likelihood estimator can be thought of as the roof of a taut tent stretched over the tent poles. The fact that the logarithm of the maximum likelihood estimator is of this ‘tent function’ form constitutes part of the proof of its existence and uniqueness.

In Section 4, we discuss the computational problem of how to adjust the $n$ tent pole heights so that the corresponding tent functions converge to the logarithm of the maximum likelihood estimator. One reason that this computational problem is so challenging in more than one dimension is the fact that it is difficult to describe the set of tent pole heights that correspond to concave functions. The key observation, discussed in Section 4, is that it is possible to minimise a modified objective function that it is convex (though non-differentiable). This allows us to apply the powerful non-differentiable convex optimisation methodology of the subgradient method (Shor, 1985) and a variant called Shor’s $r$-algorithm, which has been implemented by Kappel and Kuntsevich (2000).

As an illustration of the estimates obtained, Figure 3 presents plots of the maximum likelihood estimator, and its logarithm, for 1000 observations from a standard bivariate normal distribution. These plots were created using the LogConcDEAD package (Cule et al., 2008a) in R (R Development Core Team, 2008), which exploits the interactive surface-plotting software available in the rgl package (Adler and Murdoch, 2007).

In Section 5 we present simulations to compare the finite-sample performance of the maximum likelihood estimator with kernel-based methods. The results are striking: even when we use the theoretical, optimal bandwidth for the kernel estimator (or an asymptotic approximation to this when it is not available), we find that the maximum likelihood estimator has a rather smaller mean integrated squared error for moderate or large sample sizes, despite the fact that this optimal bandwidth depends on properties of the density that would be unknown in practice. This suggests that the maximum likelihood estimator is able to adapt to the local smoothness of the underlying density automatically.

Nonparametric density estimation is a fundamental tool for the visualisation of structure in exploratory data analysis, and has an enormous literature that includes the monographs of Devroye and Györfi (1985), Silverman (1986), Scott (1992) and Wand and Jones (1995). Our proposed method may certainly be used for this purpose; however, it may also be used as an intermediary stage in more involved statistical procedures. For instance:

Clustering problems are closely related to the classification problems described above. The difference is that, in the above notation, we do not observe $Y_{1},\ldots,Y_{n}$, and have to assign each of $X_{1},\ldots,X_{n}$ to one of the $p$ populations. A common technique is based on fitting a mixture density of the form $f(x)=\sum_{j=1}^{p}\pi_{j}f_{j}(x)$, where the mixture proportions $\pi_{1},\ldots,\pi_{p}$ are positive and sum to one. Under the assumption that each of the component densities $f_{1},\ldots,f_{p}$ is log-concave, we show in Section 6 that our methodology can be extended to fit such a finite mixture density, which need not itself be log-concave – cf. Section 2. We also illustrate this clustering algorithm on a Wisconsin breast cancer data set in Section 6, where the aim is to separate observations into benign and malignant component populations.

A functional of the true underlying density may be estimated by the corresponding functional of a density estimator, such as the log-concave maximum likelihood estimator. Examples of functionals of interest include probabilities, such as $\int_{\|x\|\geq 1}f(x)\,dx$, moments, e.g. $\int\|x\|^{2}f(x)\,dx$, and the differential entropy, $-\int f(x)\log f(x)\,dx$. It may be possible to compute the plug-in estimator based on the log-concave maximum likelihood estimator analytically, but in Section 7, we show that even if this is not possible, in many cases of interest we can sample from the log-concave maximum likelihood estimator $\hat{f}_{n}$, and hence obtain a Monte Carlo estimate of the functional. This nice feature also means that the log-concave maximum likelihood estimator can be used in a Monte Carlo bootstrap procedure for assessing uncertainty in functional estimates – see Section 7 for further details.

The fitting of a nonparametric density estimate may give an indication of the validity of a particular smaller model (often parametric). Thus, a contour plot of the log-concave maximum likelihood estimator may provide evidence that the underlying density has elliptical contours, and thus suggest that a model that exploits this elliptical symmetry.

In the univariate case, Walther (2002) describes methodology based on log-concave density estimation for addressing the problem of detecting the presence of mixing in a distribution. As an application, he cites the Pickering/Platt debate (Swales, 1985) on the issue of whether high blood pressure is a disease (in which case observed blood pressure measurements should follow a mixture distribution), or simply a label attached to people in the right tail of the blood pressure distribution. As a result of our algorithm for computing the multidimensional log-concave maximum likelihood estimator, this methodology extends immediately to more than one dimension.

There has been considerable recent interest in shape-restricted nonparametric density estimation, but most of it has been confined to the case of univariate densities, where the computational algorithms are more straightforward. Nevertheless, as was discussed above, it is in multivariate situations that the automatic nature of the maximum likelihood estimator is particularly valuable. Walther (2002), Dümbgen and Rufibach (2007) and Pal et al. (2007) have proved the existence and uniqueness of the log-concave maximum likelihood estimator in one dimension and Dümbgen and Rufibach (2007), Pal et al. (2007) and Balabdaoui et al. (2008) have studied its theoretical properties. Rufibach (2007) has compared different algorithms for computing the univariate estimator, including the iterative convex minorant algorithm (Groeneboom and Wellner, 1992; Jongbloed, 1998), and three others. Dümbgen et al. (2007) also present an Active Set algorithm, which has similarities with the vertex direction and vertex reduction algorithms described in Groeneboom et al. (2008). For univariate data, it is also well-known that there exist maximum likelihood estimators of a non-increasing density supported on $[0,\infty)$ (Grenander, 1956) and of a convex, decreasing density (Groeneboom et al., 2001).

In Section 8, we give a brief concluding discussion, and suggest some directions for future research. Finally, we present in Appendix A a glossary of terms and results from convex analysis and computational geometry that appear in italics at their first occurrence in the main body of the paper; the references are Rockafellar (1997) and Lee (1997). Proofs are deferred to Appendix B, except that the beginning of the proof of Theorem 2 is given in the main text, as the ideas and notation introduced are needed in the remainder of the paper.

Log-concave densities: examples, applications and properties

The assumption of log-concavity is a popular one in economics; Caplin and Naelbuff (1991b) show that in the theory of elections and under a log-concavity assumption, the proposal most preferred by the mean voter is unbeatable under a 64% majority rule. As another example, in the theory of imperfect competition, Caplin and Naelbuff (1991a) use log-concavity of the density of consumers’ utility parameters as a sufficient condition in their proof of the existence of a pure-strategy price equilibrium for any number of firms producing any set of products. See Bagnoli and Bergstrom (1989) for many other applications of log-concavity to economics. Brooks (1998) and Mengersen and Tweedie (1996) have exploited the properties of log-concave densities in studying the convergence of Markov chain Monte Carlo sampling procedures.

An (1998) lists many useful properties of log-concave densities. For instance, if $f$ and $g$ are (possibly multidimensional) log-concave densities, then their convolution $f\ast g$ is log-concave. In other words, if $X$ and $Y$ are independent and have log-concave densities, then their sum $X+Y$ has a log-concave density. The class of log-concave densities is also closed under the taking of pointwise limits. One-dimensional log-concave densities have increasing hazard functions, which is why they are of interest in reliability theory. Moreover, Ibragimov (1956) proved the following characterisation: a univariate density $f$ is log-concave if and only if the convolution $f\ast g$ is unimodal for every unimodal density $g$. There is no natural generalisation of this result to higher dimensions.

As was mentioned in Section 1, this paper concerns multidimensional log-concave densities, for which fewer properties are known. It is therefore of interest to understand how the property of log-concavity in more than one dimension relates to the univariate notion. Our first proposition below is intended to give some insight into this issue. It is not formally required for the subsequent development of our methodology in Sections 3 and 4, although we did apply the result when designing our simulation study in Section 5. We assume throughout that log-concave densities are with respect to Lebesgue measure on the affine hull of their support, and ‘$X$ has a log-concave density’ means ‘there exists a version of the density of $X$ that is log-concave’.

necessary that for any subspace $V$, the marginal density of $P_{V}(X)$ is log-concave and the conditional density $f_{X|P_{V}(X)}(\cdot|t)$ of $X$ given $P_{V}(X)=t$ is log-concave for each $t$

sufficient that for every $(d-1)$-dimensional subspace $V$, the conditional density $f_{X|P_{V}(X)}(\cdot|t)$ of $X$ given $P_{V}(X)=t$ is log-concave for each $t$.

The part of Proposition 1(a) concerning marginal densities is an immediate consequence of Theorem 6 of Prékopa (1973). One can regard Proposition 1(b) as saying that a multidimensional density is log-concave if the restriction of the density to any line is a (univariate) log-concave function.

It is interesting to compare the properties of log-concave densities presented in Proposition 1 with the corresponding properties of Gaussian densities. In fact, Proposition 1 remains true if we replace ‘log-concave’ with ‘Gaussian’ throughout (at least, provided that in part (b) we also assume there is a point at which $f$ is twice differentiable). These shared properties suggest that the class of log-concave densities is a natural, infinite-dimensional generalisation of the class of Gaussian densities.

Existence, uniqueness and structure of the maximum likelihood estimator

Suppose that $n\geq d+1$. Then, with probability one, a nonparametric maximum likelihood estimator $\hat{f}_{n}$ of $f_{0}$ exists and is unique.

Suppose that $f$ maximises $\psi_{n}(\cdot)$ over $\mathcal{F}$. The main part of the proof, which is completed in the Appendix, consists of showing that

there exists $M>0$ such that if $\max_{i}|\bar{h}_{y}(X_{i})|\geq M$, then \psi_{n}\bigl{(}\exp(\bar{h}_{y})\bigr{)}\leq\psi_{n}(f).

Although step (iii) above gives us a finite-dimensional class of functions to which $\log\hat{f}_{n}$ belongs, the proof of Theorem 2 gives no indication of how to find the member of this class that maximises the likelihood function. We therefore seek an iterative algorithm to compute the estimator, but first we describe the structure we see in Figure 2 in Section 1 more precisely. From now on, we assume:

$n\geq d+1$, and every subset of $\{X_{1},\ldots,X_{n}\}$ of size $d+1$ is affinely independent.

the relative interiors of the sets $\{C_{n,j}:j\in J\}$ are pairwise disjoint

In the iterative algorithm that we propose in Section 4 for computing the maximum likelihood estimator, we need to find convex hulls and triangulations at each iteration. Fortunately, these can be computed efficiently using the Quickhull algorithm of Barber et al. (1996).

Computation of the maximum likelihood estimator

As a first attempt to find an algorithm which produces a sequence that converges to the maximum likelihood estimator in Theorem 2, it is natural to try to minimise numerically the function

Although this approach might work in principle, one difficulty is that $\tau$ is not convex, so this approach is extremely computationally intensive, even with relatively few observations. Another reason for the numerical difficulties stems from the fact that the set of $y$-values on which $\tau$ attains its minimum is rather large: in general it may be possible to alter particular components $y_{i}$ without changing $\bar{h}_{y}$. Of course, we could have defined $\tau$ as a function of $\bar{h}_{y}$ rather than as a function of the vector of tent pole heights $y=(y_{1},\ldots,y_{n})$. Our choice, however, motivates the following definition of a modified objective function:

The great advantages of minimising $\sigma$ rather than $\tau$ are seen by the following theorem.

Thus Theorem 3 shows that the unique minimum $y^{*}=(y_{1}^{*},\ldots,y_{n}^{*})$ of $\sigma$ belongs to the minimum set of $\tau$. In fact, it corresponds to the element of the minimum set for which $\bar{h}_{y^{*}}(X_{i})=y_{i}^{*}$ for $i=1,\ldots,n$. Informally, then, $\bar{h}_{y^{*}}$ is ‘a tent function with all of the tent poles touching the tent’.

For each $j=(j_{1},\ldots,j_{d+1})\in J$, let $A_{j}$ be the $d\times d$ matrix whose $l$th column is $X_{j_{l+1}}-X_{j_{1}}$ for $l=1,\ldots,d$, and let $\alpha_{j}=X_{j_{1}}$. Then the affine transformation $w\mapsto A_{j}w+\alpha_{j}$ takes the unit simplex T_{d}=\bigl{\{}w=(w_{1},\ldots,w_{d}):w_{l}\geq 0,\sum_{l=1}^{d}w_{l}\leq 1\bigr{\}} to $C_{n,j}$. Letting $z_{j,l}=y_{j_{l+1}}-y_{j_{1}}$, we can then establish by a simple change of variables and induction on $d$ that if $z_{j,1},\ldots,z_{j,d}$ are non-zero and distinct, then

Further details of this calculation can be found in a longer version of this paper (Cule et al., 2008b). The singularities that occur when some of $z_{j,1},\ldots,z_{j,d}$ may be zero or equal are removable. Thus, although (4.2) is a little complicated, it allows the computation of our objective function.

2 Nonsmooth optimisation

There is a vast literature on techniques of convex optimisation (cf. Boyd and Vandenberghe (2004), for example), including the method of steepest descent and Newton’s method. Unfortunately, these methods rely on the differentiability of the objective function, and the function $\sigma$ is not differentiable. This can be seen informally by studying the schematic diagram in Figure 2 again. If the $i$th tent pole, say, is touching but not critically supporting the tent, then decreasing the height of this tent pole does not change the tent function, and thus does not alter the integral in (4.1); on the other hand, increasing the height of the tent pole does alter the tent function and therefore the integral in (4.1). This argument may be used to show that at such a point, the $i$th partial derivative of $\sigma$ does not exist.

Shor recognised, however, that the convergence of this algorithm could be slow in practice, and that although appropriate step size selection could improve matters somewhat, the convergence would never be better than linear (compared with quadratic convergence for Newton’s method near the optimum – see Boyd and Vandenberghe (2004, Section 9.5)). Slow convergence can be caused by taking at each stage a step in a direction nearly orthogonal to the direction towards the optimum, which means that simply adjusting the step size selection scheme will never produce the desired improvements in convergence rate.

One solution (Shor, 1985, Chapter 3) is to attempt to shrink the angle between the subgradient and the direction towards the minimum through a (necessarily nonorthogonal) linear transformation, and perform the subgradient step in the transformed space. By analogy with Newton’s method for smooth functions, an appropriate transformation would be an approximation to the inverse of the Hessian matrix at the optimum. This is not possible for nonsmooth problems, because the inverse might not even exist (and will not exist at points at which the function is not differentiable, which may include the optimum).

Instead, we perform a sequence of dilations in the direction of the difference between two successive subgradients, in the hope of improving convergence in the worst-case scenario of steps nearly perpendicular to the direction towards the minimiser. This variant, which has become known as Shor’s $r$-algorithm, has been implemented in Kappel and Kuntsevich (2000). Accompanying software SolvOpt is available from http://www.uni-graz.at/imawww/kuntsevich/solvopt/.

Although the formal convergence of the $r$-algorithm has not been proved, we agree with the authors’ claims that it is robust, efficient and accurate. Of course, it is clear that if we terminate the $r$-algorithm after any finite number of steps and apply the original Shor algorithm using our terminating value of $y$ as the new starting value, then formal convergence is guaranteed. We have not found it necessary to run the original Shor algorithm after termination of the $r$-algorithm in practice.

for some small $\delta,\epsilon\textrm{ and }\eta>0$. The first two termination criteria follow Kappel and Kuntsevich (2000), while the third is based on our knowledge that the true optimum corresponds to a density (Section 3). As default values, and throughout this paper, we took $\delta=10^{-8}$ and $\epsilon=\eta=10^{-4}$.

Table 1 gives approximate running times and number of iterations of Shor’s $r$-algorithm required for different sample sizes and dimensions on an ordinary desktop computer (1.8GHz, 2GB RAM). Unsurprisingly, the running time increases relatively quickly with the sample size, while the number of iterations increases approximately linearly with $n$. Each iteration takes longer as the dimension increases, though it is interesting to note that the number of iterations required for the algorithm to terminate decreases as the dimension increases. When $d=1$, we recommend the Active Set algorithm of Dümbgen et al. (2007), which is implemented in the R package logcondens (Rufibach and Dümbgen, 2006).

Finite sample performance

Our simulation study considered, for $d=2$ and $3$, the following densities:

standard normal, $\phi_{d}\equiv\phi_{d,I}$

the joint density of independent $\Gamma(2,1)$ components

the normal location mixture $0.6\phi_{d}(\cdot)+0.4\phi_{d}(\cdot-\mu)$ for (d) $\|\mu\|=1$, (e) $\|\mu\|=2$, (f) $\|\mu\|=3$. An application of Proposition 1 gives that such a normal location mixture is log-concave if and only if $\|\mu\|\leq 2$.

In Tables 2 and 3 we present, for each density and for four different sample sizes, an estimate of the mean integrated squared error (MISE) of the nonparametric maximum likelihood estimator based on 100 Monte Carlo iterations. We also show the MISE for the kernel density estimates with a Gaussian kernel and, for all of the normal and mixture of normal examples, the choice of bandwidth that minimises the MISE. In the gamma example, exact MISE calculations are not possible, so we took the bandwidth that minimises the asymptotic mean integrated squared error (AMISE). These optimal bandwidths can be computed using the formulae in Wand and Jones (1995, Sections 4.3 and 4.4). As minimisation of the expressions for both the MISE and the AMISE requires knowledge of certain functionals of the true density that would be unknown in practice, we also provide a comparison with an empirical bandwidth selector based on least squares cross validation (LSCV) (Wand and Jones, 1995, Section 4.7). The LSCV bandwidths were computed using the ks package (Duong, 2007) in R, and we used the option of constraining the bandwidth matrices to be diagonal in cases (a) and (c) where the components are independent.

We see that in cases (a)-(e) the log-concave maximum likelihood estimator has a smaller MISE than the kernel estimate with bandwidth chosen by LSCV, and at least for moderate and large sample sizes, the difference is quite dramatic. Even more remarkably, in these cases the log-concave estimator also outperforms the kernel estimate with optimally chosen bandwidth when the sample size is not too small. It seems that for small sample sizes, the fact that the convex hull of the data is rather small hinders the performance of the log-concave estimator, but that this effect is reduced as the sample size increases. The log-concave estimator copes well with the dependence in case (b), and it also deals particularly impressively with case (c), where the true density decays to zero at the boundary of the positive orthant.

In case (f), where the log-concavity assumption is violated, the performance of our estimator is not as good as the kernel estimate with the optimally chosen bandwidth, but is still comparable in most cases with the LSCV method. One would not expect the MISE of $\hat{f}_{n}$ to approach zero as $n\rightarrow\infty$ if log-concavity is violated, and in fact we conjecture that in this case the log-concave maximum likelihood estimator will converge to the density $f^{*}$ that minimises the Kullback–Leibler divergence $d(f_{0}\,\|\,f)=\int f_{0}(x)\log\frac{f_{0}(x)}{f(x)}\,dx$ over $f\in\mathcal{F}_{0}$. Such a result would be interesting for robustness purposes, because it could be interpreted as saying that provided the underlying density does not violate the log-concavity assumption too seriously, the log-concave maximum likelihood estimator is still sensible.

Clustering example

In a recent paper, Chang and Walther (2008) introduced an algorithm which combines the univariate log-concave maximum likelihood estimator with the EM algorithm (Dempster et al., 1977), to fit a finite mixture density of the form

Owing to the previous lack of an algorithm for computing the maximum likelihood estimator of a multidimensional log-concave density, Chang and Walther (2008) discuss an extension of the model in (6.1) to a multivariate context where the univariate marginal densities of each component in the mixture are assumed to be log-concave, and the dependence structure within each component density is modelled with a normal copula. Now that we are able to compute the maximum likelihood estimator of a multidimensional log-concave density, we can carry this method through to its natural conclusion. That is, in the finite mixture model (6.1) for a multidimensional log-concave density $f$, we simply assume that each of the component densities $f_{1},\ldots,f_{p}$ is log-concave. An interesting problem that we do not address here that of finding appropriate conditions under which this model is identifiable – see Titterington et al. (1985, Section 3.1) for a nice discussion.

2 Breast cancer example

We illustrate the log-concave EM algorithm on the Wisconsin breast cancer data set of Street et al. (1993), available on the UCI Machine Learning Repository website (Asuncion and Newman, 2007):

http://archive.ics.uci.edu/ml/datasets/Breast+Cancer+Wisconsin+%28Diagnostic%29.

The data set was created by taking measurements from a digitised image of a fine needle aspirate of a breast mass, for each of 569 individuals, with 357 benign and 212 malignant instances. We study the problem of trying to diagnose (cluster) the individuals based on the standard errors of two of the measurements, namely the radius of the cell nucleus (mean of distances from center to points on the perimeter, $X$) and its texture (standard deviation of grey-scale values, $Y$). The data are presented in Figure 4(a). In fact, the full data set consists of 30 measurements for each patient, representing the mean, standard error and ‘worst’ (mean of the three largest values) of 10 different features computed for each cell nucleus in the image. Since one would reasonably expect the means of each feature to be approximately normally distributed, and hence the Gaussian EM algorithm to be appropriate, we took the standard errors of the first two measurements to illustrate the log-concave EM algorithm methodology.

It is important also to note that although for this particular data set we do know whether a particular instance is benign or malignant, we did not use this information in fitting our mixture model. Instead this information was only used afterwards to assess the performance of the method, as reported below. Thus we are studying a clustering (or unsupervised learning) problem, by taking a classification (or supervised learning) data set and ‘covering up the labels’ until it comes to performance assessment.

The skewness in the data suggests that the mixture of Gaussians model may be inadequate, and in Figure 4(b) we show the contour plot and misclassified instances from this model. The corresponding plot obtained from the log-concave EM algorithm is given in Figure 4(c), while Figure 4(d) plots the fitted mixture distribution from the log-concave EM algorithm. For this example, the number of misclassified instances is reduced from 144 with the Gaussian EM algorithm to 121 with the log-concave EM algorithm.

In some examples, it will be necessary to estimate $p$, the number of mixture components. In the general context of model-based clustering, Fraley and Raftery (2002) cite several possible approaches for this purpose, including methods based on resampling (McLachlan and Basford, 1988) and an information criterion (Bozdogan, 1994). Further research will be needed to ascertain which of these methods is most appropriate in the context of log-concave component densities.

Plug-in estimation of functionals, sampling and the bootstrap

Suppose $X$ has density $f$. Often, we are less interested in estimating a density directly than in estimating some functional $\theta(f)$. Examples of functionals of interest (some of which were given in Section 1), include:

The differential entropy of $X$ (or $f$), defined by $H(f)=-\int f(x)\,\log f(x)\,dx$

For some functionals we can compute $\hat{\theta}=\theta(\hat{f}_{n})$ analytically. If this is not possible, but we can write $\theta(f)=\int f(x)g(x)\,dx$, we may approximate $\hat{\theta}$ by

for some (large) $B$, where $X_{1}^{*},\ldots,X^{*}_{B}$ are independent samples from $\hat{f}_{n}$. Conditional on $X_{1},\ldots,X_{n}$, the strong law of large numbers gives that $\hat{\theta}_{B}\stackrel{{\scriptstyle a.s.}}{{\rightarrow}}\hat{\theta}$ as $B\rightarrow\infty$. In practice, even when analytic calculation of $\hat{\theta}$ was possible, this method was found to be fast and accurate.

In order to use this Monte Carlo procedure, we must be able to sample from $\hat{f}_{n}$. Fortunately, this can be done efficiently using the following rejection sampling procedure. As in Section 4, for $j\in J$ let $A_{j}$ be the $d\times d$ matrix whose $l$th column is $X_{j_{l+1}}-X_{j_{1}}$ for $l=1,\ldots,d$, and let $\alpha_{j}=X_{j_{1}}$, so that $w\mapsto A_{j}w+\alpha_{j}$ maps the unit simplex $T_{d}$ to $C_{n,j}$. Recall that $\log\hat{f}_{n}(X_{i})=y_{i}^{*}$, and let $z_{j}=(z_{j,1},\ldots,z_{j,d})$, where $z_{j,l}=y_{j_{l+1}}^{*}-y_{j_{1}}^{*}$ for $l=1,\ldots,d$. Write

We may then draw an observation $X^{*}$ from $\hat{f}_{n}$ as follows:

Select $j^{*}\in J$, selecting $j^{*}=j$ with probability $q_{j}$

Select $w\sim\textrm{Unif}(T_{d})$ and $u\sim\textrm{Unif}()$ independently. If

accept the point and set $X^{*}=A_{j}w+\alpha_{j}$. Otherwise, repeat (ii).

2 Simulation study

In this section we illustrate some simple applications of this technique to functionals (c) and (d) above, using the Monte Carlo procedure and sampling scheme described in Section 7.1. Estimates are based on random samples from a $N_{2}(0,I)$ distribution, and we compare the performance of the LogConcDEAD estimate with that of a kernel-based plug-in estimate, where the bandwidth matrix was chosen using our knowledge of the underlying density to minimise the MISE.

For the differential entropy estimators, we find a similar pattern to that observed in Section 5: the log-concave plug-in estimator provides an improvement on the kernel-based estimator for the moderate and large sample sizes in our simulations. For the case of highest density regions, the relative performance of the log-concave estimator is better for the estimation of smaller density regions. In Figure 5, we illustrate the estimation of three highest density regions based on 500 points from a $N_{2}(0,I)$ distribution. For comparison, a kernel-based plug-in estimate (where the regions are not guaranteed to be convex) is also given.

In real data examples, we are unable to assess uncertainty in our functional estimates by taking repeated samples from the true underlying model. Nevertheless, the fact that we can sample from the log-concave maximum likelihood estimator does mean that we can apply standard bootstrap methodology to compute standard errors or confidence intervals, for example. Finally, we remark that the plug-in estimation procedure, sampling algorithm and bootstrap methodology extend in an obvious way to the case of a finite mixture of log-concave densities.

Concluding discussion

We have developed methodology that gives a fully automatic nonparametric density estimate under the condition that the density is log-concave, and shown how it may be extended to fit finite mixtures of log-concave densities. We have indicated a wide range of possible applications, including classification, clustering and functional estimation problems. The area of shape-constrained estimation is currently undergoing rapid growth, as evidenced by the many recent publications cited in the penultimate paragraph of Section 1, as well as recent workshops in Oberwolfach (November 2006), Eindhoven (October 2007) and Bristol (November 2007). We hope that this paper will stimulate further interest and research in the field.

As well as the continued development and refinement of the computational algorithms and graphical displays of estimates, and studies of theoretical performance, there remain many challenges and interesting directions for future research. These include:

Studying other shape constraints. These have received some attention for univariate data, dating back to Grenander (1956), but much less in the multivariate setting.

Developing both formal and informal diagnostic tools for assessing the validity of shape constraints.

Assessing the uncertainty in shape-constrained nonparametric density estimates, through confidence intervals/bands.

Developing analogous methodology for discrete data from shape-constrained distributions.

Examining nonparametric shape constraints in regression problems.

Studying methods for choosing the number of clusters in nonparametric, shape-constrained mixture models.

Appendix A Glossary of terms and results from convex analysis and computational geometry

which always exists (allowing $-\infty$ and $\infty$ as limits) provided $\sigma(y)$ is finite.

Appendix B Proofs

a product of log-concave functions. Thus $f_{X|P_{V}(X)}(\cdot|t)$ is log-concave for each $t$.

Thus $f$ is log-concave, as required. $\Box$

Completion of the Proof of Theorem 2 We prove each of the steps (i)–(v) outlined in Section 3 in turn. First note that if $x_{0}\in C_{n}$, then by Carathéodory’s theorem (Theorem 17.1 of Rockafellar (1997)), there exist distinct indices $i_{1},\ldots,i_{r}$ with $r\leq d+1$, such that $x_{0}=\sum_{l=1}^{r}\lambda_{l}X_{i_{l}}$ with each $\lambda_{l}>0$ and $\sum_{l=1}^{r}\lambda_{l}=1$. Thus, if $f(x_{0})=0$, then by Jensen’s inequality,

so $f(X_{i})=0$ for some $i$. But then $\psi_{n}(f)=-\infty$. This proves (i).

Note that $\log f$ has no direction of increase, because if $x\in C_{n}$, $z$ is a non-zero vector and $t>0$ is large enough that $x+tz\notin C_{n}$, then $-\infty=\log f(x+tz)<\log f(x)$. It follows by Theorem 27.2 of Rockafellar (1997) that the supremum of $f$ is finite (and is attained). Using properties (i) and (ii) as well, we may write $\int f(x)\,dx=c$, say, where $c\in(0,\infty)$. Thus $f(x)=c\bar{f}(x)$, for some $\bar{f}\in\mathcal{F}_{0}$. But then

with equality only if $c=1$. This proves (iv).

To prove (v), we may assume by (iv) that $\exp(\bar{h}_{y})$ is a density. Let $\max_{i}\bar{h}_{y}(X_{i})=M$ and let $\min_{i}\bar{h}_{y}(X_{i})=m$. We show that when $M$ is large, in order for $\exp(\bar{h}_{y})$ to be a density, $m$ must be negative with $|m|$ so large that \psi_{n}\bigl{(}\exp(\bar{h}_{y})\bigr{)}\leq\psi_{n}(f). First observe that if $x\in C_{n}$ and $\bar{h}_{y}(X_{i})=M$, then for $M$ sufficiently large we must have $M-m>1$, and then

For $\exp(\bar{h}_{y})$ to be a density, then, we require $m\leq-\frac{1}{2}e^{(M-1)/d}\mu(C_{n})^{1/d}$ when $M$ is large. But then

when $M$ is sufficiently large. This proves (v).

It is not hard to see that for any $M>0$, the function y\mapsto\psi_{n}(\exp(\bar{h}_{y})\bigr{)} is continuous on the compact set $[-M,M]^{n}$, and thus the proof of the existence of a maximum likelihood estimator is complete. To prove uniqueness, suppose that $f_{1},f_{2}\in\mathcal{F}$ and both $f_{1}$ and $f_{2}$ maximise $\psi_{n}(f)$. We may assume $f_{1},f_{2}\in\mathcal{F}_{0}$, $\log f_{1},\log f_{2}\in\mathcal{H}$ and $f_{1}$ and $f_{2}$ are supported on $C_{n}$. Then the normalised geometric mean

However, by Cauchy–Schwarz, $\int_{C_{n}}\{f_{1}(y)f_{2}(y)\}^{1/2}\,dy\leq 1$, so $\psi_{n}(g)\geq\psi_{n}(f_{1})$. Equality is obtained if and only if $f_{1}=f_{2}$ almost everywhere, but since $f_{1}$ and $f_{2}$ are continuous relative to $C_{n}$ (Theorem 10.2 of Rockafellar (1997)), this implies that $f_{1}=f_{2}$. An alternative way of proving the uniqueness of the maximum likelihood estimator may be based on the fact that \psi_{n}\bigl{(}tf_{1}+(1-t)f_{2}\bigr{)}>t\psi_{n}(f_{1})+(1-t)\psi_{n}(f_{2}) for all $t\in(0,1)$, provided $f_{1}$ and $f_{2}$ are distinct elements of $\mathcal{F}$. $\Box$

In this section, we find explicitly the set of points at which the function $\sigma$ defined in (4.1) is differentiable, and compute a subgradient of $\sigma$ at each point. For $i=1,\ldots,n$, define

Assume (A1). (a) For $y\in\mathcal{Y}$, the function $\sigma$ is differentiable at $y$ and for $i=1,\ldots,n$ satisfies

(b) For $y\in\mathcal{Y}^{c}$, the function $\sigma$ is not differentiable at $y$, but the vector $(\partial_{1}(y),\ldots,\partial_{n}(y))$ is a subgradient of $\sigma$ at $y$.

If $j_{1}=i$, then for sufficiently small $t$, we have $z_{j}^{(t)}=z_{j}-t1_{d}$, where $1_{d}$ denotes a $d$-vector of ones, so that $b^{(t)}_{j}=b_{j}-t(A_{j}^{T})^{-1}1_{d}$ and $\beta^{(t)}_{j}=\beta_{j}-t(1+\langle A_{j}^{-1}\alpha_{j},1_{d}\rangle)$

If $j_{l+1}=i$ for some $l\in\{1,\ldots,d\}$, then for sufficiently small $t$, we have $z_{j}^{(t)}=z_{j}+te_{l}^{d}$, so that $b^{(t)}_{j}=b_{j}+t(A_{j}^{T})^{-1}e_{l}^{d}$ and $\beta_{j}^{(t)}=\beta_{j}+t\langle A_{j}^{-1}\alpha_{j},e_{l}^{d}\rangle$.

where to obtain the final line we have made the substitution $x=A_{j}w+\alpha_{j}$, after taking the limit as $t\rightarrow 0$.

when $z_{1},\ldots,z_{d}$ are non-zero and distinct. In Cule et al. (2008b), it is shown that the required formula is

References