Quasi-concave density estimation

Introduction.

Our objective is to introduce a general class of shape constraints applicable to the estimation of probability densities, multivariate as well as univariate. Elements of the class are represented by restricting certain monotone functions of the density to lie in convex cones. Maximum likelihood estimation of log-concave densities constitutes an important special case; however, the wider class allows us to include a variety of other shapes. A one parameter subclass modeled on the means of order $\rho$ studied by Hardy, Littlewood and Pólya 1934 incorporates all the quasi-concave densities, that is, all densities with convex upper contour sets. Estimation methods for these densities, as described below, bring new opportunities for statistical data analysis.

Log-concave densities play a crucial role in a wide variety of probabilistic models: in reliability theory, search models, social choice and a broad range of other contexts it has proven convenient to assume densities whose logarithm is concave. Recognition of the importance of log-concavity was already apparent in the work of Schoenberg and Karlin on total positivity beginning in the late 1940s. Karlin 1968 forged a link between log-concavity and classical statistical properties such as the monotone likelihood ratio property, the theory of sufficient statistics and uniformly most powerful tests. Maximum likelihood estimation of densities constrained to be log-concave has recently enjoyed a considerable vogue with important contributions of Walther (Walther 2001, Walther 2002, Walther 2009), Pal, Woodroofe and Meyer 2007, Rufibach 2007, Dümbgen and Rufibach 2009, Balabdaoui, Rufibach and Wellner 2009, Chang and Walther 2007 and Cule, Samworth and Stewart 2010, among others.

Log-concave densities are constrained to exhibit exponential tail behavior. This restriction motivates a search for weaker forms of the concavity constraint capable of admitting common densities with algebraic tails like the $t$ and $F$ families. The $\rho$-concave densities introduced in Section 2 constitute a rich source of candidates. While it would be possible, in principle, to consider maximum likelihood estimation of such densities, duality considerations lead us to consider a more general class of maximum entropy criteria. Maximizing Shannon entropy in the dual is equivalent to maximum likelihood for the leading log-concave case, but other entropies are also of interest. Section 3 describes several examples arising in the dual from the class of Rényi entropies, each corresponding to a distinct specification of the concavity constraint, and each corresponding to a distinct fidelity criterion in the primal. The crucial advantage of adapting the fidelity criterion to the form of the concavity constraint is that it assures a convex optimization problem with a tractable computational strategy.

Quasi-concave probability densities and their estimation.

A probability density function, $f$, is called log-concave if $-{\log f}$ is a (proper) convex function on the support of $f$. We adhere to the usual conventions of Rockafellar 1970, which allow convex functions to take infinite values—although we will allow only $+\infty$, because all our convex functions will be proper. The domain of a convex (concave) function, $\operatorname{dom}g$, is then the set of $x$ such that $g(x)$ is finite. We adopt the convention $-{\log 0}=+\infty$.

Unimodality of concave functions implies that log-concave densities are unimodal. An interesting connection in the multivariate case was pointed out by Silverman 1981: the number of modes of a kernel density estimate is monotone in the bandwidth when the kernel is log-concave. However, as illustrated by the Student $t$ family, not every unimodal density is log-concave. Laplace densities, with their exponential tail behavior, are; but heavier, algebraic tails are ruled out. This prohibition motivates a relaxation of the log-concavity requirement.

for $\rho\neq 0$; the limiting case for $\rho=0$ is

In this terminology, $\log$-concave functions are $0$-concave and concave functions are $1$-concave. As $M_{\rho}(a,p)$ is monotone, increasing in $\rho$ for $a\geq 0$ and any $p\in\mathcal{S}$, it follows that if $f$ is $\rho$-concave, then $f$ is also $\rho^{\prime}$-concave for any $\rho^{\prime}<\rho$. Thus, concave functions are $\log$-concave, but not vice-versa. In the limit $-\infty$, concave functions satisfy the condition

so they are (and consequently for all $\rho$-concave functions) quasi-concave.

The hierarchy of $\rho$-concave density functions was considered in the economics literature by Caplin and Nalebuff 1991 in spatial models of voting and imperfect competition; their results reveal some intriguing connections to Tukey’s half-space depth in multivariate statistics; see Mizera 2002. Curiously, it appears that the first thorough investigation of the mathematical concept of quasi-concavity was carried out by de Finetti 1949. Further details and motivation for $\rho$-concave densities can be found in Prékopa 1973, Borell 1975 and Dharmadhikari and Joag-Dev 1988.

2 Maximum likelihood estimation of log-concave densities.

It is convenient to recast (1) in terms of $g=-{\log f}$, the estimate becoming $f={e}^{-g}$,

The objective function of (2) is equal to $+\infty$, given the convention adopted above, unless all $X_{i}$ are in the domain of $g$. As in Silverman 1982, it proves convenient to move the integral constraint into the objective function,

a device that ensures that the solution integrates to one without enforcing this condition explicitly. Apart from the multiplier $1/n$, the crucial difference between (2) and (3) is that the latter is a convex problem, while the former not.

It is well known that naïve, unrestricted maximum likelihood estimation is doomed to fail when applied in the general density estimation context: once “log-concave” is dropped from the formulation of (1), any sequence of putative maximizers is attracted to the the linear combination of point masses situated at the data points. One escape from this “Dirac catastrophe” involves regularization by introducing a roughness, or complexity, penalty; various proposals in this vein can be found in Good 1971, Silverman 1982, Gu 2002 and Koenker and Mizera 2008.

