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. This paper presents a flexible class of mixture of Beta densities that are constrained to be unimodal. We show 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. We efficiently compute the number of mixing components and associated mixing weights of the beta mixture using quadratic programming techniques. Three criterion 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., L p norms). The true underlying densities are allowed to be unimodal symmetric and skewed, with finite, infinite or semi-finite supports. 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.