arXiv:1606.06988v1 [] 22 Jun 2016

Recursive kernel density estimators under missing data
Yousri Slaoui
Université de Poitiers, Laboratoire de Mathématiques et Application,
Futuroscope Chasseneuil, France
E.mail: Yousri.Slaoui@math.univ-poitiers.fr

Abstract: In this paper we propose an automatic bandwidth selection of the recursive kernel density estimators with missing data in the context of global and local density estimation. We showed
that, using the selected bandwidth and a special stepsize, the proposed recursive estimators outperformed the nonrecursive one in terms of estimation error in the case of global estimation. However,
the recursive estimators are much better in terms of computational costs. We corroborated these
theoretical results through simulation studies and on the simulated data of the Aquitaine cohort of
HIV-1 infected patients and on the coriell cell lines using the chromosome number 11.
2010 Mathematics Subject Classification: Primary 62G07, 62L20; Secondary 65D10, 62N01.
Key words and phrases: Density estimation; Stochastic approximation algorithm; Smoothing, curve
fitting; missing data.

1

Introduction

Kernel density estimator (KDE) were introduced in Parzen (1962) and Rosenblatt (1956) have been
widely used in many applications research including clinical trials, epidemiology, genetics as an
exploratory tool.
There has been an intensive work in the literature about the bandwidth selection methods of this
estimator. See for example Duin (1976), Rudemo (1982), Scott and Terrell (1987) and Marron
(1988). However, in many situation the full data are not available and are subjects to missing
data. Several methods have been used to overcome the problem and can be found in the literature.
For example, the EM algorithm introduced in Dempster et al. (1977), can be used to obtain maximum likelihood estimates when some data are missing, and multiple imputation, which is a Monte
Carlo technique in which the missing values are replaced by simulated versions, and the results are
combined to produce estimates which incorporate missing data, see the paper of Little and Rubin
(2002). Note that all of these methods are concerned with the estimation of unknown parameters.
In the framework of nonparametric estimation, Dubnicka (2009) used a modified Horvitz-Thompsontype KDE to estimate the global density under missing data. However, several methods of local
density estimation was widely discussed and extended in many directions; see Sheather (1983,
1986), Thombs and Sheather (1992), Hazelton (1999), Chan et al. (2010) and recently Dutta (2014).
In this paper, we propose a Recursive KDE to estimate the global unknown density f of a univariate variable when the data are missing at random, we also give some comparative elements in
the case of local density estimation. Then, we present a data-driven bandwidth selection algorithm
for the proposed estimators. Data-driven bandwidth selection procedure was proposed by Slaoui
(2014a) in the framework of the recursive kernel density estimators, then, Slaoui (2014b) propose a
plug-in selection method for recursive kernel distribution estimators, Slaoui (2015) propose a plug-in
1

selection method for recursive kernel regression estimators with a fixed design setting, and Slaoui
(2016) propose a plug-in selection algorithm for the semi-recursive kernel regression estimators, all of
this works suppose that the full data are observed. Here, we developed a specific second-generation
plug-in bandwidth selection method of the Recursive KDE under missing data.
Interestingly, unlike the EM-algorithm, the proposed Recursive KDE has explicit form, which
facilitates its usage and analysis. Moreover, the CPU time using the proposed recursive KDE are
approximately two times faster than the CPU time using the nonrecursive Horvitz-Thompson-type
KDE.
For further motivation in order to help intuition, let us introduce in some detail two datasets,
the first one concerne a simulated data of the Aquitaine cohort of HIV-1 infected patients, see the
paper of Thiébault et al. (2000), the second one concerne a real data which correspond to two
array Comparative Genomic Hybridization (CGH) studies of fibroblast cell strains. Assuming that
all these measurement are corrects, we simulated various proportions of missing at random by
removing the corresponding entries before performing our algorithm. The estimated densities are
then compared to the shape of the histograms obtained using the full data. Even when 30% of the
original measurements are missing, the estimated densities using the proposed estimators remain
very accurate thus demontrating the effectiveness of our approach. Moreover, we compared our
approach with the nonrecursive Horvitz-Thompson-type KDE via real dataset as well as simulations.
Results showed that our approach outperformed other approaches in terms of estimation accuracy
and computing efficiency. The remainder of the paper is organized as follows. In the next section we
present our proposed recursive KDE. In Section 3, we state our main results. Section 4 is devoted to
our application results, first by simulation (subsection 4.1) and second using two datasets, the first
one concerne the simulated data of the Aquitaine cohort of HIV-1 infected patients and the second
one concerne the coriell cell lines using the chromosome number 11 (subsection 4.2). We conclude
the paper in Section 5, wehereas the technical details are deferred to Appendix A and Appendix B
gives the results for the Horvitz-Thompson-type KDE.

2

The estimator

Let T, T1 , . . . , Tn be independent, identically distributed random variables, and let f and F denote
respectively the probability density and the distribution of T . We suppose that the full data
T1 , . . . , Tn are not totally available and are subjects to missing data. The observed random variables
are then Xi and δi where
Xi = δi × Ti and δi = 1{Ti is observed} , 1 ≤ i ≤ n.
When some Ti are missing, we let : πi = P [δi = 1|Ti ] for 1 ≤ i ≤ n.
This probability is often called the propensity score, see Rosenbaum and Rubin (1983). Let us
recall that, in order to construct a stochastic algorithm, which approximates the function f at a
given point x, we need to define an algorithm of search of the zero of the function h : y → f (x) − y.
Following Robbins-Monro’s procedure, this algorithm is defined by setting f0 (x) ∈ R, and, for all
n ≥ 1,
fn (x) = fn−1 (x) + γn Wn ,
where Wn (x) is an observation of the function h at the point fn−1 (x), and the stepsize (γn ) is
a sequence of positive real numbers that goes to zero. In order to define Wn (x), we follow the
2

approach of Révész (1973, 1977), Tsybakov
R (1990) and of Slaoui (2013, 2014a,b), and we introduce
a kernel K (that is, a function satisfying R K(x)dx = 1), and a bandwidth (hn ) (that is, a sequence
−1
−1
of positive real numbers that goes to zero), and sets Wn (x) = h−1
n δn πn K hn [x − Xn ] − fn−1 (x).
Then, the estimator fn to recursively estimate the function f at the point x can be written as

−1
fn (x) = (1 − γn ) fn−1 (x) + γn δn πn−1 h−1
(1)
n K hn [x − Xn ] ,

where the stepsize (γn ) is a sequence of positive
Q real numbers that goes to zero. Let us underline
that we consider f0 (x) = 0 and we let Qn = nj=1 (1 − γj ), then the equation (1) can be rewritten
as follows:


n
X
x − Xk
−1 −1
,
fn (x) = Qn
Q−1
γ
δ
π
h
K
k k k
k
k
hk
k=1

which means that, we can estimate f recursively at the point x by using one of the two previous
equations.
R
Moreover, we show in the next section that the optimal bandwidth which minimize the E R [fn (x) − f (x)]2 dx
depend on the choice of the stepsizes (γn ); we show in
 particular that under some conditions of
regularity of f and using the stepsizes (γn ) = γ0 n−1 , with γ0 > 2/5, the bandwidth (hn ) must
equal


)1/5
)1/5 ( R
(
R 2
2
R K (z) dz
R f (x) dx
2−1/5 (γ0 − 2/5)1/5 R
πn−1/5 n−1/5  .

2
R
2 K (z) dz 2
(2)
z
f
(x)
f
(x)
dx
R
R
The first aim of the next section is to propose a data-driven bandwidth selection of the proposed
recursive KDE in the case of missing data, and the second aim is to give the conditions under which
the proposed estimators fn may behave more efficiently than the nonrecursive Horvitz-Thompsontype KDE see the papers Horvitz and Thompson (1952) and Dubnicka (2009), defined as


n
x − Xi
1 X
−1
e
δi πi K
.
(2)
fn (x) =
nhn
hn
i=1

In the next section we state our main results: the global density estimation is considered in Section 3.1 and the local density estimation is developed in Section 3.2.

3

Assumptions and main results

We define the following class of regularly varying sequences.
Definition 1. Let γ ∈ R and (vn )n≥1 be a nonrandom positive sequence. We say that (vn ) ∈ GS (γ)
if


vn−1
= γ.
(3)
lim n 1 −
n→+∞
vn
Condition (3) was introduced Galambos and Seneta (1973) to define regularly varying sequences,
see also the paper Bojanic and Seneta (1973) and by Mokkadem and Pelletier (2007) in the context
of stochastic approximation algorithms. Noting that the acronym GS stand for (Galambos and
Seneta). Typical sequences in GS (γ) are, for b ∈ R, nγ (log n)b , nγ (log log n)b , and so on.
In this section, we investigate asymptotic properties of the proposed estimators (1). The assumptions to which we shall refer are the following
3

(A1) K : RR → R is a continuous, bounded function satisfying
and R z 2 K (z) < ∞.

R

R

K (z) dz = 1, and,

R

R zK

(z) = 0

(A2) i) (γn ) ∈ GS (−α) with α ∈ ]1/2, 1].
ii) (hn ) ∈ GS (−a) with a ∈ ]0, 1[.
iii) limn→∞ (nγn ) ∈ ]min {2a, (α − a) /2} , ∞].