Another way to obtain a well-posed problem is by imposing shape constraints, a line of development dating back to the celebrated Grenander 1956 nonparametric maximum likelihood estimator for monotone densities. While monotonicity regularizes the maximum likelihood estimator, unimodality per se—somewhat surprisingly—does not. The desired effect is achieved only by enforcing somewhat more stringent shape constraints—for instance log-concavity, sometimes also called “strong unimodality.” An advantage of shape constraints over regularization based on norm penalties is that it is not encumbered by the need to select additional tuning parameters; on the other hand, it is limited in scope—applicable only when the shape constraint is plausible for the unknown density.

3 Quasi-concave density estimation.

Expanding the scope of our investigation, we now replace ${e}^{-g}$ in the integral of the objective function by a generic function $\psi(g)$ and define

The following conditions on the form of $\psi$ will be imposed:

The domain of $\psi$ is an open interval containing $(0,+\infty)$.

The limit, as $\tau\to+\infty$, of $\psi(y+\tau x)/\tau$ is $+\infty$ for every real $y$ and any $x<0$.

$\psi$ is differentiable on the interior of its domain.

$\psi$ is bounded from below by $0$, with $\psi(x)\to 0$ when $x\to+\infty$.

The most crucial condition is (A1) ensuring the convexity of $\Phi$. Condition (A2) assures that $\psi(x)$ is finite for all $x>0$, while (A3) is required in the proof of the existence of the estimates. The relationship between primal and dual formulations of the estimation problem is facilitated by (A4), and (A5) rules out possible complications regarding the existence of the integral $\int\psi(g)\,dx$ in (4), allowing for the convention $\psi(+\infty)=0$. In the spirit of the Lebesgue integration theory, the integral then exists, although $\psi(g)$ does not have to be summable: it is either finite [which is automatically true for any $g$ convex and $\psi(g)={e}^{-g}$] or $+\infty$. In the latter case, the objective function $\Phi(g)$ is considered to be equal to $+\infty$; $\Phi(g)$ is also $+\infty$ if $g(X_{i})=+\infty$ for some $X_{i}$, which occurs unless all $X_{i}$ lie in the domain of $g$. On the other hand, any $g$ equal to some positive constant on $\mathcal{H}(X)$ and $+\infty$ elsewhere yields $\Phi(g)<\infty$.

A rigorous treatment without assumption (A5), that is, for functions $\psi$ not bounded below, would introduce technicalities involving handling of the integrals in the spirit of singular integrals of calculus of variations, a strategy resembling the contrivance of Huber 1967 of subtracting a fixed quantity from the objective function to ensure finiteness of the integral. Although we do not believe that such formal complications are unsurmountable, we do not pursue such a development.

Careful deliberation reveals that replacing $g$ by its closure (lower semicontinuous hull) does not change the integral term in (4), and potentially only decreases the first term; this means that without any restriction of its scope, we may reformulate the estimation problem as

Unlike in (3), $\psi(g)$ is not necessarily the estimated density $f$; the relationship of $g$ to $f$ will be revealed in Section 3, together with the motivation leading to concrete instances of some possible functions $\psi$.

4 Characterization of estimates.

Any function of this type is finitely generated in the sense of Rockafellar 1970, whose Corollary 19.1.2 asserts that it is polyhedral, being the maximum of finitely many affine functions, and therefore convex. The convention $\inf\varnothing=+\infty$ used in (6) means that the domain of $g_{(X,Y)}$ is equal to $\mathcal{H}(X)$. If $h$ is a convex function such that $h(X_{i})\leq Y_{i}$, for all $i$, then $h(x)\leq g_{(X,Y)}(x)$ for all $x$; the function $g_{(X,Y)}$ is thus the maximum of convex functions with this property—the lower convex hull of points $(X_{i},Y_{i})$.

The theorem shows that it is sufficient to seek potential solutions of (5) in $\mathcal{G}(X)$; this means, due to the one–one correspondence of the latter to $\mathcal{D}(X)$, that the optimization task (5) is essentially finite dimensional. The theorem also justifies the transition to a more convenient optimization domain in the primal formulation appearing in the next section.

Duality, entropy and divergences.

The conjugate dual formulation of the primal estimation problem (5) conveys a maximum entropy interpretation and leads us to several concrete proposals for $\psi$. To conform to existing mathematical apparatus, we begin by further clarifying the optimization and constraint functional classes of our primal formulation. For definitions and general background on convex analysis, our primary references are Rockafellar (Rockafellar 1970, Rockafellar 1974) and Zeidler 1985; we may also mention Hiriart-Urruty and Lemaréchal 1993 and Borwein and Lewis 2006.

Hereafter, $\mathcal{K}(X)$ will denote the cone of closed (lower semicontinuous) convex functions on $\mathcal{H}(X)$, the convex hull of $X$. This cone is a subset of $\mathcal{C}(X)$, the collection of functions continuous on $\mathcal{H}(X)$; it is important that $\mathcal{C}(X)$ is a linear topological space, with respect to the topology of uniform convergence. Note that $\mathcal{G}(X)\subset\mathcal{K}(X)\subset\mathcal{C}(X)$. In view of Theorem 2.1, any solution of (5) is also the solution of

and conversely; thus, we will refer to (7) as our primal formulation.

2 The dual formulation.

Since $\psi$ is nonincreasing, there are no affine functions with positive slope that minorize the graph of $\psi$, hence $\psi^{\ast}(y)=+\infty$ for all $y>0$. If $\psi$ is differentiable on the (nonempty) interior of its domain, then $\psi^{\ast}$ can be obtained using differential calculus—as the Legendre transformation of $\psi$; denoting the derivative $\psi^{\prime}$ by $\chi$, we have

where $\chi^{-1}(y)$ is any solution, $z$, of the equation $\chi(z)=y$. The (topological) dual of $\mathcal{C}(X)$ is $\mathcal{C}^{\ast}(X)$, the space of (signed) Radon measures on $\mathcal{H}(X)$; its distinguished element is $P_{n}$, the empirical measure supported by the data points $X_{i}$. The polar cone to $\mathcal{K}(X)$ is

