Bernstein Polynomial Model for Grouped Continuous Data

In real world applications of statistics, many data are provided in the form of frequencies of observations in some fixed mutually exclusive intervals, which are called grouped data. Strictly speaking, all the data from a population with a continuous distribution are grouped due to rounding of the individual observations (Hall, 1982). The EM algorithm has been used to deal with grouped data (Dempster et al., 1977). McLachlan & Jones (1988) introduced the EM algorithm for fitting mixture model to grouped data (see Jones & McLachlan, 1990, also). Under a parametric model, let $f(x;\theta)$ be the probability density function (PDF) of the underlying distribution with an unknown parameter $\theta$. The maximum likelihood estimate (MLE) of the parameter $\theta$ can be obtained from grouped data and is shown to be consistent and asymptotically normal (see, for example, Lindley, 1950; Tallis, 1967). Parametric MLE is sensitive to model misspecification and outliers. The minimum Hellinger distance estimate (MHDE) of the parameter using grouped continuous data is both robust for contaminated data and asymptotically efficient (Beran, 1977a, b). Parametric methods for grouped data requires evaluating integrals which makes the computation expensive. To lower the computation cost Lin & He (2006) proposed the approximate minimum Hellinger distance estimate (AMHDE) for grouped data by the data truncation and replacing the probabilities of class intervals with the first order Taylor expansion. Clearly their idea works for MLE based on grouped data.

Under nonparametric setting, the underlying PDF $f$ is unspecified. Based on grouped data $f$ can be estimated by the empirical density, the relative frequency distribution, which is actually a discrete probability mass function. The kernel density estimation (Rosenblatt, 1956, 1971) can be applied to grouped data (see Linton & Whang, 2002; Jang & Loh, 2010; Minoiu & Reddy, 2014, for example). The effects of rounding, truncating, and grouping of the data on the kernel density estimate have been studied, maybe among others, by Hall (1982), Scott & Sheather (1985), and Titterington (1983). However, the expectation of kernel density estimate is the convolution of $f$ and the kernel scaled by the bandwidth. It is crucial and difficult to select an appropriate bandwidth to balance between the bias and variance. Many authors have proposed different methods for data-based bandwidth selection over the years. The readers are referred to a survey by Jones et al. (1996) for details and more references therein. Another drawback of the kernel density is its boundary effect. Methods of boundary-effect correction have been studied, among others, by Rice (1984) and Jones (1993).

“All models are wrong”(Box, 1976). So all parametric models are subject to model misspecification. The normal model is approximate because of the central limit theorem. The goodness-of-fit tests and other methods for selecting a parametric model introduce additional errors to the statistical inference.

Any continuous function can be approximated by polynomials. Vitale (1975) proposed to estimate the PDF $f$ by estimating the coefficients $f(i/m)$ of the Bernstein polynomial (Bernstein, 1912) $\mathbb{B}f(t)=\sum_{i=0}^{m}f(i/m){m\choose i}t^{i}(1-t)^{m-i}$ by $\hat{f}(i/m)=(m+1)\{F_{n}[(i+1)/(m+1)]-F_{n}[i/(m+1)]\}$, $i=0,1,\ldots,m$, where $F_{n}$ is the empirical distribution function of $x_{1},\ldots,x_{n}$. Since then, many authors have applied the Bernstein polynomial in statistics in similar ways (see Guan, 2014, for more references). These and the kernel methods are not model-based and not maximum likelihood method. Thus they are not efficient. The estimated Bernstein polynomial $\widehat{\mathbb{B}f(t)}=\sum_{i=0}^{m}\hat{f}(i/m){m\choose i}t^{i}(1-t)^{m-i}$ aims at $\mathbb{B}f(t)$. It is known that the best convergence rate of $\mathbb{B}f(t)$ to $f(t)$ is at most ${\cal O}(m^{-1})$ if $f$ has continuous second or even higher order derivatives on [0,1]. Buckland (1992) proposed a density estimation with polynomials using grouped and ungrouped data with the help of some specified parametric models.

