Recursive density estimators based on Robbins-Monro's scheme and using Bernstein polynomials

TL;DR

This paper introduces a recursive density estimator using Robbins-Monro's algorithm and Bernstein polynomials, addressing boundary issues and demonstrating its effectiveness through theoretical analysis and empirical validation.

Contribution

It presents a novel recursive density estimator that alleviates boundary problems, with proven asymptotic properties and validated performance against existing methods.

Findings

The estimator effectively handles boundary bias.

Theoretical asymptotic properties are established.

Simulation and real data confirm improved performance.

Abstract

In this paper, we consider the alleviation of the boundary problem when the probability density function has bounded support. We apply Robbins-Monro's algorithm and Bernstein polynomials to construct a recursive density estimator. We study the asymptotic properties of the proposed recursive estimator. We then compared our proposed recursive estimator with many others estimators. Finally, we confirm our theoretical result through a simulation study and then using two real datasets.

Tables

Table 1. Table 1: The average integrated squared error ( I ​ S ​ E 𝐼 𝑆 𝐸 ISE ) of Vitale’s estimator f ~ n subscript ~ 𝑓 𝑛 \widetilde{f}_{n} and the three recursive estimators; recursive 1 correspond to the estimator  f n subscript 𝑓 𝑛 f_{n} with the choice ( γ n ) = ( n − 1 ) subscript 𝛾 𝑛 superscript 𝑛 1 (\gamma_{n})=(n^{-1}) , recursive 2 correspond to the estimator  f n subscript 𝑓 𝑛 f_{n} with the choice ( γ n ) = ( [ 1 − a 2 ] ​ n − 1 ) subscript 𝛾 𝑛 delimited-[] 1 𝑎 2 superscript 𝑛 1 (\gamma_{n})=\left(\left[1-\frac{a}{2}\right]n^{-1}\right) ( a = 2 / 9 𝑎 2 9 a=2/9 ) and recursive 3 correspond to the estimator  f n subscript 𝑓 𝑛 f_{n} with the choice ( γ n ) = ( [ 1 − a ] ​ n − 1 ) subscript 𝛾 𝑛 delimited-[] 1 𝑎 superscript 𝑛 1 (\gamma_{n})=\left(\left[1-a\right]n^{-1}\right) ( a = 1 / 5 𝑎 1 5 a=1/5 ).

Density function	n	Vitale’s estimator	Recursive 1	Recursive 2	Recursive 3
$\left(a\right)$	50	0.085389	0.088252	0.089350	0.092738
	200	0.028167	0.025737	0.026057	0.027045
	500	0.013533	0.011398	0.011540	0.011977
$\left(b\right)$	50	0.145441	0.129385	0.130996	0.135963
	200	0.047977	0.037733	0.038202	0.039651
	500	0.023050	0.016710	0.016918	0.017560
$\left(c\right)$	50	0.074634	0.061163	0.061925	0.064273
	200	0.024620	0.017837	0.018059	0.0187441
	500	0.011828	0.007899	0.007997	0.0083012
$\left(d\right)$	50	0.108664	0.119492	0.120980	0.125567
	200	0.035845	0.034847	0.035281	0.036619
	500	0.017222	0.015433	0.015625	0.016217
$\left(e\right)$	50	0.124816	0.157438	0.159398	0.165442
	200	0.041174	0.045914	0.046485	0.048248
	500	0.019782	0.020333	0.020587	0.021367
$\left(f\right)$	50	0.084119	0.109237	0.110596	0.114790
	200	0.027749	0.031857	0.032253	0.033476
	500	0.013332	0.014108	0.014284	0.014825
$\left(g\right)$	50	0.063962	0.065611	0.066428	0.068947
	200	0.021099	0.019134	0.019372	0.020107
	500	0.010137	0.008474	0.008579	0.008904
$\left(h\right)$	50	0.046147	0.042890	0.043424	0.04507
	200	0.015222	0.012508	0.012664	0.013144
	500	0.007313	0.005539	0.005608	0.005821
$\left(i\right)$	50	0.031391	0.036611	0.037066	0.038472
	200	0.010355	0.010676	0.010809	0.011219
	500	0.004975	0.004728	0.004787	0.004968
$\left(j\right)$	50	0.071621	0.065671	0.066489	0.069010
	200	0.023626	0.019152	0.019390	0.020125
	500	0.011351	0.008481	0.008587	0.008913

Table 2. Table 2: The average integrated squared error ( I ​ S ​ E 𝐼 𝑆 𝐸 ISE ) of Leblanc estimator’s f ~ n , m , m / 2 subscript ~ 𝑓 𝑛 𝑚 𝑚 2 \widetilde{f}_{n,m,m/2} and the three estimators introduced by Kakizawa: f ~ n , m , m / b subscript ~ 𝑓 𝑛 𝑚 𝑚 𝑏 \widetilde{f}_{n,m,m/b} with b = 3 , 4 𝑏 3 4 b=3,4 , f ~ n , m , b , ε subscript ~ 𝑓 𝑛 𝑚 𝑏 𝜀 \widetilde{f}_{n,m,b,\varepsilon} and f ~ n , m , b , ε N subscript superscript ~ 𝑓 𝑁 𝑛 𝑚 𝑏 𝜀 \widetilde{f}^{N}_{n,m,b,\varepsilon} with b = 2 , 3 , 4 𝑏 2 3 4 b=2,3,4 and ε = 0.00001 𝜀 0.00001 \varepsilon=0.00001 .

			${\tilde{f}}_{n,m,b}$		${\tilde{f}}_{n,m,b,\epsilon }$, $\epsilon =0.00001$			${\tilde{f}}_{n,m,b,\epsilon }^N$, $\epsilon =0.00001$
	n	${\tilde{f}}_{n,m,m/2}$	$b=3$	$b=4$	$b=2$	$b=3$	$b=4$	$b=2$	$b=3$	$b=4$
$\left(a\right)$	50	0.08504	0.08687	0.08874	0.10597	0.10824	0.11057	0.09852	0.10064	0.10280
	200	0.02480	0.02533	0.02587	0.03090	0.03156	0.03224	0.02873	0.02935	0.02998
	500	0.01098	0.01122	0.01146	0.01368	0.01398	0.01428	0.01272	0.01299	0.01327
$\left(b\right)$	50	0.12469	0.12736	0.13010	0.13284	0.13569	0.13861	0.14319	0.14626	0.14940
	200	0.03636	0.03714	0.03794	0.03874	0.03957	0.04042	0.04176	0.04265	0.04357
	500	0.01610	0.01644	0.01680	0.01715	0.01752	0.01790	0.01849	0.01889	0.01929
$\left(c\right)$	50	0.05894	0.06020	0.06150	0.08278	0.08455	0.08637	0.08721	0.08908	0.09100
	200	0.01719	0.01755	0.01793	0.02414	0.02465	0.02518	0.02543	0.02598	0.02653
	500	0.00761	0.00777	0.00794	0.01069	0.01092	0.01115	0.01126	0.01150	0.01175
$\left(d\right)$	50	0.11515	0.11762	0.12015	0.14391	0.14699	0.15015	0.14656	0.14971	0.15292
	200	0.03358	0.03430	0.03504	0.04196	0.04286	0.04379	0.04274	0.04366	0.04459
	500	0.01487	0.01519	0.01551	0.01858	0.01898	0.01939	0.01892	0.01933	0.01975
$\left(e\right)$	50	0.15172	0.15498	0.15830	0.15311	0.15640	0.15976	0.15463	0.15795	0.16134
	200	0.04424	0.04519	0.04616	0.04465	0.04561	0.04659	0.04509	0.04606	0.04705
	500	0.01959	0.02001	0.02044	0.01977	0.02020	0.02063	0.01997	0.02040	0.02083
$\left(f\right)$	50	0.10527	0.10753	0.10984	0.10527	0.11383	0.11627	0.11420	0.11665	0.11916
	200	0.03070	0.03135	0.03203	0.03250	0.03319	0.03391	0.03330	0.03402	0.03475
	500	0.01359	0.01388	0.01418	0.01439	0.01470	0.01501	0.01475	0.01506	0.01539