Suppose that assumptions (A1) and (A2) hold. The strong (Fenchel) dual of the primal formulation (7) is

in the sense that the value, $\Phi(g)$, of the primal objective for any $g$ satisfying the constraints of (7), dominates the value, for any $f$ satisfying the constraints of (3.1), of the objective function in (3.1); the minimal value of (7) and maximal value of (3.1) coincide. Moreover, there exists $f$ attaining the maximal value of (3.1). Any dual feasible function $f$, that is, any $f$ satisfying the constraints of (3.1) and yielding finite objective function of (3.1), is a probability density with respect to the Lebesgue measure: $f\geq 0$ and $\int f\,dx=1$. If condition (A4) is also satisfied, then the dual and primal optimal solutions satisfy the relationship $f=-\psi^{\prime}(g)$.

It should be emphasized that the expression of absolute continuity in (3.1) is a requirement on $F=P_{n}-G$; the dual objective function is defined as the conjugate to the primal objective function $\Phi$, and is equal to $-\infty$ for any Radon measure that is not absolutely continuous with respect to the Lebesgue measure. This is how regularization operates here: only those $F$ qualify for which $P_{n}$ gets canceled with the discrete component of $G$. Once $F$ satisfies this requirement, its density integrates to $1$, as shown in the proof of Theorem 3.1. The nonnegativity for $f$ yielding finite dual objective function is the consequence of $\psi^{\ast}(-y)$ being infinite for $y<0$. In practical implementations, it may be prudent to enforce $f\geq 0$ in the dual explicitly as a feasibility constraint.

3 The interpretation of the dual.

An immediate consequence of Theorem 3.1 is that we can reformulate the maximum likelihood problem posed in (3) as an equivalent maximum (Shannon) entropy problem.

Maximum likelihood estimation of a log-concave density as posed in (3) has an equivalent dual formulation

whose solution satisfies the relationship $f={e}^{-g}$, where $g$ is the solution of (3). In particular, the solution of (3) satisfies $\int{e}^{-g(x)}\,dx=1$, therefore problems (2) and (3) are equivalent.

The emergence of the Shannon entropy is hardly surprising—in view of the well-established connections of maximum likelihood estimation to the Kullback–Leibler divergence and maximum entropy. Note that the dual criterion can be also interpreted as choosing the $f$ closest in the Kullback–Leibler divergence to the uniform distribution on $\mathcal{H}(X)$, from all $f$ satisfying the dual constraints.

4 Rényi entropies.

While the outcome of Corollary 3.1, the equivalence of (2) and (3), could be also shown by elementary means, it is important to emphasize that the real value of the dual connection lies in the vista of new possibilities it opens. To explore the link to potential alternatives, we consider the family of entropies originally introduced for $\alpha>0$ by Rényi (Rényi 1961, Rényi 1965),

as an extension of the limiting case for $\alpha=1$, the Shannon entropy. For $\alpha\neq 1$, maximizing (13) over $f$ is equivalent to the maximization of

The dependence of convexity/concavity properties of $y^{\alpha}$ necessitates a separate treatment of the cases with $\alpha>1$, when the conjugate pair is

where $\beta$ and $\alpha$ are conjugates in the usual sense that $1/\beta+1/\alpha=1$. See Figure 1.

The general form of the primal formulations (7) corresponding to (13) can be written, for $\alpha\neq 1$, in a unified way as

together with the relation between the dual and primal solutions, $f=|g|^{\beta-1}$. Several particular instances merit special attention.

5 Power divergences.

For $\alpha>1$, we may write $(-g)$ instead of $|g|$, and then introduce $h=-g$. The resulting primal formulation is

By Theorem 2.1, this formulation is equivalent to

After substituting $f^{1/(\beta-1)}$ for $h$, multiplying by $\beta$, and rewriting in terms of $\alpha$ we obtain a new objective function

which recalls the “minimum density power divergence estimators,” proposed, for $\alpha\geq 1$, by Basu et al. 1998 in the context of estimation in parametric families.

6 Pearson χ2\chi^{2}.

Although $\alpha=2$ is a special case of the power divergence family mentioned above, it deserves a special mention. The choice of $\alpha=2$ in the Rényi family leads to the dual formulation

The primal formulation can be written, after the application of Theorem 2.1, in a particularly simple form

which can be interpreted as a variant of the minimum Pearson $\chi^{2}$ criterion. A similar theme can be found in the dual, which can be interpreted as returning among all densities satisfying its constraints the one with minimal Pearson $\chi^{2}$ distance to the uniform density on $\mathcal{C}(X)$.

7 Hellinger.

While the form of the objective function for $\alpha=2$ has some computational advantages, its secondary consequence—constraining the density itself to be concave rather than its logarithm—is not at all appealing. Indeed, all Rényi choices with $\alpha>1$ impose a more restrictive form of concavity than log-concavity. From our perspective, it seems more reasonable to focus attention on weaker forms of concavity, corresponding to $\alpha\leq 1$. Apart from the celebrated log-concave case $\alpha=1$, a promising alternative would seem to be Rényi entropy with $\alpha=1/2$. This choice in the Rényi system leads to the dual

and primal, again after the application of Theorem 2.1,

The estimated density satisfies $f=1/g^{2}$, which means that the primal constraint, $g\in\mathcal{K}$, enforces the convexity of $g=1/\sqrt{f}$. In the terminology of Section 2, the estimated density is now required to be only $-1/2$-concave, a significant relaxation of the log-concavity constraint; in addition to all log-concave densities, all the Student $t_{\nu}$ densities with $\nu\geq 1$ satisfy this requirement. The dual problem (21) can be interpreted as a Hellinger fidelity criterion, selecting from the cone of dual feasible densities the one closest in Hellinger distance to the uniform distribution on $\mathcal{H}(X)$.

