Unimodal density estimation using Bernstein polynomials

The estimation of probability density functions is one of the fundamental aspects of any statistical inference. Many data analyses are based on an assumed family of parametric models, which are known to be unimodal (e.g., exponential family, etc.). Often a histogram suggests the unimodality of the underlying density function. Parametric assumptions, however, may not be adequate for many inferential problems. A flexible class of mixture of Beta densities that are constrained to be unimodal is presented. It is shown that the estimation of the mixing weights, and the number of mixing components, can be accomplished using a weighted least squares criteria subject to a set of linear inequality constraints. The mixing weights of the Beta mixture are efficiently computed using quadratic programming techniques. Three criteria for selecting the number of mixing weights are presented and compared in a small simulation study. More extensive simulation studies are conducted to demonstrate the performance of the density estimates in terms of popular functional norms (e.g., Lp norms). The true underlying densities are allowed to be unimodal symmetric and skewed, with finite, infinite or semi-finite supports. A code for an R function is provided which allows the user to input a data set and returns the estimated density, distribution, quantile, and random sample generating functions.