$\left(g\right)$	50	0.06323	0.06458	0.06597	0.06298	0.06433	0.06571	0.06252	0.06386	0.06524
	200	0.01844	0.01883	0.01924	0.01836	0.01876	0.01916	0.01823	0.01862	0.01902
	500	0.00816	0.00834	0.00852	0.00813	0.00830	0.00848	0.00807	0.00824	0.00842
$\left(h\right)$	50	0.04133	0.04222	0.04312	0.03872	0.03955	0.04040	0.03566	0.03643	0.03721
	200	0.01205	0.01231	0.01257	0.01129	0.01153	0.01178	0.01040	0.01062	0.01085
	500	0.00533	0.00545	0.00557	0.00500	0.00510	0.00521	0.00460	0.00470	0.00480
$\left(i\right)$	50	0.03528	0.03603	0.03681	0.03594	0.03671	0.03750	0.03589	0.03666	0.03745
	200	0.01028	0.01051	0.01073	0.01048	0.01070	0.01093	0.01046	0.01069	0.01092
	500	0.00455	0.00465	0.00475	0.00464	0.00474	0.00484	0.00463	0.00473	0.00483
$\left(j\right)$	50	0.06328	0.06464	0.06603	0.06439	0.06577	0.06718	0.06236	0.06370	0.06507
	200	0.01845	0.01885	0.01925	0.01877	0.01918	0.01959	0.01818	0.01857	0.01897
	500	0.00817	0.00817	0.00852	0.00831	0.00849	0.00867	0.00805	0.00822	0.00840

Equations

$$\hat{f}(x)=\frac{1}{n{h}_n}\sum _{i=1}^{n}K\left(\frac{x-X_i}{h_n}\right)$$

$$\mathbb{E}\left[\hat{f}(x)\right]\simeq f(x)+\frac{1}{2}{h}^2{f}^{\prime \prime }(x)\int x^{2}K(x)dx,$$

$$Var\left[\hat{f}(x)\right]\simeq (nh{)}^{-1}f(x)\int K^2(x)dx.$$

$$\mathbb{E}\left[\hat{f}(x)\right]\simeq f(x){\int }_{-\infty }^{p}K(x)dx-f^{\prime }\left(x\right){\int }_{-\infty }^{p}xK(x)dx+\frac{1}{2}{h}^2{f}^{\prime \prime }(x){\int }_{-\infty }^p{x}^{2}K(x)dx,$$

$$Var\left[\hat{f}(x)\right]\simeq (nh{)}^{-1}f(x){\int }_{-\infty }^p{K}^2(x)dx.$$

$$\sum _{k=0}^{m}u\left(\frac{k}{m}\right)b_k(m,x)\to u(x),\text{mbox}uniformlyinx\in [0,1],$$

$${\tilde{F}}_n(x)=\sum _{k=0}^m{F}_n\left(\frac{k}{m}\right)b_k(m,x),$$

$${\tilde{f}}_n(x)$$

$$T_{n,m}(x)$$

$$f_n(x)=(1-{\gamma }_n)f_{n-1}(x)+{\gamma }_n{Z}_n\left(x\right),$$

$$f_n(x)={\Pi }_n\sum _{k=1}^n{\Pi }_k^{-1}{\gamma }_k{Z}_k(x).$$

$$\underset{n\to +\infty }{\lim }n\left[1-\frac{v_{n-1}}{v_n}\right]=\gamma .$$

$${\Delta }_1(x)$$

$${\Delta }_2(x)$$

$$C_1$$

$$C_4$$

$$C_6$$

$$Bias\left[f_n(x)\right]=-m_n^{-2}\frac{2}{1-2a\xi }{\Delta }_2(x)+o\left(m_n^{-2}\right).$$

$$Bias\left[f_n(x)\right]=o\left(\sqrt{{\gamma }_n{m}_n^{1/2}}\right).$$

$$Var[f_n(x)]=C_3{\gamma }_n{m}_n^{1/2}\frac{2}{4-(2\alpha -a)\xi }f(x)\psi (x)+o\left({\gamma }_n{m}_n^{1/2}\right).$$

$$Var[f_n(x)]=o\left(m_n^{-4}\right).$$

$$MSE\left[f_n(x)\right]={\Delta }_2^2(x)m_n^{-4}\frac{4}{(1-2a\xi {)}^2}+o\left(m_n^{-4}\right).$$

$$MSE\left[f_n(x)\right]$$

$$m_n=\underset{m_n}{\mathrm{argmin}}MSE\left[f_n\left(x\right)\right],\text{mbox}forx\in (0,1).$$

$$\left(2^{2/9}{\left(1-\frac{4}{9}\xi \right)}^{-2/9}{\left[\frac{32{\Delta }_2^2(x)}{C_{3}f(x)\psi (x)}\right]}^{2/9}{\gamma }_n^{-2/9}\right),$$

$$MSE\left[f_n(x)\right]=\frac{9(32C_3^8{)}^{1/9}}{8}\frac{({\Delta }_2(x){)}^{2/9}(f(x)\psi (x){)}^{8/9}}{2^{2/9}{\left(1-\frac{4}{9}\xi \right)}^{10/9}}{\gamma }_n^{\frac{8}{9}}+o\left({\gamma }_n^{\frac{8}{9}}\right).$$

$$\left(m_n\right)=\left(2^{2/9}{\left({\gamma }_0-4/9\right)}^{-2/9}{\left[\frac{32{\Delta }_2^2(x)}{C_{3}f(x)\psi (x)}\right]}^{2/9}{n}^{2/9}\right),$$

$$MSE\left[f_n(x)\right]=\frac{9(32C_3^8{)}^{1/9}({\Delta }_2(x){)}^{2/9}(f(x)\psi (x){)}^{8/9}}{8}\frac{{\gamma }_0^2}{2^{8/9}{\left({\gamma }_0-4/9\right)}^{10/9}}{n}^{-8/9}+o\left(n^{-8/9}\right).$$

$$Bias\left[f_n(x)\right]=-m_n^{-2}\frac{2}{1-2a\xi }{\Delta }_2(x)+o\left(m_n^{-2}\right).$$

$$Bias\left[f_n(x)\right]=o\left(\sqrt{{\gamma }_n{m}_n}\right).$$

$$Var[f_n(x)]=\frac{5}{2}{\gamma }_n{m}_n\frac{1}{2-(\alpha -a)\xi }f(x)+o\left({\gamma }_n{m}_n\right).$$

Full text

In this paper, we consider the alleviation of the boundary problem when the probability density function has bounded support. We apply Robbins-Monro’s algorithm and Bernstein polynomials to construct a recursive density estimator. We study the asymptotic properties of the proposed recursive estimator. We then compared our proposed recursive estimator with many others estimators. Finally, we confirm our theoretical result through a simulation study and then using two real datasets.*

**Keywords: Density estimation; Stochastic approximation algorithm; Bernstein polynomial.

Mathematics Subject Classification: 62G07, 62L20.**

1 Introduction

There has been a considerable development of methods for smooth estimation of density and distribution functions, following the introduction of several kernel smoothing by Rosenblatt (1956) and the further advances made on kernel method by Parzen (1962). We advise the reader to see the paper of Härdle (1991) for an introduction of several kernel smoothing techniques. However, these methods have difficulties at and near boundaries when curve estimation is attempted over a region with boundaries. Moreover, it is well known in nonparametric kernel density estimation that the bias of the standard kernel density estimator

[TABLE]

is of a larger order near the boundary than that in the interior, where, $K$ is a kernel (that is, a function satisfying $\int_{\mathbb{R}}K\left(x\right)dx=1$), and $\left(h_{n}\right)$ is a bandwidth (that is, a sequence of positive real numbers that goes to zero). We suppose for simplicity that there is a single known boundary to the support of the density function $f$ which we might as well take to be at the origin, then we are dealing with positive data. For convenience, we consider a symmetric kernel (for instance, normal kernels). Away from the boundary, which means that at any $x>h_{n}$, the usual asymptotic mean and variance expressions are applied. Let us now suppose that $f$ has two continuous derivatives everywhere, and that as $n\to\infty$, $h=h_{n}\to 0$ and $nh\to 0$. Then,