8 The frontier and beyond?

Although the original Rényi system was confined to $\alpha>0$, a limiting form for $\alpha=0$ can be obtained similarly to the $\alpha=1$ case. It yields the conjugate pair

As is apparent from Figure 1, this $\psi$ violates our condition (A5), but may nevertheless deserve a brief consideration. Note first that the possible complications with existence of integrals may occur only in the formulation (5) with unbounded domain—not in (7), where all integrals are of bounded functions over a compact domain. The major technical complications with $\psi$ violating (A5) concern theorems in Section 4, and are briefly discussed there. Here we mention only that the resulting dual, adapted directly from (3.1), is

In this case $g=1/f$, and the estimate is constrained to be $-1$-concave, a yet still weaker requirement that admits all of the Student $t_{\nu}$ densities for $\nu>0$.

If we interpret the dual problem (3.1), for $\alpha=1$, as choosing a constrained $f$ to minimize the Kullback–Leibler divergence of $f$ from the uniform distribution on $\mathcal{H}(X)$, we can similarly interpret the $\alpha=0$ dual as minimizing the reversed Kullback–Leibler divergence. In parametric estimation, the latter objective is sometimes associated with empirical likelihood, while the former is associated with exponentially tilted empirical likelihood. See, for example, Hall and Presnell 1999 for related discussion in the context of kernel density estimation, and Schennach 2007.

One might try to continue in this fashion marching inexorably toward weaker and weaker concavity requirements. There appears to be no obstacle in considering $\alpha<0$; the general form (15) of the primal is still applicable. The shape constraints corresponding to negative $\alpha$ encompass a wider and wider class of quasi-concave densities, eventually arriving at the $-\infty$-concave constraint, at which point we would have sanctioned all of the quasi-concave densities. But formal complications, as well as computational difficulties dictate the more prudent strategy of restricting attention to $\alpha>0$ cases.

Existence and Fisher consistency of estimates.

Returning to our general setting, existence, uniqueness and Fisher consistency are established under mild conditions on the function $\psi$.

Theorem 2.1 not only justifies the choice of the optimization domain in (7), but also shows, due to the one–one correspondence between $\mathcal{G}(X)$ and $\mathcal{D}(X)$, that the optimization task (7) is essentially finite dimensional, parametrized by the values $Y_{i}=g(X_{i})$. This facilitates the proof of the following existence result.

Suppose that assumptions (A1), (A2), (A3) and (A5) hold, and that $\mathcal{H}(X)$ has a nonempty interior. Then the formulation (5) has a solution $g\in\mathcal{C}(X)$; if $\psi$ is strictly convex, then this solution is unique.

2 Fisher consistency.

While anything else in this direction may be viewed as speculative, Fisher consistency, a crucial prerequisite for a more detailed asymptotic theory, can be verified in a quite straightforward manner and essentially complete generality. For differentiable $\psi$, Theorem 3.1 gives the relationship between the solution $g$ of the optimization task (7) and the density estimate: $f=-\psi^{\prime}(g)$. Using the notation $\chi$ for $\psi^{\prime}$, and $\chi^{-1}$ for its inverse, as in Section 3, we may write $g=\chi^{-1}(-f)$, and subsequently rewrite the formulation (5) in terms of the estimated density $f$ (omitting, for brevity, the integration variables)

This yields a new objective function—which we nevertheless denote, slightly abusing the notation, also $\Phi$. The population version of this $\Phi$ is obtained by replacing $dP_{n}$ by $f_{0}\,dx$:

The Fisher consistency for an estimator defined by solving (4.2) requires that $\Phi_{0}(f_{0})\leq\Phi_{0}(f)$, for every $f$; however, there may be a formal problem now with the existence of the integral in (25), as $\chi^{-1}(f)$ may take both positive or negative values. A possible way of handling this obstacle is the strategy of Huber 1967, briefly mentioned in Section 2: instead of $\Phi$, we consider a modified objective function

is now better suited for the ensuing version of the Fisher consistency theorem.

In fact, Theorem 4.2 can be proved in the same manner for the unmodified $\Phi$, if $\operatorname{dom}\psi=(\omega,+\infty)$ with $\omega>-\infty$. Then the inverse of $\chi=\psi^{\prime}$, and hence the range of $\chi^{-1}$ is bounded from below by $\omega$. In such a case, $\chi(f)f_{0}\geq\omega f_{0}$, so the first term in (25) is minorized by an integrable function $\omega f_{0}$; the second term is bounded from below by $0$ by assumption (A5), so the whole integral then exists in the Lebesgue sense, being either finite or equal to $+\infty$.

If, however, assumption (A5) is not satisfied, then the existence of the integral should be assumed explicitly; we return to this point briefly at the end of the proof of Theorem 4.2. Note that, by comparing (8) and (25), existence of the integral is equivalent to assuming the integrability (summability) of

that is, the existence and finiteness of the entropy term in the dual (3.1).

Examples of practical use.

We employed two independent algorithms for solving the convex programming problems posed above: mskscopt from the MOSEK software package of Andersen 2006, and the PDCO MATLAB procedure of Saunders 2003. Both algorithms are coded in MATLAB and employ similar primal-dual, log-barrier methods. Further details regarding numerical implementation appear in Appendix B. The crux of both algorithms is a sequence of Newton-type steps that involve solving large, very sparse least squares problems, a task that is very efficiently carried out by modern variants of Cholesky decomposition. Several other approaches have been explored for computing quasi-concave density estimators that are log-concave. An active set algorithm for univariate log-concave density estimation was described in Dümbgen, Hüsler and Rufibach 2007 and implemented in the R package logcondens of Rufibach and Dümbgen 2009. Cule, Samworth and Stewart 2010 have recently implemented a promising steepest descent algorithm for multivariate log-concave estimation that may be adaptable to other quasi-concave density estimation problems.

