Semi-Parametric Estimation Using Bernstein Polynomial and a Finite Gaussian Mixture Model

The central focus of this paper is upon the alleviation of the boundary problem when the probability density function has a bounded support. Mixtures of beta densities have led to different methods of density estimation for data assumed to have compact support. Among these methods, we mention Bernstein polynomials which leads to an improvement of edge properties for the density function estimator. In this paper, we set forward a shrinkage method using the Bernstein polynomial and a finite Gaussian mixture model to construct a semi-parametric density estimator, which improves the approximation at the edges. Some asymptotic properties of the proposed approach are investigated, such as its probability convergence and its asymptotic normality. In order to evaluate the performance of the proposed estimator, a simulation study and some real data sets were carried out.


Introduction
Density estimation is a widely adopted tool for multiple tasks in statistical inference, machine learning, visualization and exploratory data analysis. Existing density estimation algorithms can be categorized into either parametric, semi-parametric, or non-parametric approaches. In the non-parametric framework, several methods have been set forward for the smooth estimation of density and distribution functions. The most popular one, called kernel method, was introduced by [1]. The advances were carried out by [2] to estimate a density function. The reader is recommended to consult the paper [3] for an introduction to several kernel smoothing techniques. However, kernel methods display estimation problems at the edges, when we have a random variable X with density function f supported on a compact interval. Moreover, if X 1 , . . . , X n is a sample with the same density f , it is well known, in non-parametric kernel density estimation, that the bias of the standard kernel density estimator is of a larger order near the boundary than that in the interior, where K is a kernel (that is, a positive function satisfying K(x)dx = 1) and (h n ) is a bandwidth (that is, a sequence of positive real numbers that goes to zero). Let us now suppose that f has two continuous derivatives everywhere and that, as n → ∞, h = h n → 0 and nh → 0. Let x = ph for p > 0. Near the boundary, the expression of the mean and the variance are indicated as and These bias phenomena are called boundary bias. Numerous authors have elaborated methods for reducing these phenomena, such as data reflection [4], boundary kernels [5][6][7], local linear estimator [8,9], use of beta and gamma kernels [10,11] and bias reduction [12,13]. For a smooth estimator of a density function f with finite known support, there have been several methods, such as Vitale's method [14], which is based on Bernstein polynomials and expressed asf 1,n,m where F n is the empirical distribution function and b k (m, x) = C k m x k (1 − x) m−k is the Bernstein polynomial. This estimator was investigated in the literatures [15][16][17][18] and, more recently, by [12,19,20].
Within the parametric framework, it is noteworthy that the Gaussian mixture model can be used to estimate any density function, without any problem of estimation on the edge. This refers to the fact that the set of all normal mixture densities is dense in the set of all density functions under the L 1 metric [21]. The investigation of mixture models stands for a full field in modern statistics. It is a probabilistic model introduced by [22] to illustrate the presence of subpopulations within an overall population. It has been developed so far by various authors, such as [23]. It is used for data classification and it provides efficient approaches of model-based clustering. The authors of [24] demonstrated that, when a Gaussian mixture model is used to estimate a density non-parametrically, the density estimator that uses the Bayesian information criterion (BIC) of [25] to select the number of components in the mixture is consistent [26].
However, we obtain the non-parametric kernel estimate of a density if we fit a mixture of n components in equal proportions 1/n, where n is the size of the observed sample. As a matter of fact, it can be inferred that mixture models occupy an interesting niche between parametric and non-parametric approaches to statistical estimation.
More recently, in the parametric context, [27] proposed a parametric model using Bernstein polynomials with positive coefficients to estimate the unknown density function f ; this estimator is defined as follows: where and p mi are the estimators of the parameters p mi , obtained by the Expectation Maximization (EM) algorithm as follows: mi , for i = 1, . . . , m. The proposed method gives a consistent estimator in L 2 distance under some conditions.
The problem at the edge does not arise for the parametric model. For this reason, the basic idea of this work is to consider a shrinkage method using Bernstein (Vitale's estimator) and Gaussian mixture estimators, to construct a shrinkage density estimator, in order to improve the approximation at the edge. A shrinkage estimator is a convex combination between estimators [28]. Basically, this implies that a naive or raw estimate is improved by combining it with other information.
The remainder of this paper is organized as follows: In the next section, we recall some intrinsic properties of the classical EM algorithm in the context of the Gaussian mixture parameter estimation. In Section 3, we introduce a new semi-parametric estimation approach based on the shrinkage method using Bernstein polynomials and Gaussian mixture densities. In Section 4, the consistency of the proposed estimator is exhibited, as well as its asymptotic normality.
Section 5 highlights a simulation study that compares the performance of the proposed approach with the Bernstein estimator, the standard Gaussian kernel estimator and Guan's estimator. The closing Section 6 crowns the whole work, wraps the conclusion and provides new perspectives for future work.

Background The Gaussian Mixture Model and Em Algorithm
Let us consider X = (X 1 , . . . , X n ), a sequence of independent and identically distributed (i.i.d.) with common Gaussian mixture density defined by where θ = (π, µ, σ) = (π 1 . . . , π K , µ 1 . . . , µ K , σ 1 , . . . , σ K ), Finally, for each observed data point X i , we associate a component label vector Z i in order to manage the data clustering. This random vector Z i = (Z ik ) 1≤k≤K is defined such that Z ik = 1 if the considered observation X i is drawn form the k kt component of the mixture and Z ik = 0 otherwise. Consequently, Z i is distributed as a multivariate Bernoulli distribution with vector parameters (π 1 , . . . , π K ) as follows: The EM algorithm is a popular tool in statistical estimation problems involving incomplete data or problems which can be posed in a similar form, such as the mixture parameters estimation [23,29]. In the EM framework, (X 1 , . . . , X n , Z 1 , . . . , Z n ) corresponds to the complete data and (Z 1 , . . . , Z n ) stand for the hidden data. Hence, the complete-data log-likelyhood is expressed by The two steps of the EM algorithm, after l iterations, are the following: (i) E-step: The conditional expectation of the complete-data log-likelyhood given the observed data, using the current fit θ (l) , is defined by ϕ θ|θ (l) = E θ (l) (L(X 1 , . . . , X n , Z 1 , . . . , Z n , θ)|X 1 , . . . , X n ).
The posterior probability that X i belongs to the jth component of the mixture at the lth iteration, is expressed as Finally, we obtain (ii) M-step: It consists of a global maximization of ϕ θ|θ (l) with respect to θ.

Proposed Approach
The proposed semi-parametric approach rests upon the shrinkage combination between the Gaussian mixture model and the Bernstein density estimators using the EM algorithm for the parameter estimations. The literature on shrinkage estimation is enormous. From this perspective, it is noteworthy to mention the most relevant contributions. The authors of [28] were the first to introduce the classic shrinkage estimator. The authors of [31] provided theory for the analysis of risk. Oman [32,33] developed estimators which shrink Gaussian density estimators towards linear subspaces. An in-depth investigation of shrinkage theory is displayed in Chapter 5 of [34].
The proposed semi-parametric approach based upon estimating the density function f relies on the same principle of Stein's works and there are two aspects along this line. The first setting is non-parametric in the sense that we do not assume any parametric form of the density. The non-parametric setting is very important as it allows us to perform statistical inference without making any assumption on the parametric form of the true density f . The second setting is to consider the Gaussian mixture model as a parametric estimator of the unknown density f .
In what follows, we consider X 1 , . . . , X n a sequence of i.i.d. random variables having a common unknown density function f supported on [0, 1]. We here develop a shrinkage method to estimate the density function, which is divided into the following three steps: Step 1 We consider the Bernstein estimator of the density function f , which is defined as Step 2 In view of (13), we consider the Gaussian mixture density as an estimator of the density function f , given bỹ where µ k , σ k and π k are estimated by the EM algorithm defined in (13).
Step 3 We consider the shrinkage density estimator f n,m form defined by and we use the EM algorithm to estimate the parameter λ ∈ [0, 1] of the proposed model.
By the same way as considered in Section 2, the two steps of the EM algorithm, after t iterations, are denoted in terms of the following:

1.
E-step: The conditional expectation of the complete-data log-likelihood given the observed data, using the current λ (t) , is provided by ) is a discrete random vector, following a multivariate Bernoulli distribution with vector parameters (λ, 1 − λ). Using Bayes's formula, we obtain the posterior probability in the tth iteration denoted by i1 .