[TABLE]

and

[TABLE]

Near the boundary, the expression of the mean and the variance are different. Let $x=ph$, we have

[TABLE]

and

[TABLE]

These bias phenomena are called boundary bias. Many authors have suggested methods for reducing this phenomena such as data reflection (Schuster (1985)), boundary kernels (Müller (1991, 1993) and Müller and Wand (1994)), the local linear estimator (Lejeune and Sarda (1992) and Jones (1993)), the use of beta and gamma kernels (Chen (1999, 2000)).

For a smooth estimate of a density function with finite known support, Vitale’s method (Vitale (1975)) based on the Bernstein polynomials, illustrated below, also has been investigated in the literature (Ghosal (2001), Babu et al. (2002), Kakizawa (2004), Rao (2005)) and more recently by Leblanc (2010) and Igarashi and Kakizawa (2014)). The idea comes from the Weierstrass’s approximation theorem that for any continuous function $u$ on the interval $[0,1]$, we have

[TABLE]

where $b_{k}(m,x)={{m}\choose{k}}x^{k}(1-x)^{m-k}$ is the Bernstein polynomial of order $m$ .

In the context of distribution function $F$ with support $[0,1]$, Vitale (1975) proposed an estimator

[TABLE]

where $F_{n}$ is the empirical distribution based on a random sample $X_{1},X_{2},\ldots,X_{n}$. Hence, an estimator for the density $f$ is given by

[TABLE]

In this paper, we propose a recursive method to estimate the unknown density function $f$. The advantage of recursive estimators is that their update from a sample of size $n$ to one of size $n+1$, requires considerably less computations. This property is particularly important, since the number of points at which the function is estimated is usually very large.

Let $X_{1},X_{2},\ldots,X_{n}$ be a sequence of i.i.d random variables having a common unknown distribution $F$ with associated density $f$ supported on $[0,1]$. In order to construct a recursive method to estimate the unknown density $f$, we first follow Leblanc (2010) and we introduce $T_{n,m}$ as follows:

[TABLE]

Here, we let $k_{n}=[mX_{n}]$, where $[x]$ denotes the largest integer smaller than $x$. Then, we use Robbins-Monro’s scheme (Robbins and Monro (1951)), and we set $f_{0}(x)\in\mathbb{R}$ and for all $n\in\mathbb{N}^{*}$, we set

[TABLE]

where $(\gamma_{n})$ is a sequence of real numbers, called a stepsize and $Z_{n}(x)=2T_{n,m}(x)-T_{n,m/2}(x)$. Then for simplicity, we suppose that $f_{0}(x)=0$ and $\Pi_{n}=\prod_{j=1}^{n}(1-\gamma_{j})$.

Then, it follows from (2), that one can estimate $f$ recursively at the point $x$ by

[TABLE]

Our first aim in this paper is to compute the bias, the variance, the mean squared error ($MSE$) and the mean integrated squared error ($MISE$) of our proposed recursive estimators. It turns out that they heavily depend on the choice of the stepsize $(\gamma_{n})$. Moreover, we give the optimal order $(m_{n})$ which minimizes the $MSE$ and the $MISE$ of the proposed recursive estimators. Further, we show that using the stepsize $(\gamma_{n})=(n^{-1})$ and the optimal order $\left(m_{n}\right)$, the proposed estimator $f_{n}$ can dominate Vitale’s estimator $\widetilde{f}_{n}$ in terms of $MISE$. Finally, we confirm our theoretical results from a simulation study.

The remainder of this paper is organized as follows. In Section 2, we list our notations and assumptions. In Section 3, we state the main theoretical results regarding bias, variance, $MSE$ and $MISE$. Section 4 is devoted to some numerical studies : first by simulation (Subsection 4.1) and second using some real datasets (Subsection 4.2). We conclude the paper in Section 5. Appendix A gives the proof of our theoretical results.

2 Assumptions and Notations

Definition 1*.*

Let $\gamma\in\mathbb{R}$ and $(v_{n})_{n\geq 1}$ be a nonrandom positive sequence. We say that $(v_{n})\in\mathcal{GS}(\gamma)$ if

[TABLE]

This condition was introduced by Galambos and Seneta (1973) to define regularly varying sequences (see also Bojanic and Seneta (1973)). Typical sequences in $\mathcal{GS}(\gamma)$ are, for $b\in\mathbb{R}$, $n^{\gamma}(\log n)^{b}$, $n^{\gamma}(\log\log n)^{b}$, and so on.

To study the estimator (3), we make the following assumptions :

$\left(A1\right)$ $f$ admits a continuous fourth-order derivative $f^{(4)}$ on $[0,1]$.

$\left(A2\right)$ $(\gamma_{n})\in\mathcal{GS}\left(-\alpha\right)$, $\alpha\in(\frac{1}{2},1]$.

$\left(A3\right)$ $(m_{n})\in\mathcal{GS}(a)$, $a\in(0,1)$.

$\left(A4\right)$ $\lim_{n\rightarrow\infty}(n\gamma_{n})\in(\min\left(2a,(2\alpha-a)/4\right),\infty)$.

• •

Assumption $(A1)$ is standard in the framework of nonparametric estimation of probability density using Bernstein polynomials (see Leblanc (2010)).

• •

Assumption $(A2)$ on the stepsize is usual in the framework of the recursive estimation for density estimation (see Mokkadem et al. (2009) and Slaoui (2013, 2014a, 2014b)). This assumption ensures that $\sum_{n\geq 1}\gamma_{n}=\infty$ and $\sum_{n\geq 1}\gamma_{n}^{2}<\infty$, which are two classical assumptions for obtaining the convergence of Robbins-Monro’s algorithm (see Duflo (1997)).

• •

Assumption $(A3)$ on $(m_{n})$ was introduced similarly to the assumption on the bandwidth used for the recursive kernel distribution estimator (see Slaoui (2014a, b)) to ensure the application of the technical lemma given in the appendix A.

• •

Assumption $(A4)$ on the limit of $\left(n\gamma_{n}\right)$ as $n$ goes to infinity is usual in the framework of stochastic approximation algorithms. This condition ensures the application of the technical lemma given in the appendix A to obtain the asymptotic bias and variance respectively.

Remark 1*.*

The intuition behind the use of such order $\left(m_{n}\right)$ belonging to $\mathcal{GS}\left(a\right)$ is that the ratio $m_{n-1}/m_{n}$ is equal to $1-a/n+o\left(1/n\right)$, then using such order and using the assumption on the stepsize, which ensures that $\gamma_{n-1}/\gamma_{n}$ is equal to $1+\alpha/n+o\left(1/n\right)$, the application of the technical lemma given in the appendix A ensures that the bias and the variance will depend only on $m_{n}$ and $\gamma_{n}$ and not on $m_{1},\ldots,m_{n}$ and $\gamma_{1},\ldots,\gamma_{n}$, then the $MISE$ will depend also only on $m_{n}$ and $\gamma_{n}$, which will be helpful to deduce an optimal order and an optimal stepsize.

Throughout this paper we will use the following notations:

[TABLE]

3 Main results

In this section we start by developing the bias, the variance and the $MSE$ of the proposed recursive estimators using Bernstein estimator first in the boundary region and then at the edges.

3.1 Within the interval $[0,1]$

The following proposition gives the bias, variance, and $MSE$ of the proposed recursive estimator $f_{n}(x)$ for $x\in(0,1)$.

Proposition 1*.*

Let Assumptions $\left(A1\right)-\left(A4\right)$ hold. For $x\in(0,1)$, we have

If $a\in\left(0,\frac{2}{9}\alpha\right]$, then

[TABLE]

If $a\in\left(\frac{2}{9}\alpha,1\right)$, then

[TABLE] 2. 2.

If $a\in\left[\frac{2}{9}\alpha,1\right)$, then

[TABLE]

If $a\in\left(0,\frac{2}{9}\alpha\right]$, then

[TABLE] 3. 3.