To illustrate the application of the foregoing methods, we briefly consider some realistic examples. Our first example features data similar to those considered by Pal, Woodroofe and Meyer 2007, the type of data where shape constraints sometimes arise in a natural manner. The two samples consist of 9092 measurements of radial and 3933 of rotational velocity for the stars from Bright Star Catalog, Hoffleit and Warren 1991. The left and right panels of Figure 2 show the results for the radial and rotational velocity samples, respectively.

The broken line in the upper panels shows kernel density estimates, each time with default MATLAB bandwith selection; the solid lines correspond to one of the norm penalized estimates proposed in Koenker and Mizera 2008: maximum likelihood penalized by the total variation of the second derivative of the logarithm of the estimated density. This is the $L^{1}$ version of the Silverman’s (Silverman 1982) estimator penalizing the squared $L^{2}$ norm of the third derivative. The smoothing parameter for the latter estimate was set quite arbitrarily at $1$; it seems that this arbitrary choice works quite satisfactorily here, providing—somewhat surprisingly, for both samples—about the same level of smoothing as the kernel estimator. For the radial velocity sample, the two estimates are essentially the same. For the rotational velocity sample, however, the right upper panel shows that the kernel density estimate differs substantially from the penalized one. Both estimators have the unfortunate effect of assigning considerable mass to negative values despite the fact that there are no negative observations. This effect is somewhat more pronounced for the kernel estimate.

Since the preliminary analyses of the upper panels indicates that the hypothesis of unimodality is plausible for both of the datasets, a natural next step is the application of a shape-constrained estimator—a move that, among other things, may relieve us of insecurities related to the arbitrary choice of smoothing parameters. The broken line in the lower panels of Figure 2 shows the log-concave maximum likelihood ($\alpha=1$), and the solid line the Hellinger ($-1/2$-concave) estimate ($\alpha=1/2$). While, as expected, there is almost no difference between the two (and in fact, among all four) estimates for the radial velocity dataset, the right lower panel reveals that the log-concave estimate yields for the rotational velocity sample a density that is monotonically decreasing—which contradicts the evidence suggested by all other methods. The Hellinger estimate, on the other hand, exhibits a subtle, but visible bump at the location of the plausible mode, thus turning out to be visually more informative about the center of the data than the tails. This is somewhat paradoxical given its original heavy-tail motivation, confirming that the real universe of data analysis can be much more subtle than that of the surrounding theoretical constructs.

2 Bivariate example: Criminal fingers.

To illustrate our approach in a simple bivariate setting, we reconsidered the well-known MacDonell 1902 data on the heights and left middle finger lengths of 3000 British criminals. This data was employed by Gosset in preliminary simulation work described in Student 1908.

Figure 3 illustrates the Hellinger ($-1/2$-concave) fit of this data. Contours are labeled in units of log-density. A notable feature of the data is the unusual observation in the middle of the upper edge. This point is highly anomalous, at least for any density with exponential tail behavior. The maximum likelihood estimate of the log-concave model in Figure 4 has very similar central contours, but the outer contours fall off much more rapidly implying that the log-concave estimate assigns much smaller probability to the region near the unusual point.

3 Some simulation evidence.

Motivated by a suggestion of one of the referees, we undertook some numerical experiments to explore performance of our shape constrained estimators and evaluate whether consistency appeared to be a plausible conjecture. For the log-concave estimator Pal, Woodroofe and Meyer (Pal, Woodroofe and Meyer 2007) report “Hellinger error” for a fully crossed design involving five target densities and five sample sizes with 500 replications per cell.

In Figure 5, we report results of our attempt to reproduce the PWM experiment expanded somewhat to consider two competing estimators: the adaptive kernel estimator of Silverman (Silverman 1986) using a Gaussian kernel, and the logspline estimator of Kooperberg and Stone (Kooperberg and Stone 1991) as implemented in the logspline R package of Kooperberg. Five target densities are considered: (standard) normal, Laplace, Gamma(3), Beta($3,2$) and Weibull($3,1$) as in PWM. Five sample sizes are studied: 50, 100, 200, 500, 1000. And two measures of performance are considered: squared Hellinger distance as in PWM in the left panel and $L_{1}$ distance in the right panel. Plotted points in these figures represent cell means. Both figures support the contention that the rates of convergence are comparable for all three estimators.

Figure 6 reports results from a similar experiment for the the $-1/2$-concave estimator described in Section 3.6. We consider five new target densities: lognormal, $t_{3}$, $t_{6}$, $F_{3,6}$ and Pareto(5), all of which fall into the $-1/2$-concave class. The same competing estimators and sample sizes are used. In a small fraction of cases for the second group of densities, less than 0.2 percent, there were problems either with the convergence of the logspline or shape-constrained estimator, or with the numerical integration required to evaluate the performance measures, so Figure 6 plots cell medians rather than cell means. Again, the figures support the conjecture that the rates of convergence for the shape-constrained estimator are competitive with those of the adaptive kernel and logspline estimators.

A concise way to summarize results from these experiments is to estimate the simple linear model

where $y_{ij}$ denotes a cell average of our two error criteria for one of our three estimators, for target density $i$ and sample size $n_{j}$. In this rather naïve framework, $\hat{\beta}$ can be interpreted as an empirical rate of convergence for the estimator. In Tables 1 and 2, we report these estimates suppressing the estimated target density

specific $\alpha_{i}$’s. In this comparison too, the shape constrained estimators perform quite well.

Extensions and conclusions.

