Abstract: 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 study develops a nonparametric methodology using Bernstein polynomials for bivariate density estimation where each of the marginal densities is restricted to unimodal.