If $\lim_{n\rightarrow\infty}(n\gamma_{n})>\max\left(2a,(2\alpha-a)/4\right)$, that implies $\lim_{n\rightarrow\infty}(n\gamma_{n})>2a$ and $\lim_{n\rightarrow\infty}(n\gamma_{n})>(2\alpha-a)/4$, in this case, we have $1-2a\xi>0$ and $4-(2\alpha-a)\xi>0$ then (4) and (6) hold simultaneously. 4. 4.

If $a\in\left(0,\frac{2}{9}\alpha\right)$, then

[TABLE]

If $a=\frac{2}{9}\alpha$, then

[TABLE]

If $a\in\left(\frac{2}{9}\alpha,1\right]$, then

[TABLE]

In order to choose the optimal order $\left(m_{n}\right)$, we choose $\left(m_{n}\right)$ such that

[TABLE]

Then, it follows from (8), that $\left(m_{n}\right)$ must be equal to

[TABLE]

and then the corresponding $MSE$ is equal to

[TABLE]

Moreover, since the optimal stepsize should be obtained by minimizing the $MSE$ of our proposed estimators $f_{n}\left(x\right)$, then $\left(\gamma_{n}\right)$ must be chosen in $\mathcal{GS}\left(-1\right)$. By considering the case when $(\gamma_{n})=\left(\gamma_{0}n^{-1}\right)$, we obtain

[TABLE]

and the corresponding $MSE$ is equal to

[TABLE]

3.2 The edges of the interval $[0,1]$

For the cases $x=0,1$, we need the following additional assumption :

$\left(A^{\prime}4\right)$ $\quad\lim_{n\rightarrow\infty}(n\gamma_{n})\in(\min\left(2a,(\alpha-a)/2\right),\infty)$.

The following proposition gives bias, variance and $MSE$ of $f_{n}(x)$, for $x=0,1$.

Proposition 2*.*

Let Assumptions $\left(A1\right)-\left(A3\right)$ and $\left(A^{\prime}4\right)$ hold, for $x=0,1$, we have

If $a\in\left(0,\frac{\alpha}{5}\right]$, then

[TABLE]

If $a\in\left(\frac{\alpha}{5},1\right)$, then

[TABLE] 2. 2.

If $a\in\left[\frac{\alpha}{5},1\right)$, then

[TABLE]

If $a\in\left(0,\frac{\alpha}{5}\right)$, then

[TABLE] 3. 3.

If $\lim_{n\rightarrow\infty}(n\gamma_{n})>\max\left(2a,(\alpha-a)/2\right)$, that implies $\lim_{n\rightarrow\infty}(n\gamma_{n})>2a$ and $\lim_{n\rightarrow\infty}(n\gamma_{n})>(\alpha-a)/2$ then $1-2a\xi>0$ and $2-(\alpha-a)\xi>0$, then, (9) and (11) hold simultaneously. 4. 4.

If $a\in\left(0,\frac{\alpha}{5}\right)$, then

[TABLE]

If $a=\frac{\alpha}{5}$, then

[TABLE]

If $a\in\left(\frac{\alpha}{5},1\right)$, then

[TABLE]

In order to choose the optimal order $\left(m_{n}\right)$, we choose $\left(m_{n}\right)$ such that

[TABLE]

Then, it follows from (13), that $\left(m_{n}\right)$ must be equal to

[TABLE]

and then the corresponding $MSE$ is equal to

[TABLE]

Moreover, since the optimal stepsize should be obtained by minimizing the $MSE$ of our proposed estimators $f_{n}\left(x\right)$, then $\left(\gamma_{n}\right)$ must be chosen in $\mathcal{GS}\left(-1\right)$. By considering the case when $(\gamma_{n})=\left(\gamma_{0}n^{-1}\right)$, we obtain

[TABLE]

and then the corresponding $MSE$ is equal to

[TABLE]

In the next subsection we give the $MISE$ of the proposed recursive estimator $f_{n}$ introduced in (2) for $x\in\left(0,1\right)$.

3.3 $MISE$ of the recursive estimator $f_{n}$

The following proposition gives the $MISE$ of the proposed recursive estimator $f_{n}$.

Proposition 3*.*

Let Assumptions $\left(A1\right)-\left(A4\right)$ hold. We have

If $a\in\left(0,\frac{2}{9}\alpha\right)$, then

[TABLE] 2. 2.

If $a=\frac{2}{9}\alpha$, then

[TABLE] 3. 3.

If $a\in\left(\frac{2}{9}\alpha,1\right)$, then

[TABLE]

The following result is a consequence of the previous proposition which gives the optimal order $\left(m_{n}\right)$ of the estimator (3) and the corresponding $MISE$.

Corollary 1*.*

Let Assumptions $\left(A1\right)-\left(A4\right)$ hold. To minimize the $MISE$ of $f_{n}$, the stepsize $(\gamma_{n})$ must be chosen in $\mathcal{GS}(-1)$ and $(m_{n})$ must be in $\mathcal{GS}(2/9)$ such that

[TABLE]

and then the corresponding $MISE$ is equal to

[TABLE]

Moreover, in the case when $(\gamma_{n})=\left(\gamma_{0}n^{-1}\right)$, the optimal order $\left(m_{n}\right)$ must be equal to

[TABLE]

and the corresponding $MISE$ is given by

[TABLE]

Remark 2*.*

The minimum of $\frac{\gamma_{0}^{2}}{\left(\gamma_{0}-4/9\right)^{10/9}}$ is reached for $\gamma_{0}=1$. Moreover, to minimize the variance of $f_{n}$, we should choose $\gamma_{0}=1-\frac{a}{2}$, with $a=2/9$ for $x\in(0,1)$ and $\gamma_{0}=1-a$ with $a=1/5$ for $x=0,1$. Therefore, in the application section, we will consider the following stepsizes $(\gamma_{n})=(n^{-1})$, $(\gamma_{n})=\left(\frac{8}{9}n^{-1}\right)$ and $(\gamma_{n})=\left(\frac{4}{5}n^{-1}\right)$.

Let us now state the following theorem, which gives the weak convergence rate of the estimator $f_{n}$ defined in (3).

Theorem 1* (Weak pointwise convergence rate).*

Let Assumption $\left(A1\right)-\left(A4\right)$ hold. For $x\in(0,1)$, we have

If $\gamma^{-1/2}_{n}m^{-9/4}_{n}\rightarrow c$ for some constant $c\geq 0$, then

[TABLE] 2. 2.

If $\gamma^{-1/2}_{n}m^{-9/4}_{n}\rightarrow\infty$, then

[TABLE]

where $\stackrel{{\scriptstyle\mathcal{D}}}{{\rightarrow}}$ denotes the convergence in distribution, $\mathcal{N}$ the Gaussian-distribution and $\stackrel{{\scriptstyle\mathbb{P}}}{{\rightarrow}}$ the convergence in probability.

Remark 3*.*

When the bandwidth $\left(h_{n}\right)$ is chosen such that $\lim_{n\to\infty}\gamma^{-1/2}_{n}m^{-9/4}_{n}=0$ and the stepsize such that $\lim_{n\to\infty}n\gamma_{n}=\gamma_{0}$, the proposed recursive estimators fulfills the following central limit theorem

[TABLE]

3.4 Results on some classical Bernstein density estimator

In the next paragraph, we recall some results on Vitale’s Bernstein density estimator $\widetilde{f}_{n}$ given in (1).

3.4.1 Vitale’s Bernstein density estimator $\widetilde{f}_{n}$

Under some classical assumption on the density, such as $f$ is continuous and admits two continuous and bounded derivatives, for $x\in[0,1]$, we have

[TABLE]

[TABLE]

Moreover, we have

[TABLE]

To minimize the $MISE$ of $\widetilde{f}_{n}$, the order $(m_{n})$ must equal to

[TABLE]

and the corresponding $MISE$

[TABLE]

Remark 4*.*