We have described a rather general approach to qualitatively constrained density estimation. Log-concave densities are an important target class, but other, weaker, concavity requirements that permit algebraic tail behavior are also of considerable practical interest. Ultimately, the approach accommodates all quasi-concave densities as a limit of the Rényi entropy family.

There are many unexplored directions for future research. As we have seen, a consequence of the variational formulation of our concavity constraints is that the estimated densities vanish off the convex hull of the data. Various treatments for this malady may be suggested. Müller and Rufibach 2009 have recently suggested applying one of several estimators of the Pareto tail index to the smoothed ordinates from the log-concave preliminary density estimator. Our inclination would be to prefer solutions that impose further regularization on the initial problem. Thus, for example, we can add a new penalty term to the primal problem, penalizing the total variation of the derivative (gradient) of $\log f$, and choosing a suitable value of the associated Lagrange multiplier to smooth the tail behavior at the boundary.

We have adhered, thus far, to the principle that the entropy choice in the fidelity criterion of the dual problem should dictate the form of the convexity constraint: likelihood thus implies log-concavity, Hellinger fidelity implies $1/\sqrt{f}$ concavity, etc. One may wish to break this linkage and consider maximum likelihood estimation of general $\rho$-concave densities. This may have some advantages from an inference viewpoint, at the cost of complicating the numerical implementation.

Embedding shape constrained density estimation of the type considered here into semiparametric methods would seem to be an attractive option in many settings. And, it would obviously be useful to consider inference for the validity of shape constraints in the larger context of penalized density estimation. We hope to pursue some of these issues in future work.

Appendix A Proofs

[Proof of Theorem 2.1] Given $h$ convex, put $Y_{i}=h(X_{i})$ and take $g=g_{(X,Y)}$, the function defined by (6). The convexity of $h$ implies that $h(x)\leq g(x)$ for every $x$; since $\psi$ is nonincreasing, we have

The definition of $g_{(X,Y)}$ implies that $h(X_{i})=g_{(X,Y)}(X_{i})=Y_{i}$; therefore, the rest of $\Phi$ remains unchanged, and (29) implies that $\Phi(g)\leq\Phi(h)$.

the resulting function is convex itself, being a sup of affine functions. For any $f\in\mathcal{C}(X)$ and any Radon measure $G$, a linear functional from $\mathcal{C}(X)^{\ast}$, we have

where $\Phi$ is the original objective function of (7) and $\Upsilon=\iota_{\mathcal{K}(X)}$ is the indicator function of $\mathcal{K}(X)$. The expression for the Fenchel dual of this type of problem follows from Rockafellar 1966; see also Rockafellar 1974, Section 5, Example 11:

(note that one of the conjugates, in both cites sources, is in the “concave” sense, which explains the negative sign of the argument in the second term, but not in the first). The conjugate of the indicator of a convex cone $\mathcal{K}(X)$ is the indicator of $-\mathcal{K}(X)^{-}$ [Rockafellar 1974, Section 3, equation (3.14)]. The term $-\Upsilon^{\ast}(-G)$ in the objective can be therefore interpreted as a constraint $-G\in-\mathcal{K}(X)^{-}$, that is, $G\in\mathcal{K}(X)^{-}$. The definition of the conjugate of $\Phi$ gives

and $\Psi^{\ast}$ is its conjugate. The form of the latter is given by Rockafellar 1971, Corollary 4A: if $G$ is absolutely continuous with respect to the Lebesgue measure, then

otherwise $\Psi^{\ast}(G)=+\infty$. These facts, and expressions (A) and (32), yield (3.1).

Rockafellar [(Rockafellar 1966), Theorem 1; see also Rockafellar (Rockafellar 1974), Section 8, Example 11′] gives also a constraint qualification for this type of problem: to prove strong duality, we need to find some $g$ where both $\Phi$ and $\Upsilon$ are finite and one of them is continuous. Such a $g$ is provided by a function constant on $\mathcal{H}(X)$, say $g(x)=1$ for all $x\in\mathcal{H}(X)$. It is convex, thus $\Upsilon(g)=0$ is finite. So is $\Phi(g)$; the topology on $\mathcal{C}(X)$ is that of uniform convergence, and $\psi$ is continuous at $1$, hence there is a neighborhood of $g$ containing only functions for which $\Phi$ is finite and $\Phi$ is continuous at $g$.

Once the constraint qualification is verified, we know that the primal and dual optimal values coincide (zero duality gap), and that the dual is attained—there is an optimal solution to the dual; see Theorem 52.B(3) of Zeidler 1985. Due to the fact that $\psi$ is decreasing, $\psi^{\ast}(-f)=+\infty$ whenever $f<0$; thus, if $f$ yields a finite dual objective function, then $f$ is nonnegative. If $G\in\mathcal{K}(X)^{-}$, then $\langle G,f\rangle\leq 0$ for every $f\in\mathcal{K}(X)$; consequently,

Therefore, $\langle G,1\rangle=0$ and for every dual feasible $f$,

That is, every dual feasible $f$ is a probability density with respect to the Lebesgue measure.

If a primal solution, $g$, exists—the fact that is established by Theorem 4.1, but here we are exploring only the consequences of such a premise—it is related to the dual solution, $f$, via extremality (Karush–Kuhn–Tucker) conditions. The form of this relationship asserted by the theorem follows from the second condition of (8.24) in Rockafellar [(Rockafellar 1974), Section 8, Example 11′], together with the form of the subgradient of $\Psi$ given by Rockafellar [(Rockafellar 1971), Corollary 4B], combined with the fact established above that the estimated density $f$ corresponds to $F=P_{n}-G$. {proof}[Proof of Theorem 4.1] By Theorem 2.1, any potential solution of (5) lies within the class $\mathcal{G}(X)$ of polyhedral functions; due to the one–one correspondence between $\mathcal{G}(X)$ and $\mathcal{D}(X)$, the set of vectors discretely convex relative to $X$, the existence proof can be carried for (7) reparametrized by $Y_{i}$, the putative values of $g(X_{i})$. In what follows, $X$ remains fixed, and $\alpha$, $\beta$ will denote generic coefficients of convex combinations: any real numbers satisfying $\alpha,\beta\geq 0$, $\alpha+\beta=1$.