(A3) f is bounded, differentiable, and f (2) is bounded.
Assumption (A2) (iii) on the limit of (nγn ) as n goes to infinity is usual inthe framework
of

−1
is finite.
stochastic approximation algorithms. It implies in particular that the limit of [nγn ]
For simplicity, we introduce the following notations:
lim (nγn )−1 ,
Z
K 2 (z) dz,
R (K) =
R
Z
z j K (z) dz,
µj (K) =
ξ =

(4)

n→∞

R

Θ (K) = R (K)4/5 µ2 (K)2/5 ,
Z
f 2 (x) dx,
I1 =
ZR 
2
f (2) (x) f (x) dx.
I2 =
R

In this section, we explicit the choice of (hn ) through a second-generation plug-in method, which
minimize the Mean Weighted Integrated Squared Error M W ISE of the proposed recursive estimators (1) in the case of global estimation and through a second-generation plug-in method, which
minimize the Mean Squared Error M SE of the proposed recursive estimators (1) in the case of
local estimation, in order to provide a comparison with the nonrecursive estimator (2). Our first
result is the following proposition, which gives the bias and the variance of fn .
Proposition 1 (Bias and variance of fn ). Let Assumptions (A1) − (A3) hold
1. If a ∈]0, α/5], then
E [fn (x)] − f (x) =
If a ∈]α/5, 1[, then


h2n
f (2) (x) µ2 (K) + o h2n .
2 (1 − 2aξ)

E [fn (x)] − f (x) = o

q

γn h−1
n



.

2. If a ∈ [α/5, 1[, then
γn −1
1
V ar [fn (x)] =
πn
f (x) R (K) + o
hn
(2 − (α − a) ξ)
If a ∈]0, α/5[, then


V ar [fn (x)] = o h4n .
4

(5)

(6)



γn
hn



.

(7)

(8)

3. If limn→∞ (nγn ) > max {2a, (α − a) /2}, then (5) and (7) hold simultaneously.
The bias and the variance of the proposed estimator fn defined by the stochastic approximation
algorithm (1) then heavily depend on the choice of the stepsize (γn ). Let us first underline that, it
follows from (7) that the stepsize which minimize the variance of fn is (γn ) = [1 − a] n−1 , using
this stepsize the variance of fn is equal to


1
1−a
f (x) R (K) + o
.
V ar [fn (x)] = πn−1
nhn
nhn

Now, using the special stepsize (γn ) = n−1 , (see Mokkadem et al. (2009) and Slaoui (2013)), the
variance of fn is equal to


1
π −1 1
f (x) R (K) + o
.
V ar [fn (x)] = n
1 + a nhn
nhn
Let us recall that under the Assumptions (A1), (A2) ii) and (A3), we have


h
i
1
πn−1
e
f (x) R (K) + o
.
V ar fn (x) =
nhn
nhn

Which shows that performing the proposed recursive estimators (1) with one of two proposed
stepsizes we get smaller variance than the nonrecursive estimator (2). Similar results was given
by Mokkadem et al. (2009) and Slaoui (2013) in the case of complet data.
Let us now state the following theorem, which gives the weak convergence rate of the estimator
fn defined in (1).
Theorem 1 (Weak pointwise convergence rate). Let Assumptions (A1) − (A3) hold
1. If there exists c ≥ 0 such that γn−1 h5n → c, then
q
γn−1 πn hn (fn (x) − f (x))


√
1
c
D
(2)
→N
f (x) µ2 (K) ,
f (x) R (K) .
2 (1 − 2aξ)
(2 − (α − a) ξ)
2. If γn−1 h5n → ∞, then
1
1
P
(fn (x) − f (x)) →
f (2) (x) µ2 (K) ,
h2n
2 (1 − 2aξ)
P

D

where → denotes the convergence in distribution, N the Gaussian-distribution and → the convergence in probability.

3.1

Global density estimation

In order to measure globally the quality of our recursive estimators (1), we use the following quantity,
Z
M W ISE [fn ] = E [fn (x) − f (x)]2 f (x) dx
Z
Z R
2
V ar (fn (x)) f (x) dx.
(E (fn (x)) − f (x)) f (x) dx +
=
R

R

The following proposition gives the M W ISE of the recursive estimators defined in (1).
5

Proposition 2 (M W ISE of fn ). Let Assumptions (A1) − (A3) hold.
1. If a ∈]0, α/5[, then
M W ISE [fn ] =


1
h4n
I2 µ22 (K) + o h4n .
2
4 (1 − 2aξ)

2. If a = α/5, then
M W ISE [fn ] =


1
γn πn−1
1
h4
I2 µ22 (K) + o h4n .
I1 R (K) + n
2
hn (2 − (α − a) ξ)
4 (1 − 2aξ)

3. If a ∈]α/5, 1[, then
M W ISE [fn ] =
3.1.1

1
γn πn−1
I1 R (K) + o
hn (2 − (α − a) ξ)



γn πn−1
hn



.

Estimating propensity scores

In order to estimate the propensity scores πi , we need to exploit the information contained in the
auxiliary variables Ti , which are related to Xi . First, we can estimate π by the empirical proportion
based on the observed data (Qi et al. (2005)):
Pn
j=1 δj 1{Xj =Xi }
.
π
bi = π
b (Xi ) = Pn
j=1 1{Tj =Ti }
Moreover, in the case when Xi contains continuous elements, Dubnicka (2009) used the NadarayaWatson (local mean) estimator (Nadaraya (1964); Watson (1964)):

Pn
−1
j=1 δj K hn [Xi − Xj ]
 ,
π
bN W i = π
bN W (Xi ) = Pn
(9)
−1
j=1 K hn [Xi − Xj ]

and in the case of binary nature of the response variable, Tibshirani and Hastie (1987) and Fan et al.
(1995) provide an estimation of the propensity score using local likelihood. However, Dubnicka
(2009) concluded that, both the Nadaraya-Watson estimates and local likelihood estimates of the
propensity scores are more flexible than ordinary binary regression estimates. In this paper we use
the recursive version of the estimator of Nadaraya-Watson (called also the semi-recursive estimator)
and defined as:


Pn
−1
−1
δ
h
K
h
[X
−
X
]
i
j
j=1 j j
j
 .
(10)
π
bRN W i = π
bRN W (Xi ) = Pn
−1
−1
j=1 hj K hi [Xi − Xj ]
For now, we simply use the propensity score estimator (9) in the case of nonrecursive HorvitzThompson-type KDE and (10) in the case of the proposed recursive KDE.
3.1.2

Bandwidth selection

In the framework of the nonparametric kernel estimators, the bandwidth selection methods studied
in the literature can be divided into three broad classes: the cross-validation techniques, the plug-in
ideas and the bootstrap. A detailed comparison of the three practical bandwidth selection can
6

be found in Delaigle and Gijbels (2004). They concluded that chosen appropriately plug-in and
bootstrap selectors both outperform the cross-validation bandwidth, and that none of the two can
be claimed to be best in all cases. In this section, we developed a plug-in bandwidth selector
that minimizing the M W ISE of the proposed recursive KDE, using the function f (x) as a weight
function.
The following corollary ensures that the bandwidth which minimize the M W ISE of fn depend on
the stepsize (γn ) and then the corresponding M W ISE depend also on the stepsize (γn ).
Corollary 1. Let Assumptions (A1) − (A3) hold. To minimize the M W ISE of fn , the stepsize (γn )
must be chosen in GS (−1), the bandwidth (hn ) must equal

(
)1/5 
1/5
2
I1
R (K)
(1 − 2aξ)

πn−1/5 γn1/5  .
(2 − (α − a) ξ) I2
µ22 (K)
Then, we have



5
4/5 1/5
(1 − 2aξ)−2/5 (2 − (α − a) ξ)−4/5 I1 I2 Θ (K) πn−4/5 γn4/5 + o πn−4/5 γn4/5 .
4

The following corollary shows that, for a special choice of the stepsize (γn ) = γ0 n−1 , which
fulfilled that limn→∞ nγn = γ0 and that (γn ) ∈ GS (−1), the optimal value for hn depend on γ0 and
then the corresponding M W ISE depend on γ0 .

Corollary 2. Let Assumptions (A1) − (A3) hold, and suppose that (γn ) = γ0 n−1 . To minimize
the M W ISE of fn , the stepsize (γn ) must be chosen in GS (−1), limn→∞ nγn = γ0 , the bandwidth
(hn ) must equal
!
 1/5 

R (K) 1/5 −1/5 −1/5
1/5 I1
−1/5
,
2
(γ0 − 2/5)
πn n
I2
µ22 (K)
M W ISE [fn ] =

and we then have
M W ISE [fn ] =



5 1
γ02
4/5 1/5
−4/5 −4/5
−4/5 −4/5
.
I
I
Θ
(K)
π
n
+
o
π
n
n
n
4 24/5 (γ0 − 2/5)6/5 1 2

Moreover, the minimum of γ02 (γ0 − 2/5)−6/5 is reached at γ0 = 1, then the bandwidth (hn ) must
equal
!
 1/5  1/5 

3
R (K) 1/5 −1/5 −1/5
I1
,
(11)
πn n
10
I2
µ22 (K)
and then the asymptotic M W ISE
AM W ISE [fn ] =

5 1
4 24/5

 6/5
5
4/5 1/5
I1 I2 Θ (K) πn−4/5 n−4/5 .
3

In order to estimate the optimal bandwidth (11), we must estimate the unknown quantites I1
and I2 . For this purpose, we used the following modified version of the kernel estimators introduced
in Slaoui (2014a):

7

Ib1 =

Ib2 =



n
Xi − Xk
Ψn X −1
−1
K
δ
π
b
Ψk βk b−1
k k RN W k b
n
bk

(12)

i,k=1

Φ2n
n

n
X

i,j,k=1
j6=k

−1
′−3 ′−3 (2)
−1 ′ ′
−1
bRN
bRN
Φ−1
j Φk βj βk δj δk π
W jπ
W k bj bk Kb′

Xi − Xj
b′j

!

(2)
Kb′




Xi − Xk
(13)
b′k

where Kb and Kb′ are a kernels, bn and b′n are respectively the associated bandwidth (called pilot
two pilot stepsizes for the estimation of I1 and I2 respectively,
bandwidth)
βn and βn′ are theQ
Qand
n
and Ψn = i=1 (1 − βi ) and Φn = ni=1 (1 − βi′ ).
In practice, we take


Q3 − Q1
−β
, β ∈ ]0, 1[
(14)
bn = n min sb,
1.349

(see Silverman (1986)) with sb the sample standard deviation, and Q1 , Q3 denoting the first and
third quartiles, respectively.
We followed the same steps as in Slaoui (2014a) and we showed that in order to minimize the
AM ISE of Ib1 the pilot bandwidth
(bn ) must belong to GS (−2/5) and the pilot stepsize (βn )

should be equal to 1.36n−1 , and in order to minimize the AM ISE of Ib2 the pilot bandwidth
(b′n )

must belong to GS (−3/14) and the pilot stepsize (βn′ ) should be equal to 1.48n−1 .
Finally, the plug-in estimator
 of the bandwidth (hn ) using the recursive estimators defined in (1)
−1
with the stepsizes (γn ) = n
is equal to


1/5
 1/5 b !1/5 
R
(K)
I
3
−1/5
1

π
bRN W n n−1/5  ,
(15)
2 (K)
b
10
µ
I2
2

and the associated plug-in AM W ISE is equal to

 6/5
5
5 1
4/5 1/5
−4/5
Ib1 Ib2 Θ (K) π
bRN W n n−4/5 .
4 24/5 3

Now, using the stepsize (γn ) = [4/5] n−1 , the stepsize which minimize the variance of the proposed
estimators defined in (1), it follows from Corollary 2 that, the plug-in estimator of the bandwidth
(hn ) in this particular case is equal to


!1/5 
1/5
b
R (K)
−1/5
5−1/5 I1
π
bRN W n n−1/5  ,
(16)
2 (K)
b
µ
I2
2
AM\
W ISE [fn ] =

and the associated plug-in AM W ISE is equal to

4/5 1/5
−4/5
AM\
W ISE [fn ] = 51/5 Ib1 Ib2 Θ (K) π
bRN W n n−4/5 .

Moreover, following similar steps as in Slaoui (2014a), we prove the following corollary

8

Corollary 3. Let theassumptions (A1) − (A3) hold, and the bandwidth (hn ) equal to (15) and the
stepsize (γn ) = n−1 . We have
h
i
E AM\
W ISE (fn )
h
 i < 1.
E AM^
W ISE fen

Then, the expectation of the estimated AM WISE of the proposed recursive estimators defined
by (1) using the special stepsize (γn ) = n−1 is smaller than the expectation of the estimated
AM W ISE of the nonrecursive Horvitz-Thompson-type KDE defined by (2) (see Appendix B for
some results on the nonrecursive Horvitz-Thompson-type KDE estimator.

3.2

Local density estimation

In order to measure locally the quality of our recursive estimators (1), we use the following quantity,
M SE [fn ] = E [fn (x) − f (x)]2

= (E (fn (x)) − f (x))2 + V ar (fn (x)) .

Following similar step as in the case of the global density estimation, we show that the value of hn
which minimizes the asymptotic mean square error of fn (x) is equal to


!1/5 
1/5
 1/5
R (K)
f (x)
 3
πn−1/5 n−1/5  .
(17)
2
2 (K)
(2)
10
µ
f (x)
2

In practice the kernel K is known but f (x) and f (2) (x) are not. Thus, we estimate f (x) by fn (x)
(2)
and f (2) (x) by fn (x), and so the plug-in bandwidth selection in the case of locally estimation is
equal to



1/5
 1/5


 3
R (K) 1/5 −1/5 −1/5 
 fbn (x) 

,
(18)
π
bRN W n n




2 (K)
2
 10

(2)
µ
2
fbn (x)

where




x − Xk
bk
k=1


n
X
−1
′−3 (2) x − Xk
−1 ′
(2)
b
fn (x) = Φn
bRN W k bk Kb′
Φk βk δk π
b′k
fbn (x) = Ψn

n
X

−1
−1
bRN
Ψ−1
W k Kb
k βk bk δk π

k=1

where Kb and Kb′ are a kernels, bn and b′n are respectively the associated bandwidth (called pilot
′
bandwidth) and βn and
stepsizes for the estimation of f (x) and f (2) (x) reQn βn are the two pilot Q
spectively, and Ψn = i=1 (1 − βi ) and Φn = ni=1 (1 − βi′ ).

9

4

Applications

The aim of our applications is to compare the performance of the recursive kernel density estimators
under missing data defined in (1) with that of the nonrecursive Horvitz-Thompson-type KDE defined
in (2).
When applying fn one need to choose four quantities:
1. The function K, we choose the Normal kernel.

2. The stepsizes (γn ) equal to c n−1 , with c ∈ [4/5, 1].

3. The propensity score (πn ) is chosen to be equal to (10).


4. The bandwidth (hn ) is chosen to be equal to (15) in the case when (γn ) = n−1 and
to (16) in the case when (γn ) = [4/5] n−1 .

(a) To estimate I1 , we use (12); The pilot bandwidth is chosen to be equal to (14) with
the choice of β = 2/5 and the pilot stepsize equal to 1.36n−1 .
(b) To estimate I2 , we use (13); The pilot bandwidth is chosen to be
 equal to (14) with
−1
the choice of β = 3/14 and the pilot stepsize equal to 1.48n
.

When applying fen one need to choose three quantities:

1. The function K, as in the recursive framework, we use the Normal kernel.
2. The propensity score (πn ) is chosen to be equal to (9).
3. The bandwidth (hn ) is chosen to be equal to (29).
(a) To estimate I1 , we use (27); The pilot bandwidth is chosen to be equal to (14) with
the choice of β = 2/5.
(b) To estimate I2 , we use (28); The pilot bandwidth is chosen to be equal to (14) with
the choice of β = 3/14.

4.1
4.1.1

Simulations
Global density estimation

In order to investigate the comparison between the proposed estimators, we consider three sample
sizes: n = 100, 200, and 500. In each case, we consider 0%, 30%, 50% and 70% of missing data,
the number of the design points was fixed to be equal to 500, which variate from the lowest value
to highest value over (N = 500 number of simulations) with equally spaced setting. We then use
the bandwidth selection method proposed in the applications Section. Moreover, we consider three
densities functions f : 1- the standard normal : X ∼ N (0, 1) (see Table 1), 2- the normal mixture
distribution : X ∼ 1/2N (2, 1) + 1/2N (−3, 1) (see Table 2), 3- the weibull distribution with shape
parameter 2 and scale parameter 1: X ∼ Weibul (2, 1) (see Table 3). We considered also the case
of one auxiliary variable and generated the pairs (X, Y ) , (X1 , Y1 ) , . . . , (Xp
n , Yn ) by first generating
Y distributed as FY and X ∗ distributed as FX ∗ , and we let X = ρY + 1 − ρ2 X ∗ , −1 < ρ < 1,
then, Y and X were correlated. Moreover, we generated X ∗ , X1∗ , . . . , Xn∗ from the standard normal
distribution and we set ρ = 0.3, 0.5 and 0.8 and we generated the responses Y, Y1 , . . . , Yn from
one of the three previous densities functions f (see Table 4), in this case we present results for the
situations in which 70% of the responses are missing. For each situations, we compute the Mean
Weighted Integrated Squared Error M W ISE (over 500 samples).

10

0%
Nonrecursive
Recursive 1
Recursive 2
30%
Nonrecursive
Recursive 1
Recursive 2
50%
Nonrecursive
Recursive 1
Recursive 2
70%
Nonrecursive
Recursive 1
Recursive 2

n = 100
MWISE CPU
3.86e−03
364
3.70e−03 194
3.85e−03
184
MWISE CPU
1.47e−02
365
1.27e−02 164
1.32e−02
167
MWISE CPU
2.23e−02
356
−02
2.11e
169
2.17e−02
176
MWISE CPU
2.89e−02
346
−02
2.63e
168
2.77e−02
177

X ∼ N (0, 1)
n = 200
MWISE CPU
2.33e−03
2448
2.30e−03 1225
2.39e−03 1169
MWISE CPU
5.93e−03
2448
5.73e−03 1206
6.17e−03
1216
MWISE CPU
1.87e−02
2567
−02
1.69e
1146
1.85e−02
1154
MWISE CPU
2.43e−02
2267
−02
2.32e
1166
2.39e−02
1175

n = 500
MWISE
CPU
−03
1.20e
14432
1.20e−03 7118
1.25e−03
7122
MWISE
CPU
2.59e−03 14134
2.56e−03 7112
2.70e−03 7043
MWISE
CPU
3.32e−03 14123
3.24e−03 7089
3.46e−03
7134
MWISE
CPU
4.82e−03 14134
4.63e−03 7156
4.86e−03 7143

Table 1: Quantitative comparison between the nonrecursive Horvitz-Thompson-type KDE (2) and
two recursive estimators; recursive 1 correspond to the proposed estimator (1) with the choice of

(γn ) = n−1 , resursive 2 correspond to the estimator (1) with the choice of (γn ) = [4/5] n−1 .
Here we consider the normal distribution X ∼ N (0, 1) with the proportion of missing at random
equal respectively to 0% (in the first block), equal to 10% (in the second block), equal to 20% (in the
third block) and equal to 30% (in the last block), we consider three sample sizes n = 100, n = 200
and n = 500, the number of simulations is 500, and we compute the Mean Weighted Integrated
Squared Error (M W ISE) and the CPU time in seconds.
Computational cost The advantage of recursive estimators on their nonrecursive 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 (1) that for all n1 ∈ [0, n − 1],


n−1
n
n
X
Y
Y

−1
−1

(1 − γj ) fn1 (x) +
(1 − γj ) γk h−1
fn (x) =
k δk πk K hk (x − Xk ) .
j=n1 +1

k=n1

j=k+1

In order to give some comparative elements with nonrecursive Horvitz-Thompson-type KDE (2),
including computational costs. We consider a 500 samples of size n1 = ⌊n/2⌋ (the lower integer part
of n/2) generated from respectively the five considered distributions, moreover, we suppose that
we receive an additional 500 samples of size n − n1 generated also from the same five considered
distributions. Performing the two methods, we report the total CPU time values for each considered
distribution and in all cases given in Tables 1, 2, 3 and 4, the CPU time is given in seconds.

11

0%
Nonrecursive
Recursive 1
Recursive 2
30%
Nonrecursive
Recursive 1
Recursive 2
50%
Nonrecursive
Recursive 1
Recursive 2
70%
Nonrecursive
Recursive 1
Recursive 2

X ∼ 21 N (2, 1) + 12 N (−3, 1))
n = 200
n = 500
n = 100
MWISE CPU MWISE CPU MWISE
CPU
9.96e−04
372
6.09e−04
2448 3.38e−04 14104
8.64e−04
174 5.50e−04 1168 3.09e−04 7016
−04
9.00e
177
5.72e−04
1189 3.23e−04 7003
MWISE CPU MWISE CPU MWISE
CPU
−03
−04
−04
1.30e
366
8.31e
2449 3.90e
14142
1.25e−03 175 7.69e−04 1246 3.82e−04 7215
8.01e−04 1238 3.97e−04
7225
1.30e−03
168
MWISE CPU MWISE CPU MWISE
CPU
1.41e−03
354
1.55e−03
2246 8.07e−04 13994
−03
−03
1.29e
166 1.31e
1185 6.52e−04 6946
1.36e−03 1174 6.79e−04
6987
1.35e−03
173
MWISE CPU MWISE CPU MWISE
CPU
1.95e−03
356
1.87e−03
2347 1.13e−03 13965
1.80e−03 168 1.69e−03 1164 1.09e−03 7064
1.72e−03 1162 1.11e−03
7086
2.02e−03
179

Table 2: Quantitative comparison between the nonrecursive Horvitz-Thompson-type KDE (2) and
two recursive estimators; recursive 1 correspond to the proposed estimator (1) with the choice of

(γn ) = n−1 , resursive 2 correspond to the estimator (1) with the choice of (γn ) = [4/5] n−1 .
Here we consider the normal mixture distribution X ∼ 12 N (2, 1) + 21 N (−3, 1), with the proportion
of missing at random equal respectively to 0% (in the first block), equal to 30% (in the second
block), equal to 50% (in the third block) and equal to 70% (in the last block), we consider three
sample sizes n = 100, n = 200 and n = 500, the number of simulations is 500, and we compute the
Mean Weighted Integrated Squared Error (M W ISE) and the CPU time in seconds.

12

0%
Nonrecursive
Recursive 1
Recursive 2
30%
Nonrecursive
Recursive 1
Recursive 2
50%
Nonrecursive
Recursive 1
Recursive 2
70%
Nonrecursive
Recursive 1
Recursive 2

n = 100
MWISE CPU
1.82e−02
374
1.70e−02
164
1.77e−04
178
MWISE CPU
2.09e−02
366
1.93e−03 177
2.01e−03
167
MWISE CPU
2.25e−02
356
−03
2.14e
166
2.19e−02
168
MWISE CPU
2.45e−02
364
−02
2.36e
174
2.41e−02
166

X ∼ Weibull (2, 1)
n = 200
MWISE CPU
1.11e−02
2357
1.09e−02 1184
1.15e−02
1194
MWISE CPU
1.22e−02
2452
1.16e−04 1154
1.23e−04 1124
MWISE CPU
1.65e−02
2354
−02
1.59e
1163
1.62e−02
1165
MWISE CPU
2.12e−02
2346
−02
2.03e
1187
2.08e−02 1184

n = 500
MWISE
CPU
−03
4.75e
14256
4.74e−04 7122
4.94e−04 7089
MWISE
CPU
−03
5.19e
14327
5.12e−03 7012
5.27e−03 6984
MWISE
CPU
1.14e−02 14232
1.07e−02 7048
1.12e−02 7044
MWISE
CPU
1.82e−02 13982
1.72e−02 7086
1.83e−02
7094

Table 3: Quantitative comparison between the nonrecursive Horvitz-Thompson-type KDE (2) and
two recursive estimators; recursive 1 correspond to the proposed estimator (1) with the choice of

(γn ) = n−1 , resursive 2 correspond to the estimator (1) with the choice of (γn ) = [4/5] n−1 .
Here we consider the weibull distribution with shape parameter 2 and scale parameter 1, X ∼
Weibull (2, 1), with the proportion of missing at random equal respectively to 0% (in the first
block), equal to 30% (in the second block), equal to 50% (in the third block) and equal to 70%
(in the last block), we consider three sample sizes n = 100, n = 200 and n = 500, the number of
simulations is 500, and we compute the Mean Weighted Integrated Squared Error M W ISE and
the CPU time in seconds.

13

0.4
0.2
0.0

0.1

Estimated Density

0.3

Nonrecursive
Recursive1
Recursive2

−2

−1

0

1

2

Figure 1: Qualitative comparison between the nonrecursive Horvitz-Thompson-type KDE (2) and
two recursive estimators; recursive 1 correspond to the proposed estimator (1) with the choice of

(γn ) = n−1 , resursive 2 correspond to the estimator (1) with the choice of (γn ) = [4/5] n−1 .
Here we consider the normal distribution N (0, 1) under 70% of the missing data and we consider
500 samples of size n = 200.

14

0.20
0.10
0.05

Estimated Density

0.15

Nonrecursive
Recursive1
Recursive2

−4

−2

0

2

4

Figure 2: Qualitative comparison between the nonrecursive Horvitz-Thompson-type KDE (2) and
two recursive estimators; recursive 1 correspond to the proposed estimator (1) with the choice of

(γn ) = n−1 , resursive 2 correspond to the estimator (1) with the choice of (γn ) = [4/5] n−1 .
Here we consider the normal mixture distribution X ∼ 12 N (2, 1) + 21 N (−3, 1) under 70% of the
missing data and we consider 500 samples of size n = 200.

15

0.4
0.0

0.2

Estimated Density

0.6

0.8

Nonrecursive
Recursive1
Recursive2

0.0

0.5

1.0

1.5

2.0

Figure 3: Qualitative comparison between the nonrecursive Horvitz-Thompson-type KDE estimator (2) and two recursive estimators;
recursive 1 correspond to the proposed estimator (1)

with the choice of (γn ) = n−1 , resursive 2 correspond to the estimator (1) with the choice of
(γn ) = [4/5] n−1 . Here we consider the weibull distribution with shape parameter 2 and scale
parameter 1, X ∼ Weibull (2, 1) under 70% of the missing data and we consider 500 samples of size
n = 200.

16

n = 100
ρ = 30%
Nonrecursive
Recursive 1
Recursive 2
ρ = 50%
Nonrecursive
Recursive 1
Recursive 2
ρ = 80%
Nonrecursive
Recursive 1
Recursive 2

MWISE
0.02882
0.02622
0.02746
MWISE
0.02866
0.02634
0.02724
MWISE
0.02894
0.02664
0.02784

ρ = 30%
Nonrecursive
Recursive 1
Recursive 2
ρ = 50%
Nonrecursive
Recursive 1
Recursive 2
ρ = 80%
Nonrecursive
Recursive 1
Recursive 2

MWISE
0.00328
0.03014
0.00329
MWISE
0.00332
0.00318
0.00334
MWISE
0.00334
0.00324
0.00338

ρ = 30%
Nonrecursive
Recursive 1
Recursive 2
ρ = 50%
Nonrecursive
Recursive 1
Recursive 2
ρ = 80%
Nonrecursive
Recursive 1
Recursive 2

MWISE
0.02342
0.02265
0.02348
MWISE
0.02354
0.02289
0.02358
MWISE
0.02372
0.02314
0.02369

n = 200
n = 500
Y ∼ N (0, 1)
CPU MWISE
CPU MWISE CPU
357
0.02584
2335
0.05042 14342
0.0244
1244 0.0484
7332
187
194
0.02542
1255
0.05044
7412
CPU MWISE
CPU MWISE CPU
366
0.02564
2534
0.05054 14542
0.02454
1254 0.04864 7353
191
187
0.02544
1243 0.05124
7366
CPU MWISE
CPU MWISE CPU
362
0.02586
2447
0.05084 14121
0.02486
1214 0.00493 7215
188
0.02564
1209 0.05164
7226
184
Y ∼ 1/2N (2, 1) + 1/2N (−3, 1)
CPU MWISE
CPU MWISE CPU
354
0.00314
2126
0.00186 14054
0.00306
1207 0.00172 7105
184
177
0.00316
1198 0.00187
7115
CPU MWISE
CPU MWISE CPU
357
0.00221
2345
0.00192 14153
0.00215
1209 0.00180 7121
192
184
0.00224
1199 0.00194
7144
CPU MWISE
CPU MWISE CPU
366
0.00244
2543
0.00194 14532
194
0.00234
1321 0.00184 7342
0.00246
1306 0.00194 7224
188
Y ∼ Weibull (2, 1)
CPU MWISE
CPU MWISE CPU
345
0.02082
2442
0.01742 14121
173 0.001944 1235 0.01638 7214
0.02105
1226 0.01744 7116
163
CPU MWISE
CPU MWISE CPU
356
0.02142
2448
0.01754 14098
0.02008
1212 0.01644 7032
174
186
0.02144
1176 0.01760
7034
CPU MWISE
CPU MWISE CPU
352
0.02152
2356 0.001758 14024
176
0.02089
1209 0.01654 7014
0.02154
1165 0.01760 7006
184

Table 4: Quantitative comparison between the nonrecursive Horvitz-Thompson-type KDE (2) and
two recursive
 estimators; recursive 1 correspond to the proposed estimator (1) with the choice
 of
−1
−1
(γn ) = n
, resursive 2 correspond to the estimator (1) with the choice of (γn ) = [4/5] n
. We
consider the proportion of missing at random equal to 70% and the auxiliary variable ρ, equal to
30% (in the first block), to 50% (in the second block) and equal to 70% (in the third block) of each
considered distribution, we consider three sample sizes n = 100, n = 200 and n = 500, the number
of simulations is 500, and we compute the Mean Weighted Integrated Squared Error (M W ISE)
and the CPU time in seconds.
17

0.4
0.3
0.2
0.1
0.0

Estimated Density

Nonrecursive
Recursive1
Recursive2

−1

1

2

0.10

0.15

Nonrecursive
Recursive1
Recursive2

0.8

−4

−2

0

2

4

0.2

0.4

0.6

Nonrecursive
Recursive1
Recursive2

0.0

Estimated Density

0

0.05

Estimated Density

0.20

−2

0.0

0.5

1.0

1.5

2.0

2.5

Observation

Figure 4: Quantitative comparison between the nonrecursive Horvitz-Thompson-type KDE (2) and
two recursive
 estimators; recursive 1 correspond to the proposed estimator (1) with the choice
 of
(γn ) = n−1 , resursive 2 correspond to the estimator (1) with the choice of (γn ) = [4/5] n−1 for
(a) normal, (b) normal mixture and (c) weibull distributions with 70% missing data and correlation
between the response and auxiliary variable of ρ = 0.8. Based on 500 simulated datasets of size
200.

From Figures 1 2, 3, 4, Tables 1, 2, 3 and 4, we conclude that

• The proposed recursive kernel density estimator (1), with the stepsize (γn ) = n−1 is closer to
the true density function as compared with the nonrecursive Horvitz-Thompson-type KDE (2).

• In all the considered densities, the proposed recursive KDE with the stepsize (γn ) = n−1
outperformed the nonrecursive Horvitz-Thompson-type KDE (2) in terms of estimation error
and computational costs.
• The estimators get closer to the true density function as sample size increase and the proportion of the missing at random decrease.
• The CPU time using the proposed recursive estimators are approximately two times faster than
the CPU time using the nonrecursive Horvitz-Thompson-type KDE.

18

Density




1
1
√1
√1
+
N
−1,
N
1,
2
2
2
2
1
2N

(−1, 1) + 12 N (1, 1)

N (0, 1)

Exp (1)

Cauchy (0, 1)

(x0 , n)

hCLP

CLP

hN RP I

NRPI

hRP I

RPI

(0, 100)

0.192

0.227

0.606

0.191

0.716

0.189

(0, 1000)
(0, 100)
(0, 1000)
(−1.282, 100)
(−1.282, 1000)
(0, 100)
(0, 1000)
(1.282, 100)
(1.282, 1000)
(0.1054, 100)
(0.1054, 1000)
(0.693, 100)
(0.693, 1000)
(2.303, 100)
(2.303, 1000)
(−3.078, 100)
(−3.078, 1000)
(0, 100)
(0, 1000)
(3.078, 100)
(3.078, 1000)

0.151
0.099
0.039
0.149
0.070
0.131
0.056
0.168
0.070
0.435
0.260
0.145
0.070
0.213
0.077
0.390
0.147
0.218
0.076
0.435
0.151

0.233
0.118
0.048
0.175
0.131
0.107
0.045
0.178
0.129
0.180
0.088
0.135
0.071
0.232
0.101
1.573
1.514
0.135
0.046
1.662
1.513

0.274
0.806
0.474
0.558
0.321
0.408
0.334
0.566
0.353
0.177
0.083
0.387
0.301
0.406
0.277
0.350
0.205
0.332
0.190
0.336
0.213

0.086
0.083
0.037
0.161
0.071
0.170
0.069
0.148
0.069
0.338
0.130
0.071
0.038
0.267
0.101
0.533
0.203
0.152
0.070
0.456
0.211

0.579
0.865
0.751
0.757
0.445
0.646
0.560
0.701
0.523
0.234
0.196
0.527
0.411
0.300
0.309
0.471
0.358
0.447
0.305
0.398
0.471

0.158
0.085
0.046
0.135
0.067
0.133
0.074
0.134
0.069
0.403
0.357
0.118
0.029
0.296
0.100
0.450
0.150
0.163
0.088
0.415
0.149

√
M SE(fn (x0 ))
and the bandwidth obtained respectively by CLP
Table 5: Monte Carlo estimates of
f (x0 )
method, the nonrecursive Horvitz-Thompson-type KDE (2) using the proposed plug-in method and
the proposed recursive estimator using the proposed plug-in method.
4.1.2

Local density estimation

In order to investigate the comparison between the proposed estimators in the case of local density
estimation, we followed the paper of Dutta (2014) and we consider two sample sizes: n = 100
and 1000, the number of simulations is N = 200. Moreover,
functions
√ we
 consider five√densities

f : 1- the normal mixture distribution : X ∼ 1/2N −1, 1/ 2 + 1/2N 1, 1/ 2 , 2- the normal
mixture distribution : X ∼ 1/2N (−1, 1) + 1/2N (1, 1), 3- the standard normal : X ∼ N (0, 1), 4the exponetial distribution with rate 1: X ∼ Exp (1), 5- the Cauchy distribution with location 0
and scale 1. We report in Table 5, in the first two columns the bandwidth and the Monte Carlo
estimates of the ratio of the square root of the MSE of the estimator proposed by Chan et al. (2010)
(see Table 1 of Dutta (2014)), then, we give the bandwidth and the Monte Carlo estimates of the
ratio of the square root of the MSE of the nonrecursive Horvitz-Thompson-type KDE using the
proposed plug-in bandwidth selection (31) and in the last two columns we give the Monte Carlo
estimates of the square ratio of the root of the MSE of the proposed recursive estimator using
the proposed √
plug-in
bandwidth√selection
(18). In Table 5, we consider the estimation of (a) the


1/2N −1, 1/ 2 + 1/2N 1, 1/ 2 density at x0 = 0, (b) the 1/2N (−1, 1) + 1/2N (1, 1) density
at x0 = 0, and the estimation of (c) N (0, 1), (d) Exp (1), and (e) Cauchy(0, 1) densities at x0 equal
to the 10th, 50th, and 90th percentiles.
19

From Table 5, we have the following observations.
(i) The proposed recursive estimator and the Horvitz-Thompson-type KDE, preformed better than
the CLP method in many situations: using the mixed normal densities ((a) and (b)), using
the standard normal density (c) and standard Cauchy (e) in the tail regions (10th and 90th
percentiles), using the exponential density (d) with x0 equal to the median.
(ii) The proposed recursive estimator and the Horvitz-Thompson-type KDE, seems to be consistent
in all the cases (its MISE seems to decrease as n increased), however the CLP estimator does
not seem to be consistent (see for example the mixed normal density (a)).
(ii) No one between the proposed recursive estimator and the Horvitz-Thompson-type KDE can
be claimed to be the best in all the considered cases.

4.2

Examples

In this section we report our analysis of the two datasets mentioned in the Introduction, the first
one concerne the simulated data of the Aquitaine cohort of HIV-1 infected patients and the second
one concerne the coriell cell lines using the chromosome number 11.
Simulated data : Aquitaine cohort of HIV-1 infected patients The aim of this part was to
compare the performance of the proposed recursive KDE to the nonrecursive Horvitz-Thompsontype KDE using the Aquitaine cohort of HIV-1 infected patients, see the paper Thiébault et al.
(2000) and receiving highly active anti-retroviral therapy.
We generated 200 samples of 100 subjects. The number of repeated measures for each subject
was randomly distributed between 2 and 7 (mean 4) and the times of measurements were uniformly
distributed between 0 and 6.
For comparisons, we simulated data according to the following linear mixed effects model:
Yi,k = A1 + A2 ∗ Xi,k + bi + εi,k


with the assumption that bi ∼ N 0, σ12 and εi,k ∼ N 0, σ 2 . We assumed that random coefficient
(bi , εi,k ) were independent of each other. The values of parameters were chosen to be similar to
those obtained from the Aquitaine cohort of HIV-1 infected patients. The parameters were A1 = 4,
A2 = −0.5, σ12 = 0.25 and σ 2 = 1.
Then, we fixed the proportion of missing at random to be equal to 30%.

In Figure 5, we qualitatively compared the proposed recursive KDE using the two discussed stepsizes to the nonrecursive Horvitz-Thompson-type KDE. We observed that even when 30% of the
original measurements are missing, the proposed recursive estimators remain very accurate thus
demontrating the effectiveness of our approach.
Moreover, our proposed recursive KDE was over two times faster than the nonrecursive HorvitzThompson-type KDE algorithm (2).
Our second application concerned real dataset, for which we compared the proposed recursive KDE
approach to the nonrecursive Horvitz-Thompson-type KDE algorithm using a special data using
essentially in the context of change point problem.

20

0.4
0.2
0.0

0.1

Density

0.3

Nonrecursive
Recursive1
Recursive2

2

4

6

8

HIV

Figure 5: Quantitative comparison between the nonrecursive Horvitz-Thompson-type KDE (2) and
two recursive estimators; recursive 1 correspond to the proposed estimator (1) with the choice of

(γn ) = n−1 , resursive 2 correspond to the estimator (1) with the choice of (γn ) = [4/5] n−1 .
Here we consider the fibroblast cell strains using the chromosome number 11 with the proportion of
missing at random equal to 30%. Running time was roughly 483s using the Horvitz-Thompson-type
KDE, 248s using recursive 1 and 253s using recursive 2

Real dataset: Coriell cell lines using the chromosome number 11 The data correspond to
two array Comparative Genomic Hybridization (CGH) studies of Coriell cell lines, which appears in
the DNAcopy package and bcp package (see the paper Erdman and Emerson (2007)), and for more
details about the data we advise the reader to see the paper of Snijders et al. (2001).
The data correspond to two array CGH studies of fibroblast cell strains. In particular in Venkatraman and Olshe
(2007), they chose the studies GM05296 and GM13330. After selecting only the mapped data
from chromosomes 1-22 and X, there are 2271 data points. Here we perform an analysis on the
GM13330 array CGH study described above. We simulated three different proportions of missing
at random.
It is interesting to note from Figure 6 that the three estimators are quite similar, and closely
follow the shape of the histogram. Moreover,
 we can claim that the proposed recursive estimators
−1
with the choice of the stepsize (γn ) = n
is closer to the shape of the histogram of the true data
as compared with the nonrecursive Horvitz-Thompson-type KDE (2).

21

4
3
2
0

1

Density

Nonrecursive
Recursive1
Recursive2

−0.8

−0.6

−0.4

−0.2

0.0

0.2

−0.2

0.0

0.2

−0.2

0.0

0.2

chrom11

2

3

Nonrecursive
Recursive1
Recursive2

0

1

Density

4

Histogram of chrom11

−0.8

−0.6

−0.4
chrom11

2

3

Nonrecursive
Recursive1
Recursive2

0

1

Density

4

Histogram of chrom11

−0.8

−0.6

−0.4
chrom11

Figure 6: Quantitative comparison between the nonrecursive Horvitz-Thompson-type KDE (2) and
two recursive estimators; recursive 1 correspond to the proposed estimator (1) with the choice of

(γn ) = n−1 , resursive 2 correspond to the estimator (1) with the choice of (γn ) = [4/5] n−1 .
Here we consider the fibroblast cell strains using the chromosome number 11 with the proportion of
missing at random equal respectively to 0% (in the first block), equal to 15% (in the second block),
and equal to 30% (in the last block). Running time was roughly 211s using the Horvitz-Thompsontype KDE, 107s using recursive 1 and 113s using recursive 2, when the proportion of missing at
random equal to 30%

22

5

Conclusion

This paper propose an automatic bandwidth selection of the recursive KDE under missing data (1).
The proposed estimators are asymptotically follows normal distribution. The proposed estimators
are compared to the nonrecursive Horvitz-Thompson-type KDE algorithm under missing data (2).
We showed that, using a recursive version of the Nadaraya-Watson estimator to estimate the propensity scores and using a specific data-driven bandwidth and some particularly stepsizes, the proposed
recursive estimators can give better results compared to the nonrecursive Horvitz-Thompson-type
estimator in terms of estimation error in the case of global estimation and quite similar results in
the case of local estimation. Noting that the two estimators ((1) and (2)) can deal effectively with
both missing and complete data. The simulation study confirms the nice features of our proposed
recursive estimators and satisfactory improvement in the CPU time in comparison to the nonrecursive
Horvitz-Thompson-type KDE estimator.
In conclusion, the proposed recursive estimators allowed us to obtain quite better results compared to the nonrecursive Horvitz-Thompson-type KDE under missing data in terms of estimation
error and much better in terms of computational costs.

A

Technical proofs for the recursive KDE

Throughout this section we use the following notation:
Qn =

n
Y

j=1

Zn (x) =

(1 − γj ) ,

−1
h−1
n δn πn K



x − Xn
hn



.

(19)

Let us first state the following technical lemma.
Lemma 1. Let (vn ) ∈ GS (v ∗ ), (ηn ) ∈ GS (−η), and m > 0 such that m − v ∗ ξ > 0 where ξ is defined
in (4). We have
lim vn Qm
n

n→+∞

n
X
k=1

−1
∗ −1
Q−m
.
k γk vk = (m − v ξ)

Moreover, for all positive sequence (αn ) such that limn→+∞ αn = 0, and all C ∈ R,
" n
#
X
−m
−1
m
lim vn Qn
Qk ηk vk αk + C = 0.
n→+∞

k=1

Lemma 1 is widely applied throughout the proofs. Let us underline that it is its application,
which requires Assumption (A2)(iii) on the limit of (nγn ) as n goes to infinity.
Our proofs are organized as follows. Propositions 1 and 2 in Sections A.1 and A.2 respectively,
Theorem 1 in Section A.3.

23

A.1

Proof of Proposition 1

In view of (1) and (19), we have
fn (x) − f (x) = (1 − γn ) (fn−1 (x) − f (x)) + γn (Zn (x) − f (x))


n
n−1
Y
X

(1 − γj ) γk (Zk (x) − f (x)) + γn (Zn (x) − f (x))
=
k=1



+

j=k+1

n
Y

j=1

= Qn

n
X
k=1



(1 − γj ) (f0 (x) − f (x))

Q−1
k γk (Zk (x) − f (x)) + Qn (f0 (x) − f (x)) .

It follows that
E (fn (x)) − f (x) = Qn

n
X
k=1

Q−1
k γk (E (Zk (x)) − f (x)) + Qn (f0 (x) − f (x)) .

Moreover, we have
E




Zkp (x)

=
=
=
=
=
=
=






x − Xk
δk K
hk



−p −p
p x − Xk
hk πk E 1{Tk =Xk } K
hk




−p −p
p x − Xk
hk πk P [Tk = Xk ] E K
| {Tk = Xk }
hk



−p −p
p x − Tk
hk πk P [δk = 1|Ti ] E K
hk



x
−
T
k
−p+1
E Kp
h−p
k πk
hk


Z
−p −p+1
p x−t
K
hk πk
f (t) dt
hk
R
Z
K p (z) f (x − zhk ) dz
h−p+1
πk−p+1
k

−p
h−p
k πk E

p

R

(20)

Then, it follows from (20), for p = 1, that
E [Zk (x)] − f (x) =
=

Z

R
h2k

2

K (z) [f (x − zhk ) − f (x)] dz
f (2) (x) µ2 (K) + ηk (x) ,

with

1 2 2 (2)
K (z) f (x − zhk ) − f (x) − z hk f (x) dz,
ηk (x) =
2
R
Z



24

(21)

and, since f is bounded and continuous at x, we have limk→∞ ηk (x) = 0. In the case a ≤ α/5, we
have limn→∞ (nγn ) > 2a; the application of Lemma 1 then gives
Z
n
X
1 (2)
2
Q−1
f (x) z 2 K (z) dzQn
E [fn (x)] − f (x) =
k γk hk [1 + o (1)]
2
R
k=1

+Qn (f0 (x) − f (x))


1
=
f (2) (x) µ2 (K) h2n + o (1) ,
2 (1 − 2aξ)
p

γn h−1
and (5) follows. In the case a > α/5, we have h2n = o
, and limn→∞ (nγn ) > (α − a) /2,
n
then Lemma 1 ensures that

q
n
X
−1
−1
γk hk
+ O (Qn )
E [fn (x)] − f (x) = Qn
Qk γk o
k=1

= o
which gives (6). Further, we have
V ar [fn (x)] = Q2n

n
X

q

γn h−1
n



.

2
Q−2
k γk V ar [Zk (x)]

k=1

= Q2n

n
X
k=1




2
2
2
.
E
Z
(x)
−
(E
(Z
(x)))
Q−2
γ
k
k
k
k

Moreover, in view of (20), for p = 2, that
Z

−1 −1
2
f (x − zhk ) K 2 (z) dz
E Zk (x) = hk πk
R Z
−1 −1
= hk πk f (x) K 2 (z) dz + νk (x) ,
R

with

νk (x) =

−1
h−1
k πk

Moreover, it follows from (21), that

with

Z

R

νek (x) =

R

(24)

K (z) [f (x − zhk ) − f (x)] dz.

Then, it follows from (22), (23) and (24), that
V ar [fn (x)] = f (x) R (K) Q2n

n
X

2
2 −1 −1
Q−2
k γk hk πk + Qn

−Q2n

n
X
k=1

n
X
k=1

2
2
Q−2
k γk − 2f (x) Qn

2 2
ek (x) .
Q−2
k γk ν

25

n
X

2
Q−2
k γk νk (x)

k=1

k=1

−f 2 (x) Q2n

(23)

K 2 (z) [f (x − zhk ) − f (x)] dz.

E [Zk (x)] = f (x) + νek (x) ,
Z

(22)

n
X
k=1

2
ek (x)
Q−2
k γk ν

Since f is bounded continuous, we have limk→∞ νk (x) = 0 and limk→∞ νek (x) = 0. In the case
a ≥ α/5, we have limn→∞ (nγn ) > (α − a) /2, and the application of Lemma 1 gives

−1
−1
V ar [fn (x)] = γn h−1
f (x) R (K) + o γn h−1
,
n πn (2 − (α − a) ξ)
n

4
which proves (7). Now, in the case a < α/5, we have γn h−1
n = o hn , and limn→∞ (nγn ) > 2a, then
the application of Lemma 1 gives
V ar [fn (x)] = Q2n
= o
which proves (2).

A.2

n
X

4
Q−2
k γk o hk

k=1

h4n ,



Proof of Proposition 2

Following similar steps as the proof of the Proposition 2 of Mokkadem et al. (2009), we proof the
Proposition 2.

A.3

Proof of Theorem 1

Let us at first assume that, if a ≥ α/5, then
q


D
γn−1 hn πn (fn (x) − E [fn (x)]) → N 0, (2 − (α − a) ξ)−1 f (x) R (K) .

(25)

In the case when a > α/5, Part 1 of Theorem 1 follows from the combination of (6) and (25). In
the case when a = α/5, Parts 1 and 2 of Theorem 1 follow from the combination of (5) and (25).
In the case a < α/5, (2) implies that
P

h−2
n (fn (x) − E (fn (x))) → 0,
and the application of (5) gives Part 2 of Theorem 1.
We now prove (25). In view of (1), we have
fn (x) − E [fn (x)] = Qn
Set

n
X
k=1

Q−1
k γk (Zk (x) − E [Zk (x)]) .

Yk (x) = Π−1
k γk (Zk (x) − E [Zk (x)]) .
The application of Lemma 1 ensures that
vn2 =

n
X

V ar (Yk (x))

k=1

=

n
X

2
Π−2
k γk V ar (Zk (x))

k=1

=

n
X

2 −1 −1
Π−2
k γk hk πk [f (x) R (K) + o (1)]

k=1

h
i
−1
−1 −1
(2
−
(α
−
a)
ξ)
f
(x)
R
(K)
+
o
(1)
.
= Q−2
γ
h
π
n n n
n
26

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

i

h
,
E |Zk (x)|2+p = O h−1−p
k
and, since limn→∞ (nγn ) > (α − a) /2, there exists p > 0 such that limn→∞ (nγn ) >
Applying Lemma 1, we get
!
n
n
i
i
h
h
X
X
−2−p 2+p
2+p
2+p
= O
Qk
γk E |Yk (x)|
E |Yk (x)|
k=1

k=1

= O

n
X
Π−2−p γ 2+p
k

k=1

= O

k

h1+p
k
!

γn1+p

1+p
Q2+p
n hn

1+p
2+p

(α − a).

!

,

and we thus obtain
n
1 X

vn2+p

k=1

i

h
p/2 
E |Yk (x)|2+p = O γn h−1
= o (1) .
n

The convergence in (25) then follows from the application of Lyapounov’s Theorem.

B
B.1

Results for Horvitz-Thompson-type KDE
Global density estimation using HT KDE

Now, let us recall that the bias and variance of the nonrecursive Horvitz-Thompson-type KDE
defined by (2) are given by

and

h

i

h2
e
E fn (x) − f (x) = n f (2) (x) µ2 (K) + o h2n ,
2
h

It follows that,

i
V ar fen (x) =
h i
AM W ISE fen =

πn−1
f (x) R (K) + o
nhn



1
nhn



.

h4
πn−1
I1 R (K) + n I2 h4n µ22 (K) .
nhn
4

Then, to minimize the AM W ISE of fen , the bandwidth (hn ) must equal to
!

 1/5 
I1
R (K) 1/5 −1/5 −1/5
,
πn n
I2
µ22 (K)
and we have
h i 5
4/5 1/5
AM W ISE fen = I1 I2 Θ (K) πn−4/5 n−4/5 .
4
27

(26)

To estimate the optimal bandwidth (30), we must estimate I1 and I2 . For this purpose, we use the
following two kernel estimators :
Ie1 =
Ie2 =



n
X
Xi − Xj
1
−1
K
δj π
bN
,
Wj b
n (n − 1) bn
bn

(27)

i,j=1
i6=j





n
X
1
(2) Xi − Xj
(2) Xi − Xk
−1
−1
bN W k Kb′
δj δk π
bN W j π
Kb′
.
n3 b′6
b′n
b′n
n

(28)

i,j,k=1
j6=k

where Kb and Kb′ are a kernels, bn and b′n are respectively the associated bandwidth given in (14).
We showed that in order to minimize the AM ISE of Ie1 respectively of Ie2 , the pilot bandwidth (bn )
respectively (bn ) must belong to GS (−2/5), respectively to GS (−3/14).
Then the plug-in estimator of the bandwidth (hn ) using the nonrecursive estimator (2), is given by


!1/5 
1/5
e
R (K)
−1/5
 I1
π
bN W n n−1/5  ,
(29)
2 (K)
e
µ
I2
2

and the plug-in of the AM W ISE of the nonrecursive estimator (2), is given by

B.2

h i 5
4/5 1/5
−4/5
bN W n n−4/5 .
AM^
W ISE fen = Ie1 Ie2 Θ (K) π
4

Local density estimation using the HT KDE

Now, using the nonrecursive Horvitz-Thompson-type KDE and in order to minimize the AM SE of
fen , the bandwidth (hn ) must equal to


!1/5 
1/5
f (x)
R (K)

πn−1/5 n−1/5  .
(30)
2
µ22 (K)
f (2) (x)
Moreover, since the f (x) and f (2) (x) are not known, we estimate f (x) and f (2) (x) by
fen (x) =

fen(2) (x) =



n
x − Xi
1 X
−1
δi πN W i Kb
,
nbn
bn
i=1


n
1 X
(2) x − Xi
−1
δi πN W i Kb′
.
nb′3
b′n
n
i=1

where Kb and Kb′ are a kernels, bn and b′n are respectively the associated bandwidth given in (14).
Then the plug-in estimator of the bandwidth (hn ) using the nonrecursive estimator (2), to estimate
locally the density f at a point x is given by


1/5


 fen (x) 
R (K) 1/5 −1/5 −1/5 
 
.
π
bN W n n
(31)


2 (K)
2


(2)
µ
2
fen (x)
28

Acknowledgements
We are grateful to two reviewers and the editor for their helpful comments, which have led to this
substantially improved version of the paper.

References
Bojanic, R. and Seneta, E. (1973). A unified theory of regularly varying sequences. Math. Z., 134,
91–106.
Chan, N. H. Lee, T. C. M. and Peng, L. (2010). On nonparametric local inference for density
estimation. Comput. Statist. Data Anal., 54, 509–515.
Delaigle, A. and Gijbels, I. (2004). Practical bandwidth selection in deconvolution kernel density
estimation. Comput. Statist. Data Anal., 45, 249–267.
Dempster, A. P. Laird, N. M. and Rubin, D. B. (1977). Maximum likelihood from incomplete data
via the EM algorithm. J. R. Statist. Soc. Ser. B Statist. Methodol., 39, 1–22.
Dubnicka, S. R. (2009). Kernel density estimation with missing data and auxiliary variables. Aust.
N. Z. J. Stat., 25, 1175–1179.
Duin, R. P. (1976). On the choice of smoothing parameters for parzen estimators of probablility
density functions. IEEE Trans. Comput., 25, 1175–1179.
Dutta, S. (2014). Local smoothing using the bootstrap. Comm. Statist. Simulation Comput.,
43, 378–389.
Erdman, C. and Emerson, J. W. (2007). An R package for performing a bayesian analysis of change
point problem. J. Statist. Software, 23, 1–13.
Fan, J., Heckmann, N. E. and Wand, M.P. (1995). Local polynomial kernel regression for generalized
linear models and quasi likelihood functions. J. Amer. Statist. Assoc., 90, 141–150.
Galambos, J. and Seneta, E. (1973). Regularly varying sequences. Proc. Amer. Math. Soc., 41,
110–116.
Hazelton, M. (1999). An optimal local bandwidth selector for kernel density estimation. J. Statist.
Plann. Inference, 77, 37–50.
Horvitz, D. G. and Thompson, D. J. (1952). A generalization of sampling without replacement
from a finite universes. J. Amer. Statist. Assoc., 47, 663–685.
Little, R. J. A (2002). Statistical analysis with missing data. 2nd. edn. Hoboken, John Wiley &
Sons.
Marron, J. S. (1988). Automatic smoothing parameter selction: a survey. empec., 13, 187–208.
Mokkadem, A. and Pelletier, M. (2007). A companion for the Kiefer-Wolfowitz-Blum stochastic
approximation algorithm. Ann. Statist., 35, 1749–1772.

29

Mokkadem, A. Pelletier, M. and Slaoui, Y. (2009). The stochastic approximation method for the
estimation of a multivariate probability density. J. Statist. Plann. Inference, 139, 2459–2478.
Nadaraya, E. A. (1964). On estimating regression. Theory Probab. Appl., 10, 186–190.
Parzen, E. (1962). On estimation of a probability density and mode. Ann. Math. Statist., 33, 1065–
1076.
Qi, L., Wang, C. Y. and Prentice, R. L. (2005). Weighted estimators for proportional hazards
regression with missing covariates. J. Amer. Statist. Assoc., 100, 1250–1263.
Révész, P. (1973). Robbins-Monro procedure in a Hilbert space and its application in the theory of
learning processes I. Studia Sci. Math. Hung., 8, 391–398.
Révész, P. (1977). How to apply the method of stochastic approximation in the non-parametric
estimation of a regression function. Math. Operationsforsch. Statist., Ser. Statistics., 8, 119–126.
Rosenblatt, M. (1956). Remarks on some nonparametric estimates of a density function. Ann.
Math. Statist., 27, 832–837.
Rosenbaum, P.R. and Rubin, D.B. (1983). The central role of the propensity score in observational
studies for causal effects. Biometrika, 70, 41–55.
Rudemo, M. (1982). Empirical choice of histograms and kernel density estimators. Scand. J. Stat.,
9, 65–78.
Scott, D.W. and Terrell, G.R. (1987). Biased and unbiased cross-validation in density estimation.
J. Amer. Math. Soc., 82, 1131–1146.
Sheather, S. J. (1983). A data based algorithm for choosing the window when estimating the density
at a point. Comput. Statist. Data Anal., 1, 229–239.
Sheather, S. J. (1986). An improved data based algorithm for choosing the window when estimating
the density at a point. Comput. Statist. Data Anal., 4, 61–65.
Silverman, B. W. (1986). Density estimation for statistics and data analysis. Chapman and Hall,
London.
Snijders AM, Nowak N, Segraves R, Blackwood S, Brown N, Conroy J, Hamilton G, Hindle AK,
Huey B, Kimura K, Law S, Myambo K, Palmer J, Ylstra B, Yue JP, Gray JW, Jain AN, Pinkel
D, , Albertson DG (2001). Assembly of Microarrays for Genome-wide Measurement of DNA Copy
Number. Nature Genetics, 29, 263–246.
Slaoui, Y. (2013). Large and moderate principles for recursive kernel density estimators defined by
stochastic approximation method. Serdica Math. J., 39, 53–82.
Slaoui, Y. (2014a). Bandwidth selection for recursive kernel density estimators defined by stochastic
approximation method. J. Probab. Stat, 2014, ID 739640, doi:10.1155/2014/739640.
Slaoui, Y. (2014b). The stochastic approximation method for the estimation of a distribution
function. Math. Methods Statist., 23, 306–325.
Slaoui, Y. (2015). Plug-In Bandwidth selector for recursive kernel regression estimators defined by
stochastic approximation method. Stat. Neerl., 69, 483–509.
30

Slaoui, Y. (2016). Optimal bandwidth selection for semi-recursive kernel regression estimators. Stat.
Interface, 9, 375–388.
Thiébault, R., Morlat, P., Jacqmin-Gadda, H., Neau, D., Mercié, P., Dabis, F., and Chêne, G. for
the GECSA (2000). Clinical progression of HIV-1 infected patients according to the viral response
during the first year of anti-retroviral treatment. AIDS, 14: 971–978.
Thombs, L. A., Sheather, S. J. (1992). Local bandwidth selection for density estimation. In: Proceedings of the 22nd Symposium on the Interface. New York: Springer, 111–116.
Tibshirani, R. and Hastie, T. (1987). Local likelihood estimation. J. Amer. Statist. Assoc., 82,
559–567.
Tsybakov, A. B. (1990). Recurrent estimation of the mode of a multidimensional distribution. Probl.
Inf. Transm., 8, 119–126.
Venkatraman, E.S. and Olshen, A.B. (2007). A faster circular binary segmentation algorithm for
the analysis of array CGH data. Bioinformatics, 23, 657–663.
Watson, G. S. (1964). Smooth regression analysis. Sankhya A., 26, 359–372.

31