If we let $h=m^{-1}$ be the bandwidth of the Vitale’s estimator $\widetilde{f}_{n}$, it is clear that the bias of $\widetilde{f}_{n}$ is $O\left(h^{-1}\right)$ which is larger than the classical kernel density estimators which is $O(h^{2})$ except near the boundaries. To reduce this bias, Leblanc (2010) suggest a new Bernstein estimator using the method of bias correction. This methodology was used first by Politis and Romano (1995) in the context of spectral density estimation and is linked with the work of Schucany, Gray and Owen (1971) and Schucany and Sommers (1977) on bias reduction in estimation.

3.4.2 Bias correction for Bernstein density $\widetilde{f}_{n,m,m/2}$

The bias correction for Bernstein density proposed by Leblanc Leblanc (2010) is defined by

[TABLE]

where $\widetilde{f}_{n,m}$ and $\widetilde{f}_{n,m/2}$ are the Bernstein density estimators introduced by Vitale with order $m$ and $m/2$ respectively which defined in (1). Let us now recall the characteristics of the estimator $\widetilde{f}_{n,m,m/2}$. Under the Assumption $(A1)$, we have

[TABLE]

[TABLE]

and then

[TABLE]

To minimize the $MISE$ of $\widetilde{f}_{n,m,m/2}$, the order $(m_{n})$ must equal to

[TABLE]

and then the corresponding $MISE$

[TABLE]

Igarashi and Kakizawa (2014) have generalized the estimator proposed by Leblanc $\widetilde{f}_{n,m,m/2}$ and defined in (17)

[TABLE]

where $b=2,3,\ldots$ and $\widetilde{f}_{n,m}$ and $\widetilde{f}_{n,m/b}$ are the Vitale’s density estimators defined in (1) with order $m$ and $m/b$ respectively. Under the Assumption $\left(A1\right)$, we have

[TABLE]

[TABLE]

[TABLE]

Moreover, we have

[TABLE]

To minimize the $MISE$ of $\widetilde{f}_{n,m,m/b}$, the order $(m_{n})$ must equal to

[TABLE]

and then the corresponding $MISE$

[TABLE]

Remark 5*.*

The equation (21) indicates that the choice $b=2$ is the best choice in terms of the $MISE$ for the estimator $\widetilde{f}_{n,m,m/b}$, since the factor $\left(b\lambda^{4}_{1}(b)\right)^{2/9}$ is increasing in $b=2,3,\ldots$.

This method of bias correction reduces the bias of Bernstein estimator from $O\left(m^{-1}\right)$ to $O\left(m^{-2}\right)$, but it loses the non-negativity. As an additive bias correction to the logarithm of estimator, Terrell and Scott (1980) originally developed a multiplicative bias correction that enables keeping the non-negativity. This method was adopted by Hirukawa (2010) for the beta kernel estimator introduced by Chen (1999).

3.4.3 Multiplicative bias-correction Bernstein density estimator $\widetilde{f}_{n,m,b,\varepsilon}$

This estimator was considered by Igarashi and Kakizawa (2014) by applying a multiplicative bias correction method to the Bernstein estimator

[TABLE]

for some $\varepsilon=\varepsilon(m)>0$, converting to zero at a suitable rate. Under the Assumption $(A1)$, for $x\in[0,1]$ such as $f(x)>0$, with $m=O\left(n^{\eta}\right)$ and $\varepsilon\approx m^{\tau}$ where $\eta\in(0,1)$ and $\tau>2$ we have

[TABLE]

and

[TABLE]

Moreover, we have

[TABLE]

To minimize the $MISE$ of $\widetilde{f}_{n,m,b,\varepsilon}$, the order $\left(m_{n}\right)$ must equal to

[TABLE]

and the corresponding $MISE$ is equal to

[TABLE]

Remark 6*.*

Note that the estimator $\widetilde{f}_{n,m,b,\varepsilon}$ retains non-negativity, but it is not a genuine density. In fact, $\widetilde{f}_{n,m,b,\varepsilon}$ does not generally integrate to unity. To solve this problem, Igarashi and Kakizawa (2014) proposed the normalized bias-corrected Bernstein estimator.

3.4.4 Normalized bias-corrected Bernstein estimator $\widetilde{f}^{N}_{n,m,b,\varepsilon}$

The normalized bias-corrected Bernstein estimator is given by:

[TABLE]

Under the Assumption $(A1)$, for $x\in[0,1]$ such as $f(x)>0$, with $m=O\left(n^{\eta}\right)$ and $\varepsilon\approx m^{\tau}$ where $\eta\in(0,1)$ and $\tau>2$ we have

[TABLE]

and

[TABLE]

and then, we have

[TABLE]

To minimize the $MISE$ of $\widetilde{f}^{N}_{n,m,b,\varepsilon}$, the order $\left(m_{n}\right)$ must equal to

[TABLE]

and the corresponding $MISE$ is equal to

[TABLE]

4 Applications

When using the Bernstein polynomials, we must consider a density on $[0,1]$. For this purpose, we need to make some suitable transformations in different cases (we list below) :

Suppose that $X$ is concentrated on a finite support $[a,b]$. Then we work with the sample values $Y_{1},\ldots,Y_{n}$ where $Y_{i}=\frac{X_{i}-a}{b-a}$. Denoting by $g_{n}(x)$ the estimated density function of $Y_{1},\ldots,Y_{n}$, we compute the estimated density $f_{n}$ of $X$

[TABLE] 2. 2.

For the densities functions concentrated on the interval $(-\infty,+\infty)$, we can use the transformed sample $Y_{i}=\frac{1}{2}+\frac{1}{\pi}\arctan(X_{i})$ which transforms the range to the interval $(0,1)$. Hence

[TABLE] 3. 3.

For the support $[0,\infty)$, we can use the transformed sample $Y_{i}=\frac{X_{i}}{X_{i}+1}$ which transforms the range to the interval $(0,1)$. Hence

[TABLE]

Computational cost

As mentioned in the introduction, the advantage of the proposed recursive estimators on their non-recursive version is that their update, from a sample of size $n$ to one of size $n+1$, require less computations. This property can be generalized, one can check that it follows from (2) that for all $n_{1}\in\left[0,n-1\right]$,

[TABLE]

where $\alpha_{1}=\prod_{j=n_{1}+1}^{n}\left(1-\gamma_{j}\right)$ and $\beta_{k}=\gamma_{k}\prod_{j=k+1}^{n}\left(1-\gamma_{j}\right)$. It is clear, that the proposed estimators can be viewed as a linear combination of two estimators, which improve considerably the computational cost. Moreover, in order to give some comparative elements with Vitale’s estimator defined in (1), including computational coasts, we consider $500$ samples of size $500$ generated from the beta distribution $\mathcal{B}\left(3,5\right)$; moreover, we suppose that we receive an additional $500$ samples of size $500$ generated also from the beta distribution $\mathcal{B}\left(3,5\right)$. Performing the two methods, the running time using our proposed recursive estimator defined by algorithm (2) with stepsize $\left(\gamma_{n}\right)=\left(n^{-1}\right)$ and the order $\left(m_{n}\right)$ according to the minimization of Least Squares Cross-Validation (LSCV) described below was roughly $25$ minutes on the author’s workstation, while the running time using the non-recursive estimator (1) was roughly more than $4$ hours on the author’s workstation.

The aim of this paragraph is to compare the performance of Vitale’s estimator $\widetilde{f}_{n}$ defined in (1), Leblanc’s estimator $\widetilde{f}_{n,m,m/2}$ given in (17), the generalized estimator $\widetilde{f}_{n,m,m/b}$, defined in (19), the multiplicative bias corrected Bernstein estimator $\widetilde{f}_{n,m,b,\varepsilon}$ given in (22) and the normalized estimator $\widetilde{f}^{N}_{n,m,b,\varepsilon}$ given in (24) with that of the proposed estimator (3).

• (1)

When applying $f_{n}$, one needs to choose two quantities:

• $\bullet$

The stepsize $(\gamma_{n})$ is chosen to be equal to $\left(\nu n^{-1}\right)$, with $\nu\in\left\{\frac{4}{5},\frac{8}{9},1\right\}$.

• $\bullet$

The order $(m_{n})$ is chosen to be equal to (14).