As $\psi$ is nonincreasing and convex, we obtain

as was required. Note that the integral is also finite whenever $Y$ has all components in the domain of $\psi$, due to the polyhedral character of $g_{(X,Y)}(x)$ and the fact that $\psi$ is nonincreasing and $\mathcal{H}(X)$ is bounded. Otherwise, it may be equal only to $+\infty$; hence $\Psi_{\mathcal{D}}$ is a proper convex function.

Note first that by the definition, $g_{(X,Y)}(x)=y$ identically on $\mathcal{H}(X)$ for constant $Y$; likewise, for every $x\in\mathcal{H}(X)$,

Choose a real number $y$ lying in the domain of $\psi$, and set $Y=(y,\ldots,y)$. According to Lemma A.1, $g_{(X,Y)}(x)=y$ for every $x\in\mathcal{H}(X)$; the function constant on $\mathcal{H}(X)$ is convex, hence $Y\in\mathcal{D}(X)$. Then

when $\tau\to+\infty$. For the integral over $E^{0}$, the limit is obviously zero. If $\psi$ satisfies assumptions (A1) and (A5), then $\psi$ is nonincreasing and converging to $0$ when $\tau\to+\infty$; for every $x\in E^{+}$ then $\psi(y+\tau g(x))/\tau$ monotonically decreases with increasing $\tau$, hence the limit of the integral over $E^{+}$ is zero as well. Finally, if $\psi$ satisfies also (A3), then for every $x\in E^{-}$ the limit of $\psi(y+\tau g(x))/\tau$ is $+\infty$; at the same time, the expression is bounded from below by $0$, due to (A4). The application of the Fatou lemma then gives that the limit of the integral over $E^{-}$ is $+\infty$, whenever $E^{+}$ has positive Lebesgue measure.

The proof of the theorem is then finished by the examination of possible alternatives. If the first term in (A), the mean of the $Z_{i}$’s, is positive, then the theorem is proved, as all other terms in (A) converge either to $0$ or $+\infty$. If the first term of (A) is negative or equal to zero, then there must be some $Z_{i}<0$ (the case with all $Z_{i}=0$ is excluded). That means that $g(x)=g_{(X,Z)}(x)$ is negative for some $x$; this implies that $E^{-}$ has positive Lebesgue measure, and then the limit of the second term in (A), the integral over $E^{-}$, and thus of the whole expression (A) is $+\infty$. This proves the theorem.

Under the strict convexity of $\psi$, the strict convexity of $\Phi_{\mathcal{D}}$ follows (for appropriate $\alpha$, $\beta$) from the second inequality in (A), which becomes sharp—this is due to the fact that the sharp inequality holds pointwise for all $x$, and thus for the integrals as well. The strict convexity of $\Phi_{\mathcal{D}}$ then implies the uniqueness of the solution.

Finally, functions $\psi$ satisfying (A1)–(A3), but not necessarily (A5) require some special considerations. For the integral over $E^{-}$, we have to observe that for every $\tau>0$ we have $\psi(y+\tau g_{(X,Z)}(x))\geq\psi$; if $\psi(y)<0$, then $\psi(y)/\tau\geq\psi(y)$ for every $\tau\geq 1$, if $\psi(y)\geq 0$ then $\psi(y)/\tau\geq 0$ for every $\tau>0$. Hence we have also in this case an integrable constant minorant [due to the fact that $\mathcal{H}(X)$ has finite Lebesgue measure]; this justifies the limit transition via the Fatou lemma. Finally, for the integral over $E^{+}$, we need to assume that the limit of the integrand is $0$ for every $x$, and find an integrable minorant; this may be related to the existence of the integral of the second term in (5). {proof}[Proof of Theorem 4.2] The proof relies on the application of what is called Fenchel’s inequality by Rockafellar [(Rockafellar 1970), page 105] or the (generalized) Young inequality by Hardy, Littlewood and Pólya [(Hardy, Littlewood and Pólya 1934), Section 4.8], or Zeidler [(Zeidler 1985), Section 51.1]. The inequality says says that for arbitrary $x,y$ and convex function $\psi$

Applied pointwise to $x=-f_{0}$ and $y=\chi^{-1}(-f)$, the inequality yields

which is equivalent to the nonnegativity of the integrand in (27). For $f=f_{0}$, the equality (28) implies that

For functions not satisfying (A5), integrability of $f_{0}\chi^{-1}(-f_{0})$ is no longer equivalent to that of $-\psi^{\ast}(-f_{0})$. However, if we assume the integrability of the latter, then the proof can be carried through in the same way.

Appendix B Computational details

Our computational objective is to provide a unified algorithmic strategy for solving the entire class of problems described above. Interior point methods designed for general convex programming and capable of exploiting the inherently sparse structure of the resulting linear algebra offer a powerful, general approach. We have employed two such implementations throughout our development process: the PDCO algorithm of Saunders 2003, and the MOSEK implementation of Andersen 2006.

Our generic primal problem (5) involves minimizing an objective function consisting of a linear component, representing likelihood or some generalized notion of fidelity, plus a nonlinear component, representing the integrability constraint. Minimization is then subject to a cone constraint imposing convexity. We will first describe our procedure for enforcing convexity, and then turn to the integrability constraint.