Thanks to a result of Lorentz (1963) there exists a Bernstein (type) polynomials $f_{m}(t;\bm{p})\equiv\sum_{i=0}^{m}p_{mi}\beta_{mi}(t)$, where $p_{mi}\geqslant 0$, $\beta_{mi}(t)=(m+1){m\choose i}t^{i}(1-t)^{m-i}$, $i=0,\ldots,m$, whose rate of convergence to $f(t)$ is at least ${\cal O}(m^{-r/2})$ if $f$ has a continuous $r$-th derivative on [0,1] and $r\geqslant 2$. This is called a polynomial with “positive coefficients” in the literature of polynomial approximation. Guan (2014) introduced the Bernstein polynomial model $f_{m}(t;\bm{p})$ as a globally valid approximate parametric model of any underlying continuous density function with support $[0,1]$ and proposed a change-point method for selecting an optimal degree $m$. It has been shown that the rate of convergence to zero for the mean integrated squared error(MISE) of the maximum likelihood estimate of the density could be nearly parametric, ${\cal O}(n^{-1+\epsilon})$, for all $\epsilon>0$. This method does not suffer from the boundary effect.

If the support of $f$ is different from [0,1] or even infinite, then we can choose an appropriate (truncation) interval $[a,b]$ so that $\int_{a}^{b}f(x)dx\approx 1$ (see Guan, 2014). Therefore, we can treat $[a,b]$ as the support of $f$ and we can use the linearly transformed data $y_{i}=(x_{i}-a)/(b-a)$ in $[0,1]$ to obtain estimate $\hat{g}$ of the PDF $g$ of $y_{i}$’s, respectively. Then we estimate $f$ by $\hat{f}(x)=\hat{g}\{(x-a)/(b-a)\}/(b-a)$. In this paper, we will assume that the density $f$ has support $[0,1]$.

This Bernstein polynomial model $f_{m}(t;\bm{p})$ is a finite mixture of the beta densities $\beta_{mi}(t)$ of beta$(i+1,m-i+1)$, $i=0,\ldots,m$, with mixture proportions $\bm{p}=(p_{m0},\ldots,p_{mm})$. It has been shown that the Bernstein polynomial model can be used to fit a ungrouped dataset and has the advantages of smoothness, robustness, and efficiency over the traditional methods such as the empirical distribution and the kernel density estimate (Guan, 2014). Because these beta densities and their integrals are specified and free of unknown parameters, this structure of $f_{m}(t;\bm{p})$ is convenient. It allows the grouped data to be approximately modeled by a mixture of $m+1$ specific discrete distributions. So the infinite dimensional “parameter” $f$ is approximately described by a finite dimensional parameter $\bm{p}$. This and the nonparametric likelihood are similar in the sense that the underlying distribution function is approximated by a step function with jumps as parameters at the observations.

Due to the closeness of $f_{m}(t;\bm{p})$ to $f(t)$, by the acceptance-rejection argument for generating pseudorandom numbers, almost all the observations in a sample from $f(t)$ can be used as if they were from $f_{m}(t;\bm{p})$. It will be shown in this paper that the maximizer of the likelihood based on the approximate model $f_{m}(t;\bm{p})$ targets $\bm{p}_{0}$ which makes $f_{m}(t;\bm{p}_{0})$ the unique best approximation of $f$. This acceptance-rejection argument can be used to prove other asymptotic results under an approximate model assumption.

In this paper we shall study the asymptotic properties of the Bernstein polynomial density estimate based on grouped data and ungrouped raw data as a special case of grouping. A stronger result than that of Guan (2014) about the rate of convergence of the proposed density estimate based on ungrouped raw data will be proved using a different argument. We shall also compare the proposed estimate with those existing methods such as the kernel density, parametric MLE, and the MHDE via simulation study.

The paper is organized as follows. The Bernstein polynomial model for grouped data is introduced and is proved to be nested in Section 2. The EM algorithm for finding the approximate maximum likelihood estimates of the mixture proportions is derived in this section. Some asymptotic results about the convergence rate of the proposed density estimate are given in Section 3. The methods for determining a lower bound for the model degree $m$ based on estimated mean and variance and for choosing the optimal degree $m$ are described in Section 4. In Section 5, the proposed methods are compared with some existing competitors through Monte Carlo experiments, and illustrated by the Chicken Embryo Data. The proofs of the theorems are relegated to the Appendix.

In this section we shall state results about the convergence rate of the density estimates which will be proved in the Appendix. Unlike most asymptotic results about maximum likelihood method which assume exact parametric models, we will show our results under the approximate model $f_{m}(t;\bm{p})=\sum_{j=0}^{m}p_{mj}\beta_{mj}(t)$. For a given $\bm{p}_{0}$, we define the norm

The squared distance between $\bm{p}$ and $\bm{p}_{0}$ with respect to norm $\|\cdot\|_{B}$ is