• (2)

When applying $\widetilde{f}_{n}$, one needs to choose the order $\left(m_{n}\right)$ to be equal to (16).

• (3)

When applying $\widetilde{f}_{n,m,m/2}$, one needs to choose the order $\left(m_{n}\right)$ to be equal to (18).

• (3)

When applying $\widetilde{f}_{n,m,m/b}$, one needs to choose the order $\left(m_{n}\right)$ to be equal to (20) and $b=3,4$.

• (4)

When applying $\widetilde{f}_{n,m,b,\varepsilon}$, one needs to choose the order $\left(m_{n}\right)$ to be equal to (23), $b=2,3,4$ and $\varepsilon=0.00001$.

• (5)

When applying $\widetilde{f}^{N}_{n,m,b,\varepsilon}$, one needs to choose the order $\left(m_{n}\right)$ to be equal to (25), $b=2,3,4$ and $\varepsilon=0.00001$.

4.1 Simulations

We consider the following ten density functions :

(a) the beta density $\mathcal{B}(3,5)$, $f(x)=\frac{x^{2}(1-x)^{4}}{B(3,5)}$,

(b) the beta density $\mathcal{B}(1,6)$, $f(x)=\frac{(1-x)^{5}}{B(1,6)}$,

(c) the beta density $\mathcal{B}(3,1)$, $f(x)=\frac{x^{2}}{B(3,1)}$,

(d) the beta mixture density $1/2\mathcal{B}(3,9)+1/2\mathcal{B}(9,3)$, $f(x)=0.5\frac{x^{2}(1-x)^{8}}{B(3,9)}+0.5\frac{x^{8}(1-x)^{2}}{B(9,3)}$,

(e) the beta mixture density $1/2\mathcal{B}(3,1)+1/2\mathcal{B}(10,10)$, $f(x)=0.5\frac{x^{2}}{B(3,1)}+0.5\frac{x^{9}(1-x)^{9}}{B(10,10)}$,

(f) the beta mixture density $1/2\mathcal{B}(1,6)+1/2\mathcal{B}(3,5)$, $f(x)=0.5\frac{(1-x)^{5}}{B(1,6)}+0.5\frac{x^{2}(1-x)^{4}}{B(3,5)}$,

(g) the beta mixture density $1/2\mathcal{B}(2,1)+1/2\mathcal{B}(1,4)$, $f(x)=0.5\frac{x}{B(2,1)}+0.5\frac{(1-x)^{3}}{B(1,4)}$,

(h) the truncated exponential density $\mathcal{E}_{[0,1]}(1/0.8)$, $f(x)=\frac{\exp(-x/0.8)}{0.8\left\{1-\exp(-1/0.8)\right\}}$,

(i) the truncated normal density $\mathcal{N}_{[0,1]}(0,1)$, $f(x)=\frac{\exp(-x^{2}/2)}{\displaystyle\int_{0}^{1}\exp(-t^{2}/2)dt}$,

(j) the truncated normal mixture density $1/4\mathcal{N}_{[0,1]}(2,1)+3/4\mathcal{N}_{[0,1]}(-3,1)$,

$f(x)=0.25\frac{\exp(-(x-2)^{2}/2)}{\displaystyle\int_{0}^{1}\exp(-(t-2)^{2}/2)dt}+0.75\frac{\exp(-(x+3)^{2}/2)}{\displaystyle\int_{0}^{1}\exp(-(t+3)^{2}/2)dt}$.

For each density function and sample of size $n$, we approximate the average integrated squared error ($ISE$) of the estimator using $N=500$ trials of sample size $n$; $\overline{ISE}=\frac{1}{N}\displaystyle\sum_{k=1}^{N}ISE\left[\hat{f}_{k}\right]$, where $\hat{f}_{k}$ is the estimator computed from the $k$th sample, and, $ISE[\hat{f}]=\displaystyle\int_{0}^{1}\left\{\hat{f}(x)-f(x)\right\}^{2}dx$.

From Table 1 and 2 we conclude that :

• •

In all the considered cases, the average $ISE$ of our density estimator (3) is smaller than that of Vitale’s estimator defined in (1), except the cases $(e)$ and $(f)$ of the Beta mixture and the cases $(a)$, $(d)$, $(g)$ and $(h)$ for the small sample size $n=50$ and in the case $(i)$ for the size $n=50$ and $n=200$.

• •

In all the cases, the average $ISE$ of our recursive density estimator (3) is slightly larger then that of Leblanc’s estimator $\widetilde{f}_{n,m,m/2}$ given in (17).

• •

In all the cases, the average $ISE$ of our density estimator (3) is smaller than that of the generalized estimator $\widetilde{f}_{n,m,m/b}$, with $b=4$ (see (19)).

• •

In all the cases, the average $ISE$ of our density estimator (3) is smaller than that of the multiplicative bias corrected Bernstein estimator $\widetilde{f}_{n,m,b,\varepsilon}$ given in (22) and the normalized estimator $\widetilde{f}^{N}_{n,m,b,\varepsilon}$ given in (24), except the cases $(g)$ and $(h)$.

• •

The average $ISE$ of the generalized estimator $\widetilde{f}_{n,m,m/b}$, (19) increase when $b$ increase, then the optimal choice is $b=2$ which corresponds to Leblanc’s estimator $\widetilde{f}_{n,m,m/2}$ given in (17).

• •

The average $ISE$ of the multiplicative bias corrected Bernstein estimator $\widetilde{f}_{n,m,b,\varepsilon}$ given in (22) and the normalized estimator $\widetilde{f}^{N}_{n,m,b,\varepsilon}$ defined in (24) increase when $b$ increase

• •

The average $ISE$ of the multiplicative bias corrected Bernstein estimator $\widetilde{f}_{n,m,b,\varepsilon}$ given in (22) and the normalized estimator $\widetilde{f}^{N}_{n,m,b,\varepsilon}$ given in (24) are smaller than that of the estimator $\widetilde{f}_{n,m,m/b}$ defined in (19), in the cases $(g)$ and $(h)$ and are larger in the other cases.

• •

The average $ISE$ of the multiplicative bias corrected Bernstein estimator $\widetilde{f}_{n,m,b,\varepsilon}$ (22) is larger than that of the normalized estimator $\widetilde{f}^{N}_{n,m,b,\varepsilon}$ given in (24) in the cases $(a)$, $(g)$, $(h)$, $(i)$ and $(j)$ and is smaller in the other cases.

• •

The average $ISE$ decreases as the sample size increases.

From figure 1, we can observe that :

• •

The proposed recursive estimator $f_{n}$ given in (3) using the stepsize $\left(\gamma_{n}\right)=\left(n^{-1}\right)$ is closer to the true density function compared to Vitale’s estimator $\widetilde{f}_{n}$ given in (1), the generalized estimator $\widetilde{f}_{n,m,m/4}$ given in (19), Kakizwa’s estimator $\widetilde{f}_{n,m,b,\varepsilon}$ defined in (22) and the normalized estimator $\widetilde{f}^{N}_{n,m,b,\varepsilon}$ defined in (24) with $b=2$ and $\varepsilon=0.00001$.

• •

Within the interval $[0,1]$, our density estimator (3) using the stepsize $\left(\gamma_{n}\right)=\left(n^{-1}\right)$ is closer to the true density function compared to Leblanc’s estimator $\widetilde{f}_{n,m,m/2}$ defined in (17).

• •

When the sample size increases, we get closer estimation of the true density function.

4.2 Real dataset

In any practical situation, to estimate an unknown density function $f$, it is essential to specify the order $m$ to be used for the estimator. One way is to use least squares cross-validation ($LSCV$) method to obtain a data-driven choice of $m$.

4.2.1 Order selection method

First, we recall that the $LSCV$ method is based on minimizing the integrated square error between the estimated density function $\hat{f}$ and the true density function $f$

[TABLE]

From this, Silverman (1986) derived the score function

[TABLE]

where $\hat{f}_{-i}$ is the density estimate without the data point $X_{i}$. The smoothing parameter is chosen by minimizing $LSCV(m)$ ($m=\arg\min_{m}LSCV(m)$).