In dimension one convexity of piecewise linear functions can be imposed easily by enforcing linear inequality constraints on a set of function values, $\gamma_{i}=g(\xi_{i})$ at selected points $\xi_{1},\xi_{2},\ldots,\xi_{m}$. For ordered $\xi_{i}$’s, the convex cone constraint can be written as $D\gamma\geq 0$ for a tridiagonal matrix $D$ that does second differencing, adapted to the possible unequal spacing of the $\xi_{i}$’s.

A superior choice, one that circumvents the difficulties of the finite-element, fixed triangulation approach, relies on finite differences. Convexity is imposed directly at points on a regular rectangular grid using finite-differences to compute the discrete Hessian:

Convexity is then enforced by imposing positive semidefiniteness. These constraints—convexity at each of the grid points $(\xi_{1i},\xi_{2i})$—produce a semi-definite programming problem. In the bivariate setting the semi-definiteness of each $H$ can be reformulated as a rotated quadratic cone constraint; we need only constrain the signs of the diagonal elements of $H$ and its determinant. This simplifies the implementation of the Hellinger estimator in MOSEK. For the relatively fine grid used for Figure 3 solution requires about 25 seconds, considerably quicker than the log-concave estimate of Figure 4 computed with the implementation of Cule, Gramacy and Samworth 2009.

B.2 The integrability constraint.

For certain special $\psi$, one can evaluate the integral term $\int\psi(g(x))\,dx$ in the objective function of (5) explicitly—as was done for $\psi(g)=e^{-g}$ by Cule, Samworth and Stewart 2010. While such a strategy may also be possible for certain other specific $\psi$, we adopt a more pragmatic approach based on a straightforward Riemannian approximation

Here, $s_{i}$ are weights derived from the configuration of $\xi_{i}$. Of course, with only a modest number of $\xi_{i}$’s such an approximation may be poor; in dimension one we therefore augment the initial collection of the observed data $X_{1},\ldots,X_{n}$ by filling the gaps between their order statistics by further grid points, to ensure that the resulting grid (not necessarily uniformly spaced) and consisting of the observed data points as well as the new grid points, provides a sufficiently accurate approximation (36). The $s_{i}$’s are then simply the averages of the adjacent spacings between the ordered $\xi_{i}$’s. Given the size of problems modern optimization software can successfully handle, it is no problem to add an abundance of new points in dimension one.

In dimension two, the approximation (36) is based on the uniformly spaced grid of the points used in the finite-difference approach described in the previous subsection. As the original data points $X_{i}$ may no longer lie among the grid points $\xi_{i}$, we have to modify the fidelity component of the objective function: instead of obtaining $g(X_{i})$ directly, we obtain it via linear interpolation from the values of $g$ at the vertices of the rectangles enclosing $X_{i}$. As long as the grid is sufficiently fine, the difference is minimal. We use this approach often also in dimension one, as it provides better numerical stability especially for fine grids and large data sets.

B.3 Discrete duality.

Adopting the procedures described above, we can write the finite-dimensional version of the primal problem as

where $\Psi(\gamma)$ denotes now the $m$-vector with typical element $\psi(g(\xi_{i}))=\psi(\gamma_{i})$, $L$ is an “evaluation operator” which either selects the data elements from $\gamma$, or performs the appropriate linear interpolation from the neighboring ones, so that $L\gamma$ denotes the $n$-vector with typical element, $g(X_{i})$, and $w$ is an $n$-vector of observation weights, typically $w_{i}\equiv 1/n$.

Associated with the primal problem ( P ) is the dual problem

Here, $\eta$ is an $m$-vector of dual variables and $\phi$ is an $m$-vector of function values representing the density evaluated at the $\xi_{i}$’s, and $S=\operatorname{diag}(s)$. The vector $\Psi^{\ast}$ is the convex conjugate of $\Psi$ defined coordinate-wise with typical element $\psi^{\ast}(y)=\sup_{x}\{yx-\psi(x)\}$. Problems ( P ) and ( D ) are strongly dual in the sense of the following result, which may viewed as the discrete counterpart of Theorem 3.1.

If $\psi$ is convex and differentiable on the interior $\mathcal{I}$ of its domain, then the corresponding solutions of ( P ) and ( D ) satisfy

whenever the elements of g are from $\mathcal{I}$ and the elements of $f$ are from the image of $\mathcal{I}$ under $\psi^{\prime}$.

For $\Psi(x)$ with typical element $\psi(x)=e^{-x}$ we have $\Psi^{\ast}$ with elements $\psi^{\ast}(y)=-y\log y+y$, so the dual problem corresponding to maximum likelihood can be interpreted as maximizing the Shannon entropy of the estimated density subject to the constraints appearing in ( D ). Since $g$ was interpreted in ( P ) as $\log f$, this result justifies our interpretation of solutions of ( D ) as densities provided that they satisfy our integrability condition. This is easily verified and thus justifies the implicit Lagrange multiplier of one on the integrability constraint in ( P ), giving a discrete counterpart of Theorem 3.1.

Let $\iota$ denote an $m$-vector of ones, and suppose in ( P ) that $w^{\top}L\iota=1$ and $D\iota=0$. Then solutions $\phi$ of ( D ) satisfy $s^{\top}\phi=1$ and $\phi\geq 0$.

The crucial element of the proof is that the differencing operator D annihilates the constant vector and therefore the result extends immediately to other norm-type penalties as well as to the other entropy objectives that we have discussed. Indeed, since the second difference operator representing our convexity constraint annihilates any affine function it follows by the same argument that the mean of the estimated density also coincides with the sample mean of the observed $X_{i}$’s.

Acknowledgments.

We are grateful to Lutz Dümbgen, Kaspar Rufibach, Guenther Walther and Jon Wellner for sending us preprints of their work, to the referees for their very constructive comments, and to Mu Lin for help with the bright star example and Fisher consistency proof.

References