2.
M-step: It consists of a global maximization of Q(λ|λ (t) ) with respect to λ.
The updated estimate of λ is indicated by The estimation of λ is obtained from by iterating the EM algorithm until convergence.
Therefore, the proposed estimator of the density function f is defined by Basically, it is a shrinkage estimator that shrinks the Bernstein estimator towards the Gaussian mixture density by a specified amount of λ. If λ = 1, the estimator f n,m reduces to the Bernstein estimatorf 1,n,n .

Convergence
In this section, we derive some asymptotic properties of the proposed estimator f n,m when the sample size tends to infinity. First, we assume that λ and K are fixed. The following proposition gives the probability convergence of the proposed estimator f n,m .
where L 2 denotes the mean quadratic convergence L 2 .
The proof of this lemma is reported in [35].
Proof of Proposition 1. First, using Lemma 1 and following the same steps as the proof of Second, based on Theorem 3.1 in [16], we obtaiñ In addition, referring to (18) and (19) and grounded on the application of Slutsky's Theorem, we conclude the proof.
According to [21], the density f (x) is a close approximation to the mixture density f 2 (x). Thus, the estimator f n,m (x) provides an approximation to the true density f (x).
To study the asymptotic normality of the estimator f n,m given by (17), we set forward the following assumptions in [36].
The following corollary is a consequence of the previous proposition which gives an asymptotic confidence interval of the density f , for a risk α ∈]0, 1[.

Corollary 1.
The 100(1 − α)% asymptotic confidence interval of f (x) is given by In the next section, we study the performance of the proposed estimator in estimating different distributions by comparing it to the performances of the Bernstein estimator and of the Gaussian kernel estimator.

Comparison Study
In this section, we investigate the performance of the proposed estimator given in (17), through estimating different densities by comparing it to the performance of the Bernstein estimator defined in (2), the standard Gaussian kernel estimator defined in (1) and the Guan's estimator defined in (3). We apply the Bernstein estimator when the sample is concentrated on the interval [0, 1]. For this purpose, we need to make some suitable transformations in the different cases that are listed as follows:

1.
Let us suppose that X is concentrated on a finite support [a, b]; then, we work with the sample values Y 1 , . . . , Y n , where Y i = (X i − a)/(b − a).

2.
For the density functions concentrated on R, we can use the transformed sample Y i = 1/2 + π −1 arctan(X i ), which transforms the range to the interval (0, 1).

3.
For the support R + , we can use the transformed sample Y i = X i /(1 + X i ), which transforms the range to the interval (0, 1). If the support is infinite, say, R = (−∞, ∞), we can consider [x 1 , x t ] ⊂ [a, b] as the finite support of f , where x 1 and x t are the minimum and the maximum order, respectively. We choose a and b such that F(a) and 1 − F(b) are of O(n −1 ), a < x 1 , and b > x t , where F is the distribution function [27]. Then, we can use the transformed sample, which transforms [x 1 , x t ] to the interval [0, 1] mentioned in the case 1.
In the simulation study, three sample sizes were considered, n = 50, n = 100, and n = 200, as well as the following density functions: (a) The beta mixture density 0.5B (3,9)  Our sample was decomposed into a learning sample of a size of 2/3 of the considered sample, on which the various statistical methods were constructed, and a second sample of a size of 1/3 of the considered sample, on which the predictive performance of the three methods were tested. For each density function f and sample size n, we computed the integrated squared error (ISE), the integrated absolute error (I AE) and the Kullback-Leibler divergence (KL) of the estimator f n,m over N = 500 trials.
where f k is the estimator computed from the kth sample and Indeed, it is advised to consider a learning sample bigger than a testing sample. In this work, our sample was decomposed into a learning sample of a size of 2/3 of the considered sample, on which the various statistical methods were constructed, and a second sample of a size of 1/3 of the considered sample, on which the predictive performances of the three methods were tested. Each run of the proposed estimator performed the following steps:

-
We first generated a random sample (X i ) 1≤i≤n of size n from the models' density (a) − ( f ). -We then split the generated data into a training set of a size of 2/3 of the considered sample and a test set of a size of 1/3 of the considered sample.

-
We applied the proposed estimator, using the observed data X i only from the training set, in order to estimate the density function.

-
The test set was then used to compute the estimation errors ISE, I AE and KL.
To select the optimal parameter K, we used the Gap Statistics algorithm [37]. We considered a Monte Carlo experiment to select the optimal choice of the degree m of the Bernstein polynomial and the bandwidth h of the kernel estimator, for each point x ∈ [0, 1].
We determined the parameters m (for 1 ≤ m ≤ 300) and h (for h = i/1000 with 1 ≤ i ≤ 300), which minimized the ISE, which was approximated by the ISE.
We considered N = 500 random samples of sizes n = 50, n = 100 and n = 200.  Departing from Tables 1-3 and Figure 1, we deduce the following:

-
The results displayed in Tables 1-3 show that the ISE, I AE and KL decreased as the sample size increased.

-
Using the proposed estimator, we obtained better results than those given by the other estimators in a large part of the cases. When we changed the parameters of the gamma mixture density in the sense that we had a smaller bias, our estimator was more competitive than the other approaches and we obtained better results.

-
Almost in all considered cases, the average KL of the density estimator (17) was smaller than that obtained by the Bernstein estimator defined in (2), that of the kernel estimator defined in (1) and that of Guan's estimator.

-
In the considered distribution 0.5B (3,9) + 0.5B (9,3), by choosing the appropriate m, the curve of the proposed distribution estimator (3.4) was closer to the true distribution than that of Guan's estimator (1.3), even when the sample size was very large.
Referring to Figures 2 and 3, we infer the following: -None of the estimators for the gamma mixture density 0.5G(6, 1) + 0.5G(1, 6) had good approximations near x = 0. However, the ISE of the proposed estimator was closer to zero than that of the Bernstein estimator and the kernel estimator, especially near the edge x = 1.

-
Guan's estimator and the kernel estimator for the normal mixture density 0.25N (2, 1) + 0.75N (−3, 1) had good approximations near x = 0. However, the ISE of the pro-posed estimator was closer to zero than that of the other estimators, especially near the two edges.
Therefore, we note that, for difficult distributions that diverge at the boundaries, the proposed method would fail, but not as badly as the standard methods without shrinkage. In addition, the performed simulations revealed that, on average, the proposed approach could lead to satisfactory estimates near the boundaries, better than the classical Bernstein estimator.

COVID-19 Data
In this subsection, we consider the COVID-19 data displayed in the INED website https://dc-covid.site.ined.fr/fr/donnees/france/ (accessed on 16 February 2022). These data concern the numbers of deaths due to COVID-19 in France (daily) from 21 March 2021, for 454 days. These data are such that min i (x i ) = 605 and max i (x i ) = 0. Then, it is convenient to assume that the density of the numbers of deaths is defined on the interval [0, 605] and transform the data into the interval unit. The Monte Carlo procedure was performed and resulted in h = 0.07659612 for the standard kernel estimator defined in (1.1), m 1 = 20 for the Bernstein estimator defined in (1.2), the proposed estimator, and m 2 = 12 for Guan's estimator. These estimators are exhibited in Figure 1 (right panel) along with a histogram of the data. All the estimators are smooth and seem to capture the pattern highlighted by the histogram. We record that the proposed estimator outperformed the other estimators near the boundaries.

Tuna Data
The last example concerns the tuna data reported in [38]. The data are derived from an aerial line transect survey of Southern Bluefin Tuna in the Great Australian Bight. An aircraft with two spotters on board flew randomly over allocated line transects. These data correspond to the perpendicular sighting distances (in miles) of 64 detected tuna schools to the transect lines. The survey was conducted in summer when tuna data tend to stay on the surface. The data are such that min i (x i ) = 0.19 and max i (x i ) = 16.26. The Monte Carlo procedure was performed and resulted in h = 0.1079 for the standard kernel estimator defined in (1), m 1 = 13 for the Bernstein estimator defined in (2) and the proposed estimator, and m 2 = 6 for Guan's estimator. These estimators are illustrated in Figure 4 (left panel) along with a histogram of the data. All the estimators are smooth and seem to capture the pattern highlighted by the histogram. We assert that the proposed estimator outperformed the other estimator, especially near the boundaries.

Conclusions
In this paper, we propose a shrinkage estimator of a density function based on the Bernstein density estimator and using a finite Gaussian mixture density. This method rests on three steps. The first step consists of considering the Bernstein estimatorf 1,n,m . The second relies upon the Gaussian Mixture densityf 2,n as an estimator of the unknown density f . The last step consists of considering the shrinkage form λf 1,n,m + (1 − λ)f 2,n and EM algorithm in order to estimate the parameter λ. The asymptotic properties of this estimator were established. Afterwards, we demonstrate the effectiveness of the proposed method using some simulated and real data. We clarify how it can lead to very satisfactory estimates near the boundaries and in terms of ISE, I AE and KL. Eventually, we would simply assert that our research work is a step that may be taken further, extended and built upon as it lays the ground and paves the way for future works to elaborate a semiparametric regression estimator using the shrinkage method. We also plan to work on the case where λ is a random variable. Another future research direction would be to extend our findings to the setting of serially dependent observations.