For our proposed recursive estimator $f_{n}$ given in (3), we make the following choice of $(\gamma_{n})=(n^{-1})$, and then we choice the order $m$ in order to minimize the following criterion :

[TABLE]

We define integer sequences $p_{i}=[m_{i}X_{i}]$ and $q_{i}=[m_{i}X_{i}/2]$, we then have $X_{i}\in\left(\frac{p_{i}}{m_{i}},\frac{p_{i}+1}{m_{i}}\right]$ and $X_{i}\in\left(\frac{2q_{i}}{m_{i}},\frac{2(q_{i}+1)}{m_{i}}\right]$. Then we obtain

[TABLE]

then we conclude that

[TABLE]

Note that, the $LSCV$ function for Vitale’s estimator (1), is written as

[TABLE]

We define integer sequences $k_{i}=[mX_{i}]$, we then have $X_{i}\in\left(\frac{k_{i}}{m},\frac{k_{i}+1}{m}\right]$, consequently

[TABLE]

then we conclude

[TABLE]

The $LSCV$ function for the estimator $\widetilde{f}_{n,m,m/b}$ defined in (19) with $b=2,3,4,...$, is written as

[TABLE]

We define integer sequences $k_{i}=[mX_{i}]$ and $r_{i}=[mX_{i}/b]$, we then have $X_{i}\in\left(\frac{k_{i}}{m},\frac{k_{i}+1}{m}\right]$ and $X_{i}\in\left(\frac{br_{i}}{m},\frac{b(r_{i}+1)}{m}\right]$. Then we obtain

[TABLE]

Using the Kakizawa’s estimators $\widetilde{f}_{n,m,b,\varepsilon}$ defined in (22) and $\widetilde{f}^{N}_{n,m,b,\varepsilon}$ defined in (24), the $LSCV$ function is written as in (26).

4.2.2 Old Faithful data

In this subsection, we consider the well known Old Faithful data given in Tab.2.2 of Silverman (1986). These data consist of the eruption lengths (in minutes) of $107$ eruptions of the Old Faithful geyser in Yellowstone National Park, U.S.A. The data are such that $\min_{i}(x_{i})=1.67$ and $\max_{i}(x_{i})=4.93$, then it is convenient to assume that the density of eruption times is defined on the interval $[1.5,5]$ and transform the data into the unit interval.

The $LSCV$ procedure was performed and resulted in $m=104$ for Vitale’s estimator $\widetilde{f}_{n}$ defined in (1), $(m_{n})=(n^{0.987})$ for our proposed estimator $f_{n}$ defined in (3), $m=66$ for Leblanc’s estimator $\widetilde{f}_{n,m,m/2}$ defined in (17), $m=52$ for the estimator $\widetilde{f}_{n,m,m/4}$ defined in (19), $m=66$ for the multiplicative bias corrected Bernstein estimator $\widetilde{f}_{n,m,2,0.00001}$ defined in (22) and $m=66$ for the normalized estimator $\widetilde{f}^{N}_{n,m,2,0.00001}$ defined in (24). These estimators are shown in Figure 2 along with an histogram of the data and a Gaussian kernel density estimate using the $LSCV$-based bandwidth $h=0.3677$. All the estimators are smooth and seem to capture the pattern highlighted by the histogram. We observe that our recursive estimator outperformed the others estimators near $x=1.5$.

4.2.3 Tuna data

Our last example concerne the tuna data given in Chen (1996). The data come from an aerial line transect survey of Southern Bluefin Tuna in the Great Australian Bight. An aircraft with two spotters on board flies randomly allocated line transects. The data are the perpendicular sighting distances (in miles) of $64$ detected tuna schools to the transect lines. The survey was conducted in summer when tuna tend to stay on the surface. We analyzed the transformed data divided by $w=18$ (the data are such that $\min_{i}(x_{i})=0.19$ and $\min_{i}(x_{i})=16.26$).

The $LSCV$ procedure was performed and resulted in $m=14$ for Vitale’s estimator $\widetilde{f}_{n}$ defined in (1), $(m_{n})=(n^{0.633})$ for our proposed estimator $f_{n}$ defined in (3), $m=4$ for Leblanc’s estimator $\widetilde{f}_{n,m,m/2}$ defined in (17), $m=4$ for the estimator $\widetilde{f}_{n,m,m/4}$ defined in (19), $m=8$ for the multiplicative bias corrected Bernstein estimator $\widetilde{f}_{n,m,2,0.00001}$ defined in (22) and $m=4$ for the normalized estimator $\widetilde{f}^{N}_{n,m,2,0.00001}$ defined in (24). These estimators are shown in Figure 3 along with an histogram of the data and a Gaussian kernel density estimate using the $LSCV$-based bandwidth $h=1.291$. All the estimators are smooth and seem to capture the pattern highlighted by the histogram. We can observe that our proposed recursive estimator outperformed the other estimators near the boundaries.

5 Conclusion

In this paper, we propose a recursive estimator of a density function based on a stochastic algorithm derived from Robbins-Monro’s scheme and using Bernstein polynomials. We first study its asymptotic properties. We show that our proposed estimator of density function have a good boundary properties. Moreover, the bias rate of the proposed estimator is of $m^{-2}$, which is better than the Vitale’s estimator with a bias rate of $m^{-1}$. For almost all the cases, the average $ISE$ of the proposed estimator (3) with a stepsize $\left(\gamma_{n}\right)=\left(n^{-1}\right)$ and the corresponding order $\left(m_{n}\right)$ is smaller than that of Vitale’s estimator and than that of the multiplicative bias-corrected estimator defined by Kakizawa. Furthermore, our proposed recursive density estimator has a slightly larger average $ISE$ compared to Leblanc’s estimator. In addition, a major advantage of our proposal is that its update, when new sample points are available, require less computational cost than the non recursive estimators. Our proposed estimator always integrates to unity, but is not necessarily non negative. However, we found that truncation and renormalisation may solve this issue. Finally, through simple real-life examples (Old Faithful data and tuna data) and a simulation study, we demonstrated how the recursive Bernstein polynomial density estimators can lead to very satisfactory estimates. In conclusion, using the proposed recursive estimator $f_{n}$, we can obtain better results compared to those given by Vitale’s estimator, Leblanc’s estimator and the multiplicative bias-corrected estimator defined by Kakizawa especially near the boundaries.

Appendix A Proofs

In this section, we present proofs for the results presented in Section 3. We need the following technical lemma, which is proved in Mokkadem et al. (2009).

Lemma 1*.*

Let $(v_{n})\in\mathcal{GS}(v^{*}),(\gamma_{n})\in\mathcal{GS}(-\alpha)$, and $l>0$ such that $l-v^{*}\xi>0$. We have

[TABLE]

Moreover, for all positive sequence $(\alpha_{n})$ such that $\lim_{n\rightarrow\infty}\alpha_{n}=0$, and all $\delta\in\mathbb{R}$, we have

[TABLE]

Lemma 1 is widely applied throughout the proofs. Let us underline that it is its application, which requires Assumption $\left(A4\right)$ on the limit of $\left(n\gamma_{n}\right)$ as $n$ goes to infinity. Let us mention that, $\lim_{n\rightarrow\infty}(n\gamma_{n})<\infty$ only if $\alpha=1$, the condition $\lim_{n\rightarrow\infty}(n\gamma_{n})\in(\min\left(a,(2\alpha+a)/4\right),\infty]$ in $\left(A4\right)$, which appears throughout our proofs, is equivalent to the condition $\lim_{n\rightarrow\infty}(n\gamma_{n})\in(\min\left(a,(2+a)/4\right),\infty]$. Similarly, since $\xi\neq 0$ only if $\alpha=1$, we can consider $\alpha=1$ in all the results given in this paper.

A.1 Proof of Proposition 1 and Proposition 2

In view of (3), we have

[TABLE]

Leblanc showed that, for $x\in[0,1]$, we have (see Theorem 6 in Leblanc (2010))

[TABLE]

Substituting this result into (27) leads to

[TABLE]

For $x\in(0,1)$, we obtain