With the aid of the acceptance-rejection argument for generating pseudorandom numbers in the Monte Carlo method we have the following lemma which may be of independent interest.

For density estimation based on the raw data we have the following result.

Note that (11) is a stronger result than (12) which is an almost parametric rate of convergence for MISE. Guan (2014) showed a similar result under another set of conditions. The best parametric rate is ${\cal O}(n^{-1})$ that can be attained by the parametric density estimate under some regularity conditions.

For $\theta_{mi}(\bm{p})=\sum_{j=0}^{m}p_{mj}\{\mathcal{B}_{mj}(t_{i})-\mathcal{B}_{mj}(t_{i-1})\}$, we define norm

The squared distance between $\bm{p}$ and $\bm{p}_{0}$ with respect to norm $\|\cdot\|_{G}$ is

By the mean value theorem, we have

where $\Delta t_{i}=t_{i}-t_{i-1}$ and $t_{mij}^{*}\in[t_{i-1},t_{i}]$. Thus $\|\bm{p}-\bm{p}_{0}\|_{G}$ is a Riemann sum

where

For grouped data we have the following.

For the relationship between the norms $\|\bm{p}-\bm{p}_{0}\|_{B}^{2}$ and $\|\bm{p}-\bm{p}_{0}\|_{G}^{2}$, we have the following result.

For a grouped data based estimate $\tilde{\bm{p}}_{G}$, the rate of convergence of $\|\tilde{\bm{p}}_{G}-\bm{p}_{0}\|_{G}^{2}$ to zero is ${\cal O}((\log n)^{2}/n)$. However the rate of convergence of $\|\tilde{\bm{p}}_{G}-\bm{p}_{0}\|^{2}_{B}$ to zero depends on that of $\max_{i}\Delta t_{i}$. For equal-width classes, $\Delta t_{i}=1/N$, and $N=n^{1/2+2/k}$, we have $\|\tilde{\bm{p}}_{G}-\bm{p}_{0}\|_{B}^{2}={\cal O}((\log n)^{2}/n)+{\cal O}(m^{4}/n^{1+4/k})$. Thus $\|\tilde{\bm{p}}_{G}-\bm{p}_{0}\|_{B}^{2}={\cal O}((\log n)^{2}/n)$ if $m={\cal O}(n^{1/k})$. If $k$ is large, then $N\approx\sqrt{n}$.

Guan (2014) showed that the model degree $m$ is bounded below approximately by $m_{b}=\max\left\{1,\left\lceil\mu(1-\mu)/\sigma^{2}-3\right\rceil\right\}$. Based on the grouped data, the lower bound $m_{b}$ can be estimated by $\tilde{m}_{b}=\max\left\{1,\left\lceil\tilde{\mu}(1-\tilde{\mu})/\tilde{\sigma}^{2}-3\right\rceil\right\}$, where

Due to overfitting the model degree $m$ cannot be arbitrarily large. With the estimated $\tilde{m}_{b}$, we choose a proper set of nonnegative consecutive integers, $M=\{m_{0},m_{0}+1,\ldots,m_{0}+k\}$ such that $m_{0}<\tilde{m}_{b}$. Then we can estimate an optimal degree $m$ using the method of change-point estimation as proposed by Guan (2014). For each $m_{i}=m_{0}+i$ we use the EM algorithm to find the MBLE $\tilde{\bm{p}}_{m_{i}}$ and calculate $\ell_{i}=\ell(\tilde{\bm{p}}_{m_{i}})$. Let $y_{i}=\ell_{i}-\ell_{i-1}$, $i=1,\ldots,k$. The $y_{i}$’s are nonnegative because the Bernstein polynomial models are nested. Guan (2014) suggested that $y_{1},\ldots,y_{\tau}$ be treated as exponentials with mean $\mu_{1}$ and $y_{\tau+1},\ldots,y_{k}$ be treated as exponentials with mean $\mu_{0}$, where $\mu_{1}>\mu_{0}$, so that $\tau$ is a change point and $m_{\tau}$ is the optimal degree and use the change-point detection method (see Section 1.4 of Csörgő & Horváth, 1997) for exponential model to find a change-point estimate $\hat{\tau}$. Then we estimate the optimal $m$ by $\hat{m}=m_{\hat{\tau}}$. Specifically, $\hat{\tau}=\arg\max_{1\leqslant\tau\leqslant k}\{R(\tau)\}$, where the likelihood ratio of $\tau$ is