• •

In the case ${a\in\left(0,\frac{2}{9}\alpha\right]}$, we have $\lim_{n\rightarrow\infty}(n\gamma_{n})>2a$. Then the application of lemma 1 gives (4).

• •

In the case ${a\in\left(\frac{2}{9}\alpha,1\right)}$, we have $m^{-2}_{n}=o\left(\sqrt{\gamma_{n}m^{1/2}_{n}}\right)$, since $2a>(2\alpha-a)/4$. Then the application of lemma 1 gives (5).

For $x=0,1$, we obtain

• •

In the case ${a\in\left(0,\frac{\alpha}{5}\right]}$, we have $\lim_{n\rightarrow\infty}(n\gamma_{n})>2a$. Then the application of lemma 1 gives (9).

• •

In the case ${a\in\left(\frac{\alpha}{5},1\right)}$, we have $m^{-2}_{n}=o\left(\sqrt{\gamma_{n}m_{n}}\right)$ and $2a>(\alpha-a)/2$. Then the application of lemma 1 gives (10).

On the other hand, we have

[TABLE]

Leblanc derived that (see Theorem 6 in Leblanc (2010))

[TABLE]

It follows that

[TABLE]

Now, using this result, for $x\in(0,1)$, we obtain

• •

In the case $a\in\left[\frac{2}{9}\alpha,1\right)$, we have $\lim_{n\rightarrow\infty}(n\gamma_{n})>(2\alpha-a)/4$, the application of lemma 1 then gives

[TABLE]

and (6) follows.

• •

In the case $a\in\left(0,\frac{2}{9}\alpha\right)$, we have $2a<(2\alpha-a)/4$ and $\gamma_{n}m_{n}^{1/2}=o\left(m^{-4}_{n}\right)$, the application of lemma 1 ensures that

[TABLE]

which gives (7).

Similarly, for $x=0,1$, we have

• •

In the case $a\in\left[\frac{\alpha}{5},1\right)$, we have $\lim_{n\rightarrow\infty}(n\gamma_{n})>(\alpha-a)/2$, the application of lemma 1 gives (11).

• •

In the case $a\in\left(0,\frac{\alpha}{5}\alpha\right)$, we have $2a<(\alpha-a)/2$ and $\gamma_{n}m_{n}=o\left(m^{-4}_{n}\right)$, the application of lemma 1 gives (12).

A.2 Proof of Proposition 3

Proposition 1 ensures that

• $\bullet$

In the case $a\in\left(0,\frac{2}{9}\alpha\right)$, we have

[TABLE]

• $\bullet$

In the case $a\in\left(\frac{2}{9}\alpha,1\right)$, we have

[TABLE]

Moreover, we note that

[TABLE]

Since we have (see the proof of Theorem 7 in Leblanc (2010))

[TABLE]

then we get

[TABLE]

• $\bullet$

In the case $a\in\left[\frac{2}{9}\alpha,1\right)$, we have $\lim_{n\to+\infty}\left(n\gamma_{n}\right)>(2\alpha-a)/4$. Then the application of lemma 1 gives

[TABLE]

• $\bullet$

In the case $a\in\left(0,\frac{2}{9}\alpha\right)$, we have $\gamma_{n}m_{n}^{1/2}=o\left(m_{n}^{-4}\right)$ and $\lim_{n\to+\infty}\left(n\gamma_{n}\right)>2a$.Then the application of lemma 1 gives

[TABLE]

Part 1 of Proposition 3 follows from the combination of (28) and (31), Part 2 from that of (28) and (30), Part 3 from (29) and (30).

A.3 Proof Corollary 1

Set

[TABLE]

It follows from Proposition 3 that

[TABLE]

If $a=\frac{2}{9}\alpha$, $\left(K_{2}(\xi)m_{n}^{-4}+K_{1}(\xi)\gamma_{n}m_{n}^{1/2}\right)\in\mathcal{GS}\left(-\frac{8}{9}\alpha\right)$ with $\frac{8}{9}\alpha=4a$. If ${a\in\left(0,\frac{2}{9}\alpha\right)}$, $\left(K_{2}(\xi)m_{n}^{-4}\right)\in\mathcal{GS}(-4a)$ with $-4a>-\frac{8}{9}\alpha$, and, if ${a\in\left(\frac{2}{9}\alpha,1\right)}$, $\left(K_{1}(\xi)\gamma_{n}m_{n}^{1/2}\right)\in\mathcal{GS}\left(-(2\alpha-a)/2\right)$ with $-(2\alpha-a)/2>-4a$. It follows that, for a given $\alpha$, to minimize the $MISE\left[f_{n}\right]$, the parameter $a$ must be chosen equal to $\frac{2}{9}\alpha$. Moreover, in view of (32) the parameter $\alpha$ must be equal to 1. We conclude that to minimize the $MISE\left[f_{n}\right]$, the stepsize $(\gamma_{n})$ must be chosen in $\mathcal{GS}(-1)$ and the order $(m_{n})$ in $\mathcal{GS}(2/9)$.

Since the function $x\mapsto K_{2}(\xi)x^{-4}+K_{1}(\xi)\gamma_{n}x^{1/2}$ attains its minimum at the point $\left[\frac{8K_{2}(\xi)}{\gamma_{n}K_{1}(\xi)}\right]^{2/9}$, to minimize $MISE\left[f_{n}\right]$, the order $(m_{n})$ must be equal to $\left(\frac{8K_{2}(\xi)}{K_{1}(\xi)}\right)^{2/9}\gamma_{n}^{-2/9}$. For such a choice, the $MISE$ of $f_{n}$ becomes

[TABLE]

Now, we assume that $\left(\gamma_{n}\right)=\left(\gamma_{0}n^{-1}\right)$ for some $\gamma_{0}\in\left(0,\infty\right)$. In this case, note that $\xi=\lim_{n\to\infty}\left(n\gamma_{n}\right)^{-1}=\gamma_{0}^{-1}$. Then, (14) (hence (1)) is a consequence of the standard trade-off argument.

A.4 Proof Theorem 1

To prove theorem 1, we will use the fact that if $x\in(0,1)$ and $a\in\left[\frac{2}{9}\alpha,1\right]$, then

[TABLE]

which will be proved later. In the case $a\in\left(0,\frac{2}{9}\alpha\right)$, we have

[TABLE]

then if $\gamma^{-1/2}_{n}m^{-9/4}_{n}\rightarrow c$, for some $c\geq 0$, Part 1 of theorem 1 follows from (33). In the case when $a=\frac{2}{9}\alpha$, Parts 1 and 2 of Theorem 1 follow from the combination of (4) and (33). In the case $a\in\left(\frac{2}{9}\alpha,1\right]$, Part 1 of Theorem 1 follows from the combination of (5) and (33).

Now in the case $\gamma_{n}^{-1/2}m_{n}^{-1/4}\rightarrow\infty$ and $a\in\left(0,\frac{2}{9}\alpha\right)$, we have

[TABLE]

Noting that the equation (7) implies

[TABLE]

then, we obtain Part 2 of Theorem 1. We now prove (33). We have

[TABLE]

We set

[TABLE]

The application of Lemma 1 ensures that

[TABLE]

On the other hand, for all $p>0$, we have

[TABLE]

and, since $\lim_{n\rightarrow\infty}(n\gamma_{n})>(\alpha-\frac{a}{2})/2=(2\alpha-a)/4$, there existe $p>0$ such that $\lim_{n\rightarrow\infty}(n\gamma_{n})>\frac{1+p}{2+p}(\alpha-\frac{a}{2})=\frac{(1+p)\alpha-\frac{1+p}{2}a}{p+2}>\frac{(1+p)\alpha-(2+p)a}{p+2}$, then the application of lemma 1 gives

[TABLE]

hence

[TABLE]

Then the convergence in (33) follows from the application of Lyapounov’s Theorem.

Funding

This work benefited from the financial support of the GDR 3477 GeoSto.

Acknowledgments

The authors would like to thank the Editor and the referees for their very helpful comments, which led to considerable improvement of the original version of the paper and a more sharply focused presentation.