Fabienne Comte

COMPARISON OF SOME RECENT METHODS IN ADAPTIVE DENSITY ESTIMATION FOR BIASED DATA

Keywords: Adaptive density estimation, Biased data, Bandwidth selection, Fluorescence lifetimes, January 9, 2015

may come

Introduction

In various application settings, functional estimation can be difficult because the observed data are not a sample from the distribution of interest: this may be due to noise, missing data, censored or truncated observations. In this paper biased data models are considered where the observed distribution is the result of a (known) nonlinear distortion of the distribution of interest.

More precisely, we observe a sample Z 1 , . . . , Z n of independent identically distributed (i.i.d.) random variables with probability density function (pdf) g and cumulative distribution function (cdf) G. The observed distribution G is related to the distribution of interest, say F , by some known link function H by the following relation

G(z) = H • F (z) , z ∈ R . (1)
The aim is to recover the pdf f of the distribution of interest F in a nonparametric context using an i.i.d. sample Z 1 , . . . , Z n with distribution G and known link function H.

We have in mind the case where every Z i is the minimum of a random number N of i.i.d. random variables Y 1 , . . . , Y N with distribution F . An example in physics is the arrival time of the fasted of a random number of emitted photons [START_REF] Rebafka | A corrected likelihood approach for the pile-up model with application to fluorescence lifetime measurements using exponential mixtures[END_REF], or in biostatistics the time to the observation of a tumor originated from a clonogenic cell in the presence of a random number of competing clonogens [START_REF] Tsodikov | Estimation of survival based on proportional hazards when cure is a possibility[END_REF]. Various extensions and other examples may be considered. For example, changing the random minimum in a random maximum corresponds in actuarial science to modelling the largest claim received by an insurer in a given time interval [START_REF] Li | Preservation of stochastic orders for random minima and maxima, with applications[END_REF], or in transportation theory to the modelling of the maximal accident-free distance of a shipment of, say, explosives, with a random number of defective explosives which may explode and cause an accident during transport [START_REF] Shaked | Stochastic comparisons of random minima and maxima[END_REF].

(1) : MAP5, UMR 8145 CNRS, Sorbonne Paris Cité, Paris Descartes University, France, email: fabienne.comte@parisdescartes.fr (2) : LPMA, University of Paris 6, UPMC, France, email: tabea.rebafka@upmc.fr. The authors wish to thank PicoQuant GmbH, Berlin, Germany for kindly providing the TCSPC data.

From a methodological point of view, we propose projection and kernel estimators associated with model and bandwidth selection devices. Two properties that hold in the model given by (1) give rise to two different ways of correcting the bias in the data. Hence two different estimation strategies can be used to construct both kernel and projection estimators and they are worth being compared. The mean-square risks of the estimators are studied and oracle-type risk bounds are provided. Adaptive projection estimators correspond to methods originally described by [START_REF] Barron | Risk bounds for model selection via penalization[END_REF] and applied to survival analysis and biased data by Efromovich (2004a,b) and Brunel et al. (2005); more recently, wavelet projection estimators have been studied by [START_REF] Chesneau | Wavelet block thresholding for density estimation in the presence of bias[END_REF], [START_REF] Cutillo | Wavelet density estimation for weighted data[END_REF]. For the bandwidth selection of the kernel estimators the recent approach of [START_REF] Goldenshluger | Bandwidth selection in kernel density estimation: oracle inequalities and adaptive minimax optimality[END_REF] is applied to our model and studied from both a pointwise and a global point of view, and our results on this side, namely finite sample risk bounds for adaptive kernel estimators, are new.

It is worth mentioning that our model can be related to other biased data contexts, which have been studied from other or specific point of view by several authors: strategies for estimating cumulative distribution functions are proposed by [START_REF] Gill | Large sample theory of empirical distributions in biased sampling models[END_REF], [START_REF] Wu | Minimax kernels for density estimation with biased data[END_REF], [START_REF] Wu | A cross-validation bandwidth choice for kernel density estimates with selection biased data[END_REF], [START_REF] Efromovich | Distribution estimation for biased data[END_REF], El [START_REF] El Barmi | Transformation-based density estimation for weighted distributions[END_REF]; the specific case of length-biased sampling has been studied in a lot of papers, see [START_REF] Vardi | Nonparametric estimation in the presence of length bias[END_REF], [START_REF] Jones | Kernel density estimation for length biased data[END_REF], de Uña-Álvarez ( 2004), de Uña-Álvarez and [START_REF] De Uña-Álvarez | Comparing nonparametric estimators for length-biased data[END_REF] [START_REF] Asgharian | Length-biased sampling with right censoring: an unconditional approach[END_REF], among others.

A simulation study is performed to calibrate and compare all those methods. Several questions are in order: What is a good choice of the penalty constants? Can we adapt the existing proposals for model selection to the kernel methods? How do the pointwise and the global strategy compare in specific examples? As each estimator involves a bias correction, is there one that outperforms the others? In our simulation study we focus on the so-called pile-up model, which is used for fluorescence lifetime measurements [START_REF] O'connor | Time-correlated single photon counting[END_REF] and presented in detail in the experimental section of the paper. An application to real fluorescence data is also provided.

The paper is organized as follows. Section 2 presents the model. In Section 3 adaptive kernel and projection estimators are proposed. Section 4 gives theoretical results on risk bounds of the different estimators. In the simulation study (Section 5) different aspects of the estimators are compared. Finally, Section 6 presents the proofs for the theoretical results of the paper.

Model and Assumptions

2.1. Notations. For two functions u and v, we denote by u

• v the function x → u • v(x) := u(v(x)).
If u is a one-to-one map, we denote by u -1 the inverse of the function u, that is the function such that (u

-1 •u)(x) = (u•u -1 )(x) = x for all x.
The derivative of u is denoted by u and the second-order derivative by ü, provided that they exist. The standard convolution product is given by u * v(x) = u(t)v(xt)dt. Furthermore, we denote by . p the L p -norm given by u 1) is necessarily increasing and surjective so that H • F is a cdf for any cdf F . Note that even if H is not injective, G given by (1) may still be a cdf. However, our goal is to recover the density of distribution F using a sample from G. Hence, we need H to be bijective to ensure identifiability of the model. Indeed, when H is a one-to-one map, then

p p = |u(x)| p dx and by . ∞ the L ∞ -norm, u ∞ = sup x∈R |u(x)|.

Model and assumptions. The link function

H : [0, 1] → [0, 1] in relation (
F (z) = H -1 • G(z) , z ∈ R . (2)
Furthermore, we assume that H is differentiable, since in our context both F and G are supposed to be absolutely continuous with pdf f and g respectively. Then, deriving relation (1) implies that the densities f and g are related by

g(z) = Ḣ • F (z) f (z) , z ∈ R . (3) Combining (2) and (3) gives f = g/ Ḣ • F = g/ Ḣ • H -1 • G. Define the weight function w by w(u) = 1 Ḣ • H -1 (u) , u ∈ [0, 1] ,
then we obtain

f (z) = w • G(z) g(z) , z ∈ R . (4)
The weight function w is well defined under the assumption that Ḣ is bounded away from zero. Furthermore, we shall require that w is Lipschitz. This is ensured if there exist finite constants a, b > 0 such that

(5) Ḣ(u) ≥ a , | Ḧ(u)| ≤ b , u ∈ [0, 1] ,
and a finite constant d > 0 such that

(6) Ḣ(u) ≤ d , u ∈ [0, 1] . Indeed then 1/d ≤ w(u) ≤ 1/a and ẇ(u) = -Ḧ • H -1 (u)/[ Ḣ • H -1 (u)] 3
, so that the Lipschitz constant of w, say c w , is such that c w ≤ b/a 3 . Moreover, we may possibly require that f or g is bounded. Note that f and g are either both bounded or both unbounded, since

a f ∞ ≤ g ∞ ≤ d f ∞ .
From relation (4) follows a fundamental property that holds for any measurable bounded function ψ,

(7) E[ψ(Y )] = E [ψ(Z) w • G(Z)] ,
where Y has distribution F and Z is distributed as G.

(z) = 1 n n i=1 1{Z i ≤ z} yields a natural estimator of E[ψ(Y )] given by (8) L = 1 n n i=1 ψ(Z i ) w • Ĝn (Z i ) .
This is a useable estimator since the link function H is supposed to be known and so is the weight function w.

Denote by Z (i) the i-th order statistic associated with (Z 1 , . . . , Z n ) satisfying

Z (1) ≤ • • • ≤ Z (n) . Note that w • Ĝn (Z (i) ) = w(i/n). Then, we can rewrite L as (9) L = 1 n n i=1 ψ(Z (i) ) w i n .
We see that L takes the form of a so-called L-statistic, i.e. a linear combination of order statistics.

Adaptive density estimators

3.1. Estimation strategies. Here several strategies to estimate f are presented, namely kernel and projection estimators. All estimators make use of relation ( 4). The first strategy consists in first estimating g from the data, and then multiplying this estimate by w( Ĝn (z)) to correct the bias. This estimator is referred to as the plug-in estimator. As an estimator of g we use a kernel estimator, for which a bandwidth selection method is provided as well.

Concretely, let ĝh be the standard kernel estimate of g given by ( 10)

ĝh (x) = 1 n n i=1 K h (x -Z i ) ,
where K is a kernel, that is, an integrable function such that K(u)du = 1, h is a bandwidth parameter and

K h (u) = h -1 K (u/h) . Then plugging ĝh into relation (4) yields the estimator f (1) h (x) of f (x) defined by (11) f (1) h (x) = ĝh (x)w Ĝn (x) .
The second method is a kernel estimator as well, however we directly estimate f by using property ( 7) and an L-statistic of the form of (8). By taking ψ = K h (x -•) in ( 8), a kernel estimator of the target density f is obtained by

(12) f (2) h (x) = 1 n n i=1 w • Ĝn (Z i )K h (x -Z i ) = 1 n n i=1 w i n K h (x -Z (i) ) .
Indeed, the plug-in method as well as property (7) can also be used to construct projection estimators. Here the approach consists in approximating the orthogonal projection of g (or f ) onto some function space. More precisely, suppose that the restriction of g (resp. f ) on some interval A is square integrable, that is

g1 A ∈ L 2 (A) (resp. f 1 A ∈ L 2 (A)). Let (ϕ j ) 0≤j≤2m be the trigonometric basis (ϕ j ) 0≤j≤2m on A = [a, b] defined by ϕ j = (b - a) -1/2 ϕ 0 j ((x -a)/(b -a)) and ϕ 0 0 (x) = 1 [0,1] (x), ϕ 0 2j+1 (x) = √ 2 cos(2πjx)1 [0,1] (x) for j ≥ 0, ϕ 2j (x) =
√ 2 sin(2πjx) for j ≥ 1. Then define the subspace S m = Span(ϕ j , j = 0, 1, . . . , 2m) and D m = dim(S m ) = 2m + 1. Then, the orthogonal projection g m (resp. f m ) in the L 2 -sense of g (resp. f ) on S m is given by g m = 2m j=0 a

(1) j ϕ j with a

(1) j = g, ϕ j (resp.

f m = 2m j=0 a
(2) j ϕ j with a

(2) j = f, ϕ j ). Using the plug-in method, we obtain the projection-type estimate f (1) m defined by

f (1) m (x) = ĝm (x)w Ĝn (x) , with ĝm = 2m j=0 â(1) j ϕ j and â(1) j = 1 n n i=1 ϕ(Z i ) .
By using property (7), a second projection-type estimator of f is given by

(13) f (2) m = 2m j=0 â(2) j ϕ j , with â(2) j = 1 n n i=1 ϕ j (Z (i) )w(i/n) .
Next, we propose different data-driven bandwidth and model selection methods.

3.2. Bandwidth selection. For a data-driven choice of bandwidth h, consider a finite collection H of bandwidths given by ( 14)

H = h k , k = 1, . . . , H n , 1 n ≤ h k ≤ 1 , with H n ≤ n .
First an adaptive pointwise estimator of f (x 0 ) for some fixed x 0 is presented. As suggested in [START_REF] Goldenshluger | Bandwidth selection in kernel density estimation: oracle inequalities and adaptive minimax optimality[END_REF], define the estimators ĝh,h ′ and f (2) h,h ′ depending on two bandwidths by ĝh,h

′ (x) = K h ′ * ĝh (x) and f (2) h,h ′ (x) = K h ′ * f (2) h (x)
. Notice the symmetry of both estimators in h and h ′ . Next we denote, for i = 1, 2,

V (i) 0 (h) = κ (i) 0 C i K 2 1 K 2 2 log n nh , C 1 = g ∞ , C 2 = f ∞ /a , A (i) 0 (h, x 0 ) = sup h ′ ∈H θ(i) h,h ′ (x 0 ) - θ(i) h ′ (x 0 ) 2 -V (i) 0 (h ′ ) + , θ (1) = g , θ (2) = f (2) , (15) ĥ(i) (x 0 ) = arg min h∈H A (i) 0 (h, x 0 ) + V (i) 0 (h) .
The numerical constants κ (i) 0 are typically calibrated by simulation. The other constants are known, except g ∞ or f ∞ which in practice are replaced by some estimators (see Section 5.2). Now two adaptive pointwise estimators of f (x 0 ) for some fixed x 0 are obtained by ( 16) f (1) ĥ(1) (x 0 )

(x 0 ) = ĝĥ(1) (x 0 ) (x 0 )w( Ĝn (x 0 )) and f (2) ĥ(2) (x 0 ) (x 0 ) .
The risk bounds associated with these estimators are given in Theorem 4.1. In a similar way, a procedure for global bandwidth selection is developed. Denote for i = 1, 2,

V (i) (h) = κ (i) 1 D i max( K 2 1 , 1) K 2 2 nh , D 1 = 1 , D 2 = 1 a , (17)
A (i) (h) = sup h ′ ∈H θ(i) h,h ′ - θ(i) h ′ 2 -V (i) (h ′ ) + , with θ (1) = g , θ (2) = f (2) , (18) ĥ
(i) = arg min h∈H A (i) (h) + V (i) (h) , (19)
where κ (i) 1 are numerical constants. Then define the adaptive estimators f (1) ĥ(1) (x) = ĝĥ(1) (x)w( Ĝn (x)) and f (2) ĥ(2) , for which a risk bound is given in Theorem 4.2.

Model selection.

For the projection estimators the classical penalization approach by [START_REF] Barron | Risk bounds for model selection via penalization[END_REF] can be applied. In Section 4 we show that a bias-variance trade-off is achieved. Define penalty terms pen (i) (m) by ( 20)

pen (1) (m) = κ (1) 2 D m n , pen (2) (m) = κ (2) 2 w 2 2 D m n ,
with appropriate constants κ

m 2 + pen (2) (m) ,
and consider the density estimates f (i) m(i) , for i = 1, 2.

Theoretical results on the estimators

In this section we first provide results for the kernel estimators with super-index (2) . Then the flavor of results for kernel estimators with super-index (1) are discussed. Lastly, we study the projection estimators. All proofs are relegated to Section 6. 4.1. Study of the pointwise and integrated risk of kernel estimators. Classically, for f or g we consider Hölder classes Σ(β, C) defined as

Σ(β, C) = {f : T → T, f (ℓ) exists for ℓ = ⌊β⌋and |f (ℓ) (x)-f (ℓ) (x ′ )| ≤ C|x-x ′ | β-ℓ , ∀x, x ′ ∈ T},
where T ⊂ R. The set T may be an interval or the entire set R depending on the support of the kernel. For instance, if K has compact support [-1, 1], the set T for a given point x 0 under study and bandwidth less than 1 can be taken equal to [x 0 -1,

x 0 + 1].
The following assumptions are useful to provide risk bounds. (H1) f belongs to the Hölder class Σ(β, C).

(H2) K is a kernel of order ℓ = ⌊β⌋ satisfying |u| β |K(u)|du < ∞. Furthermore, u 2 K(u)du <
∞ and K ∞ < ∞. We recall that a kernel of order ℓ satisfies x k K(x)dx = 0 for k = 1, . . . , ℓ. Assumption (H2) is classical.

The following proposition is easily shown, and proved in Section 6.

Proposition 4.1. If f is bounded and ( H1)-( H2), ( 5) and ( 6) are fulfilled, then, for any x 0 , the estimator f (2) h defined by ( 12) satisfies

E[( f (2) h (x 0 ) -f (x 0 )) 2 ] ≤ 3C 2 1 h 2β + C 2 nh ,
where

C 1 = C |x| β |K(u)|du/ℓ! and C 2 = 3 1/a + 5db 2 /(4a 6 ) f ∞ K 2 2 .
An upper bound of the mean integrated squared error (MISE) of fh can be easily derived for densities f belonging to the Nikol'ski class. Let β > 0 and L > 0. Define the Nikol'ski class N (β, L) as

N (β, L) = f : R → R : f (ℓ) (x + t) -f (ℓ) (x) 2 dx 1/2 ≤ L|t| β-ℓ , ∀t ∈ R .
Proposition 4.2. Let β > 0 and L > 0. Assume that f is a square integrable density function belonging to the Nikol'ski class N (β, L) and that ( 2) and ( 5) hold. Then the MISE of f (2) h defined by ( 12) satisfies

E f (2) h -f 2 2 ≤ 3C 2 3 h 2β + C 4 nh ,
where

C 3 = [L/(ℓ -1)!] |u| β |K(u)|du and C 4 = 3(b 2 /a 6 + 1/a) K 2 2 .
Risk bounds for the plug-in kernel estimators are a consequence of usual density estimation results or of the above bounds in the particular case of w ≡ 1, and of the following inequality.

f (1) h (x 0 ) -f (x 0 ) 2 = w • Ĝn (x 0 )(ĝ h (x 0 ) -g(x 0 )) + (w • Ĝn (x 0 ) -w • G(x 0 ))g(x 0 ) 2 ≤ 2 a 2 [ĝ h (x 0 ) -g(x 0 )] 2 + 2b 2 a 6 g 2 (x 0 ) Ĝn (x 0 ) -G(x 0 ) 2 . (22) Since E([ Ĝn (x 0 ) -G(x 0 )] 2 ) ≤ 1/n, E [ f (1) h (x 0 ) -f (x 0 )] 2 ≤ 2 a 2 E [ĝ h (x 0 ) -g(x 0 )] 2 + 2b 2 g 2 (x 0 ) a 6 1 n .
Thus, Proposition 4.1 holds for f (1) h , for different constants C 1 , C 2 , under assumptions on g. Analogously, in the integrated case, we get

E f (1) h -f 2 2 ≤ 2 a 2 E ĝh -g 2 2 + 2 b 2 g 2 2 a 6 1 n .
It is noteworthy that here the regularity parameter involved in the bounds is related to g instead of f . Clearly, both pointwise and integrated risks admit a bias-variance decomposition with standard terms. If the bandwidth could be chosen of order n -1/(2β+1) , where β is the Hölder or the Nikol'ski regularity index, then the resulting rate of the estimators would be of order n -2β/(2β+1) . As β is unknown, this compromise cannot be performed in that naive way. For this reason, we propose data-driven methods for bandwidth selection. 4.2. Adaptive pointwise kernel estimation. Now we consider bandwidth collection satisfying the following assumption. (H3) The collection H of bandwidths is a finite set given by ( 23)

H = h k , k = 1, . . . , H n , 1 n ≤ h k ≤ 1 , with H n ≤ n ,
and there exists some finite constant S (independent of n) such that

(24) 1 n h∈H 1 h ≤ S .
The definition of the bandwidth collection H via ( 23) is very general, only condition (24) requires some comments and illustration. We give some examples satisfying (24).

(C1) For a ∈ (0, 1), the bandwidth collection H = {h k = (k/n) a , k = 1, . . . , n} satisfies (H3) with H n = n and S = 1/(1a). Note that the case a = 1 does not fulfill (H3), but S = log(n) would be admissible provided that log(n)/n is replaced by log

2 (n)/n in inequality (25) below. (C2) The collection H = {h k = 2 -k , k = 1, . . . , [log 2 (n)]} is another example satisfying (H3) with H n = [log 2 (n)] and S = 2. (C3) The collection H = {h k = 1/k, k = 1, . . . , [ √ n]} satisfies (H3) with H n = [ √ n] and S = 1. Now the following result holds for the estimator f (2) ĥ(2) (x 0 ) (x 0 ).
Theorem 4.1. Assume that ( 5), ( 6) and assumptions (H1)-( H3) hold. Then, there are constants C * , C > 0 depending only on a, b, f ∞ , S, K 1 , K 2 , K ∞ and on the Hölder class parameters β and C such that

(25) E f (2) ĥ(2) (x 0 ) (x 0 ) -f (x 0 ) 2 ≤ C * inf h∈H h 2β + V (2) 0 (h) + C log n n .
Assumption (H1) requires the regularity of the density f , however for our estimation procedure the order β and the constant C need not to be known.

The right-hand side of ( 25) is of order (n/ log(n)) -2β/(2β+1) provided that H contains bandwidths h k of order n -1/(2β+1) . This is the case for collections [C1] and [C2], and also for [C3] if β ≥ 1/2. Then (25) implies an almost optimal compromise between the two terms of the bound given in Proposition 4.1: the bias-variance trade-off is realized, with a loss of order log(n), which is classical for pointwise adaptive procedures. In density estimation (corresponding to w ≡ 1), such a loss is known to be unavoidable and thus adaptive minimax (see [START_REF] Butucea | Two adaptive rates of convergence in pointwise density estimation[END_REF]).

4.3. Adaptive global kernel estimator. For the integrated risk the following result holds.

Theorem 4.2. Assume that ( 5) and assumptions ( H2)-( H3) hold, and that the bandwidth collection H is such that for any c > 0, there exists a finite constant

Σ(c) (independent of n) such that (26) h∈H e -c/h ≤ Σ(c) . Denote f h = K h * f . Then there exists a constant κ (2) 1 such that E f (2) ĥ(2) -f 2 2 ≤ C inf h∈H f -f h 2 2 + V (h) + C ′ log n n , where C is a numerical constant and C ′ is a constant depending on a, b, f ∞ , S, K 1 , K 2 , K ∞ .
It is worth emphasizing that this inequality proves that a bias-variance trade-off is achieved in a nonasymptotic way and without any regularity assumption on f . Moreover, there is no additional log(n) factor in the definition of V (h), contrary to V 0 (h). As a consequence, if f ∈ N (β, L), where β and L need not to be known, then ff h 2 2 ≤ Ch 2β , and thus inf

h∈H f -f h 2 2 + V (h) = O n -2β/(2β+1) ,
provided that H contains bandwidths of order n -1/(2β+1) , which is the case in example (C2) as well as in (C3) if β ≥ 1/2. Therefore, the best bias-variance trade-off is automatically achieved by the procedure.

4.4. Plug-in kernel estimators. The results for the plug-in kernel estimators are a consequence of the previous ones. It follows from ( 22) and the bound

E([ Ĝn (x 0 )-G(x 0 )] 2 ) ≤ 1/n that E ( f (1) ĥ(1) (x 0 ) (x 0 ) -f (x 0 )) 2 ≤ 2 a 2 E (ĝ ĥ(1) (x 0 ) (x 0 ) -g(x 0 )) 2 + 2b 2 g 2 (x 0 ) a 6 1 n .
Clearly, inequality (25) holds for w ≡ 1. Furthermore, if g satisfies the assumptions of Theorem 4.1 and belongs to a Hölder class with regularity parameter β * , the following bound holds

(27) E f (1) ĥ(1) (x 0 ) (x 0 ) -f (x 0 ) 2 ≤ C * inf h∈H h 2β * + V (1) 0 (h) + C′ log n n .
Therefore, the risk bound on f (1) ĥ(1) (x 0 ) (x 0 ) is an automatic compromise related to the regularity of g and is optimal if f and g belong to the same Hölder space.

Analogously, in the integrated case, we get

E f (1) ĥ(1) -f 2 2 ≤ 2 a 2 E ĝĥ(1) -g 2 2 + 2 b 2 g 2 2 a 6 1 n .
The result of Theorem 4.2 can be obtained for f (1) ĥ( 1) , but with squared bias term related to g instead of f . Thus, it is optimal if f and g belong to the same Nikol'ski space. 4.5. Projection estimators. The following risk bound holds for the projection estimator.

Proposition 4.3. Consider the estimator (13), then

(28) E f (2) m -f 1 A 2 2 ≤ f 1 A -f m 2 2 + C D m n ,
where

C = 2(1/a + b 2 /a 6 ).
To evaluate the bias of a projection estimator, it is common to consider regularity spaces that are different from those used in kernel estimation. Let Barron et al., 1999, Lemma 12). Thus, choosing D m * = O(n 1/(2α+1) ) in inequality (28) yields that the mean square risk satisfies 2α+1) ). This rate is known to be optimal in the minimax sense for density estimation for direct observations [START_REF] Donoho | Density estimation by wavelet thresholding[END_REF].

f 1 A = f A belong to a ball of some Besov space B α,2,∞ (A) with r + 1 ≥ α. Then for f A α,2,∞ ≤ L we have f A -f m 2 2 ≤ C(α, L)D -2α m (
E( fm * -f A 2 2 ) ≤ O(n -2α/(
Using this approach, the following result can be shown.

Theorem 4.3. Assume that m n ≤ O( √ n) and that f is bounded on A, i.e. f ∞ < ∞.
Then there exists a numerical constant κ

(2)

2 such that we have

(29) E f - f (2) m(2) 2 2 ≤ C inf m∈Mn f -f m 2 2 + w 2 2 D m n + K log 2 (n) n ,
where C is a numerical constant and K depends on a, b, f ∞ and the basis.

Risk bounds of the form (29) are often called oracle inequality. Note that the last term c log 2 (n)/n is clearly negligible with respect to the order of the infimum (in particular, in all Besov cases described above). This result can easily be generalized to other bases, such as piecewise polynomials or wavelets.

The proof of the theorem relies on Talagrand's inequality and follows the line of the proof of Theorem 4.2 in Brunel and Comte (2005). Therefore, only a sketch of the proof is provided in Section 6.

Lastly, using inequality (22) yields that f (1) m(1) leads to an optimal bias-variance trade-off with respect to density g, and is optimal if f and g belong to the same Besov space.

Experimental study

We have at hand six estimators with (nearly-)optimal rates corresponding to different statistical methods that are intrinsically interesting to compare. To this end this section provides a simulation study and numerical results on a real data example in the specific case of the pile-up model. 5.1. Pile-up model. In the pile-up model the minimum of a random number of variables is observed. More precisely, let {Y k , k ≥ 1} be a sequence of i.i.d. random variables with pdf f and cdf F . Let N be a random variable taking its values in

N * = {1, 2, . . . } independently of this sequence. Let Z = min{Y 1 , . . . , Y N }. The cdf of Z denoted by G is related to F by G(z) = H•F (z) , z ∈ R , with H(u) = 1-M (1-u) and M (u) = E u N , u ∈ [0, 1] .
The estimation problem consists in extracting density f from an i.i.d. sample Z 1 , . . . Z n from the pile-up distribution G, when the distribution of N is known and of Poisson type. This problem arises with data in time-resolved fluorescence, where only the shortest arrival time of a group of photons can be observed [START_REF] O'connor | Time-correlated single photon counting[END_REF]. Here it is assumed that the size N of a group of photons follows the Poisson distribution restricted to N * with known parameter µ. The renormalized probability masses are given by P(N = k) = µ k /k!/(e µ -1). As µ is known, the link function H is known as well with M (u) = (e µu -1)/(e µ -1). We see that the pile-up model is a special case of model ( 1) and model assumptions ( 5) and ( 6) are fulfilled with a = µ/(e µ -1), d = µe µ /(e µ -1) and b = µ 2 /(1e -µ ). Furthermore, the weight function is given by w(u) = (1e -µ )/[µ(u(e -µ -1) + 1)].

Methods with plug-in approach

5.2. Computational issues and calibration. We implemented the adaptive pointwise kernel estimators f (i) ĥi (x 0 ) (x 0 ), i = 1, 2 defined by ( 16) with ĥi (x 0 ) given by ( 15), the adaptive global kernel estimators f (i) ĥ(i) , i = 1, 2 given by ( 17)-( 19) as well as the adaptive projection estimators f (i) m(i) , i = 1, 2 described in Section 3.3. In the following simulations the bandwidth collection (C2) is used for the kernel estimators. Indeed, collections (C1) and (C3) are much larger, without leading to proportionally better results. For the kernel estimators we used the gaussian kernel, i.e. K(u) = e In the quantities We noted that the projection estimators are much improved by normalizing f (i) m(i) , i = 1, 2 such that its integral equals one. However, normalization is only appropriate when the interval where the density is estimated covers the main support of the density. For the kernel estimators normalization does not seem to be necessary. In fact, the property is almost automatic if K is a density because then fh (x

V (i) 0 (h), i = 1, 2 the value of f ∞ resp. g ∞ is replaced by estimators. More precisely, g ∞ is approximated by the 95th percentile of {max x 0 ĝh (x 0 ), h ∈ H}. Likewise, f ∞ is approximated by the 95th percentile of {max x 0 f (2) h (x 0 ), h ∈ H}. The terms θ(i) h,h ′ - θ(i) h ′ 2 , i = 1,
)dx = (1/n) n i=1 w(i/n) ≃ 1 0 w(u)du = 1. The different penalty constants κ (i)
j are calibrated via simulation. For every estimator the mean integrated squared error ff 2 is approximated on a grid of κ-values for different distributions. The aim is to choose κ such that the MISE is minimized simultaneously for all distributions.

Figure 1 represents the results for the following five distributions for f : normal N (10, 1), Laplace L(10, 3), Beta B(5, 1), Beta B(5, 2) and exponential E(1) distribution. The Poisson parameter is set to µ = 0.8. For every point of the grid of κ-values, 1000 datasets of sample size 500 are generated, and the 6 estimators and the associated MISE are evaluated. The mean values of the MISE are represented in Figure 1. For the sake of readability, the MISE curves are rescaled such that they are contained in the interval [0, 1]. Recall that only the point matters where the MISE attains the minimum, and not the value of the minimum.

One observes that the MISE curves for the methods based on L-statistics (first row) and the plug-in method (second row) are always quite similar. Hence, the same κ may be used for both methods, i.e. κ

(1) j = κ

(2) j for j = 1, 2, 3. However, the MISE curves are quite different from one estimation strategy to another.

The global kernel estimators as well as the projection estimators are rather robust with regard to the choice of κ, as there exists an interval of κ-values where all MISE curves are quite flat and achieve there minimum. In the following we set κ • five different distributions: normal N (10, 1), Laplace L(10, 3), Beta B(5, 1), Beta B(5, 2), exponential E(1), giving a total of 30 settings.

To evaluate the performance of the estimators we proceed exactly as in the calibration study. For each setting the estimators and their MISE are evaluated on 1000 datasets. The boxplots in Figure 2 represent the corresponding results when the sample size is 500. Table 1 shows all means of the MISE in the different settings.

We now analyze the impact of the different factors on the performance of the estimators. Therefore, we study how the mean MISE evolves when one of the factors changes.

Impact of the sample size. As usual, increasing the sample size results in a decrease of the MISE. Interestingly, depending on the estimation strategy and on the type of distribution, this decrease can be more or less pronounced. The increase of the sample size is more beneficial for the projection estimators than for the kernel estimators. The improvement is by far more important for the normal distribution than for the Beta B(5, 1) or the exponential distribution.

Impact of the Poisson parameter.

In the pile-up model the Poisson parameter µ is related to the degree of distortion. In other words, it represents the amount of bias in the model. The three levels of µ considered here correspond to a very low (µ = 0.1), medium (µ = 0.8) and high (µ = 2) degree of distortion. As increasing µ results in a more difficult estimation problem, we observe increasing MISE values in almost all settings, see Figure 2. However, there are some exceptions where the MISE decreases when µ increases. We do not have any explanation for this phenomenon.

Comparison of estimation strategies: Global kernel, pointwise kernel or projection strategy?

Depending on the underlying target distribution there are significant differences in the performance of the different estimation strategies. Both projection methods are best for the normal distribution and the Beta B(5, 2)distribution. For the Beta B(5, 1)-distribution they perform similarly good as the pointwise kernel methods. However, in the exponential case there is a stark difference between the projection estimator based on L-statistics, which is best in all settings, and the projection plug-in estimator being by far the worst method for large µ.

The global kernel estimators are best in the Laplace case, and do generally well, except for the Beta distributions.

Pointwise kernel estimators are traditionally known to be well adapted to capture peaks like in the exponential or Laplace distribution. Here, we cannot observe this property, since in almost all settings global strategies (kernel and projection) are doing better.

Consequently, there is no estimation strategy that outperforms the others in all set-ups, but a preference should be given to projection strategies and global kernel estimators.

Comparison of bias correction approaches: L-statistics or plug-in strategy?

When µ = 0.1 there is almost no bias in the model. Consequently, no significant difference between the methods is observed. However, for larger µ, differences in the MISE appear for the different bias correction methods. Increasing n amplifies the difference between the methods.

Depending on the type of distribution, one or the other correction approach is preferable. In the case of the normal, Laplace and the exponential distribution, all estimators using Lstatistics always lead to better results than the corresponding plug-in estimators. However, in the case of the Beta distributions, the plug-in versions mostly give better results. 5.4. Application to real data. We now apply our statistical methods to real fluorescence lifetime measurements. The data are supposed to be generated by the pile-up model with Poisson parameter µ = 0.166. The target density f is known to be the convolution of an exponential density f E with mean 2.54 ns and some positive noise f η coming from the measuring instrument, i.e. f = f η ⊗ f E . The noise density f η can be observed separately, so that we assume that f η is a known function. Figure 3(a) gives a histogram of an independent noise sample of size 259, 386. Note that the same dataset has already been analyzed in ?, but from a deconvolution point of view, that is the aim was the recovery of the exponential density f E . Here we are interested in the estimation of f = f η ⊗ f E .

Figure 3(b) shows the data in form of a histogram with very fine bins and target density f . The sample size is 17, 402. Here we clearly observe the pile-up effect, that is the histogram is biased in the sense that mass is shifted to the origin compared to the original density f . On this dataset the pointwise and global kernel estimators coincide for both strategies (Lstatistics and plug-in). The two projection estimators differ only slightly. For illustration, Figure 3(c) displays the plug-in projection estimator and the global kernel estimator based on L-statistics compared to density f . We can see that the projection estimator gives a very good recovery of the target density f , whereas the kernel estimator seems to do too much bias correction. Indeed, the corresponding squared errors f -f 2 are given by kernel estimator using L-statistics (pointwise and global) 6.056 10 -3 kernel estimator by plug-in (pointwise and global) 6.214 10 -3 projection estimator using L-statistics 0.433 10 -3 projection estimator by plug-in 0.431 10 -3 . Consequently the plug-in projection estimator gives the best approximation.

6. Appendix 6.1. Proof of Proposition 4.1. Let x 0 be a fixed point. Denote by fh the pseudoestimator of f given by fh

(x) = 1 n n i=1 w • G(Z i )K h (x -Z i ). We write (30) f (2) h (x 0 ) -f (x 0 ) = f (2) h -fh (x 0 ) + fh -E[ fh ] (x 0 ) + E[ fh ] -f (x 0 ) .
First, we state that by property (7) we have

(31) E[ fh (x 0 )] = E[K h (x 0 -Y i )] = K h * f (x 0 ) =: f h .
The last term in (30) is a standard bias term in kernel density estimation [START_REF] Tsybakov | Introduction à l'estimation non-paramétrique[END_REF]

). Denote b(x) = E[ fh ] -f (x). As K(u)du = 1, b(x 0 ) = K h (x 0 -y)f (y)dy -f (x 0 ) = K(u) [f (x 0 -uh) -f (x 0 )] du .
By a standard Taylor expansion of f , we get (see e.g. Tsybakov ( 2004))

(32) |b(x 0 )| ≤ C |x| β |K(x)|dx ℓ! h β = C 1 h β .
To study the second term of (30), we successively apply property (7), 0 ≤ w ≤ 1/a and the fact that K is square-integrable to obtain

E fh -E( fh ) 2 (x 0 ) = 1 n Var (w • G(Z 1 )K h (x 0 -Z 1 )) ≤ 1 n E {w • G(Z 1 )K h (x 0 -Z 1 )} 2 ≤ 1 an E K h (x 0 -Y 1 ) 2 = 1 anh 2 K 2 x 0 -y h f (y)dy ≤ 1 anh f ∞ K 2 2 . ( 33
)
For the first term in decomposition (30), the Lipschitz property of w implies that

E f (2) h -fh 2 (x 0 ) ≤ c 2 w n n i=1 E ( Ĝn -G) 2 (Z i )K 2 h (x 0 -Z i ) = c 2 w n n i=1 E Ĝn,i - n -1 n G 2 (Z i )K 2 h (x 0 -Z i ) + E 1 n (1 -G(Z i )) 2 K 2 h (x 0 -Z i ) ,
since the cross product term is centered, where Ĝn,i (x

) = n -1 n j=1,j =i 1 Z j ≤x . Then E Ĝn,i - n -1 n G 2 (Z i )K 2 h (x 0 -Z i ) = E E Ĝn,i - n -1 n G 2 (Z i )K 2 h (x 0 -Z i ) Z i = E n -1 n 2 G(Z i )(1 -G(Z i ))K 2 h (x 0 -Z i ) ≤ 1 4n E K 2 h (x 0 -Z i ) ≤ g ∞ K 2 2 4nh .
Therefore, as g ∞ ≤ d f ∞ and c w ≤ b/a 3 , we obtain that

E f (2) h -fh 2 (x 0 ) ≤ 5db 2 4nha 6 f ∞ K 2 2 . (34)
Gathering ( 32), ( 33) and (34) yields the result. 6.2. Proof of Proposition 4.2. By (30) it follows that

E f (2) h -f 2 2 ≤ 3 E f (2) h -fh 2 2 + E fh -E[ fh ] 2 2 + E[ fh ] -f 2 2 . (35)
Concerning the last term of (35), proceeding as in [START_REF] Tsybakov | Introduction à l'estimation non-paramétrique[END_REF], we get

(36) E[ fh ] -f 2 2 ≤ Lh β (ℓ -1)! |K(u)||u| β du 2 .
For the first right-hand-side term of (35) we obtain

E f (2) h -fh 2 2 = E   1 n n i=1 w • Ĝn (Z i ) -w • G(Z i ) K h (x -Z i ) 2   dx ≤ E w • Ĝn (Z 1 ) -w • G(Z 1 ) 2 K 2 h (x -Z 1 ) dx ≤ c 2 w h K 2 2 E Ĝn (Z 1 ) -G(Z 1 ) 2 ≤ b 2 nha 6 K 2 2 , (37)
where we used that E[( Ĝn (Z 1 )-G(Z 1 )) 2 ] ≤ 1/n. This property can be shown by proceeding as in the pointwise case and using that G(Z 1 ) has uniform distribution.

For the second term of (35) we obtain

E fh -E[ fh ] 2 2 = 1 n E w • G(Y 1 ) (K h (x -Y 1 )) 2 dx ≤ 1 an (K h (x -y)) 2 dxf (y)dy = 1 anh K 2 2 . (38)
Combining ( 36), ( 37) and (38) completes the proof. 6.3. Proof of Theorem 4.1. For the sake of readability, super-indices (2) are omitted in the whole proof. For any h ∈ H,

fĥ (x 0 ) -f 2 (x 0 ) ≤ 3 fĥ (x 0 ) -fh, ĥ(x 0 ) 2 (x 0 ) + fh, ĥ(x 0 ) -fh 2 (x 0 ) + fh -f 2 (x 0 ) ≤ 3 A 0 (h, x 0 ) + V 0 ( ĥ(x 0 )) + A 0 ( ĥ(x 0 ), x 0 ) + V 0 (h) + fh -f 2 (x 0 ) ≤ 6A 0 (h, x 0 ) + 6V 0 (h) + 3 fh (x 0 ) -f (x 0 ) 2 , ( 39
)
where the second inequality holds by the definition of A 0 , i.e. for all h, h ′ ∈ H we have

A 0 (h, x 0 ) + V 0 (h ′ ) ≥ fh,h ′ (x 0 ) -fh ′ (x 0 ) 2 .
The last inequality holds by the definition of ĥ(x 0 ), that is A 0 ( ĥ(x 0 ), x 0 ) + V 0 ( ĥ

(x 0 )) ≤ A 0 (h, x 0 ) + V 0 (h) for all h ∈ H. The term E[( fh (x 0 ) -f (x 0 )) 2 ] is controlled by Proposition 4.1. Hence, it is sufficient to study the term E[A 0 (h, x 0 )]. We state that (40) A 0 (h, x 0 ) = sup h ′ ∈H fh,h ′ (x 0 ) -fh ′ (x 0 ) 2 -V 0 (h ′ ) + ≤ 5 (D 1 + D 2 + D 3 + D 4 + D 5 ) ,
where

D 1 = sup h ′ ∈H fh,h ′ (x 0 ) -fh,h ′ (x 0 ) 2 , D 2 = sup h ′ ∈H fh,h ′ (x 0 ) -E fh,h ′ (x 0 ) 2 - V 0 (h ′ ) 10 + , D 3 = sup h ′ ∈H E fh,h ′ (x 0 ) -E fh ′ (x 0 ) 2 , D 4 = sup h ′ ∈H E fh ′ (x 0 ) -fh ′ (x 0 ) 2 - V 0 (h ′ ) 10 + ,
and

D 5 = sup h ′ ∈H fh ′ (x 0 ) -fh ′ (x 0 ) 2 , with fh,h ′ = K h ′ * fh .
We start with term

D 3 . Recall that E[ fh (x 0 )] = K h * f (x 0 ) by (31). Likewise, by property (7), E fh,h ′ (x 0 ) = K h ′ * K h * f (x 0 ). In general we have s * r ∞ ≤ s ∞ r 1 and K h 1 = K 1 , yielding E fh,h ′ (x 0 ) -E fh ′ (x 0 ) = |K h ′ * (K h * f -f )(x 0 )| ≤ K h * f -f ∞ K 1 . Now remark that (K h * f -f ) (x) = b(x)
is the pointwise bias term considered in the proof of Proposition 4.1 Hence, (32) yields ( 41)

D 3 ≤ C 2 1 K 2 1 h 2β . Concerning term D 4 we note that E[D 4 ] ≤ h∈H E fh (x 0 ) -E fh (x 0 ) 2 - V 0 (h) 10 + = h∈H ∞ 0 P fh (x 0 ) -E fh (x 0 ) 2 - V 0 (h) 10 + > x dx = h∈H ∞ 0 P fh (x 0 ) -E fh (x 0 ) > V 0 (h) 10 + x dx .
The probability in the last term can be bounded by the Bernstein inequality. To this end we introduce the random variables

S i = w•G(Z i )K h (x 0 -Z i ). Obviously, |S i | ≤ K ∞ /(ah) =:
M almost surely and by property ( 7)

Var (S i ) ≤ E w 2 • G(Z 1 )K 2 h (x 0 -Z 1 ) = E w • G(Y 1 )K 2 h (x 0 -Y 1 ) ≤ K 2 2 f ∞ ah =: v .
Hence, the Bernstein inequality implies for any x > 0

P fh (x 0 ) -E fh (x 0 ) ≥ V 0 (h) 10 + x = P 1 n n i=1 (S i -E[S i ]) ≥ V 0 (h) 10 + x ≤ 2 max exp - n 4v V 0 (h) 10 + x , exp - nα 4M V 0 (h) 10 exp - n(1 -α) 4M √ x ,
for any α ∈ [0, 1] as √ • is a concave function. By the definition of V 0 (h)

n 4v V 0 (h) 10 = κ 0 K 2 1 log n 40 ≥ p log n , for κ 0 ≥ 40p, since K 2 1 ≥ 1. Furthermore, nα 4M V 0 (h) 10 = α K 2 K 1 4 K ∞ κ 0 a f ∞ hn log n 10 = ρ α hn log n ≥ p log n ,
for nh ≥ (p 2 /ρ 2 α ) log n. We can choose α ∈ [0, 1] sufficiently close to 0 such that p/ρ α > 1. Then the inequality in the last display holds under (H3), yielding

E[D 4 ] ≤ h∈H ∞ 0 2n -p max exp - nhax 4 K 2 2 f ∞ , exp - (1 -α)nha √ x 4 K ∞ dx ≤ 2n -p h∈H ∞ 0 max e -τ 1 nhx , e -τ 2 nh √ x dx ≤ 2n -p h∈H max 1 τ 1 , 2 τ 2 2 ≤ C ′ n -p+1 ,
as h ≥ 1/n and the cardinality of H verifies #H ≤ n. Finally, we choose p = 2 to get

(42) E[D 4 ] ≤ C ′ n .
Term D 2 can be treated in exactly the same way as D 4 . More precisely, instead of S i use

T i = w • G(Z i )K h * K h ′ (Z i -x 0 ) verifying fh,h ′ (x 0 ) -E fh,h ′ (x 0 ) = 1 n n i=1 T i -E[T i ] ,
and

|T i | ≤ K ∞ K 1 /(ah ′ ) =: M and Var(T 1 ) ≤ f ∞ K 2 1 K 2 2 /(ah ′ ) =: v. Hence, the Bernstein inequality yields (43) E[D 2 ] ≤ C ′′ n .
To study the terms D 5 and D 1 we first prove the following property.

Lemma 6.1. Under the assumptions (H2) and ( 5), for any set Ω and for all t ∈ R,

E fh (t) -fh (t) 2 1 Ω c ≤ c 2 w K 2 ∞ n 2 P (Ω c
) and

E fh ′ ,h (t) -fh ′ ,h (t) 2 1 Ω c ≤ c 2 w K 2 ∞ K 2 1 n 2 P (Ω c ) .
Proof. By using Ĝn -G ∞ ≤ 1, we have

E fh (t) -fh (t) 2 1 Ω c = E   1 n n i=1 (w • Ĝn (Z i ) -w • G(Z i ))K h (t -Z i ) 2 1 Ω c   ≤ c 2 w n 2 E   n i=1 |K h (t -Z i )| 2 1 Ω c   ≤ c 2 w E K 2 h (t -Z 1 )1 Ω c ≤ c 2 w K h 2 ∞ E [1 Ω c ] = c 2 w h 2 K 2 ∞ P (Ω c ) ≤ c 2 w K 2 ∞ n 2 P (Ω c ) ,
as 1/h ≤ n. In the same way, we show the second statement of the Lemma, by using

K h ′ * K h ∞ ≤ K h ′ ∞ K h 1 ≤ n K ∞ K 1 .
Now let Ω = {ω : Ĝn -G ∞ ≤ s} for some constant s > 0. Then (see [START_REF] Massart | The tight constant in the dvoretzky-kiefer-wolfowitz inequality[END_REF]), ( 44)

P(Ω c ) = P( Ĝn -G ∞ > s) ≤ e -2ns 2 ,
by the Dvoretzky-Kiefer-Wolfowitz inequality. This implies that

E sup h∈H fh (x 0 ) -fh (x 0 ) 2 1 Ω c ≤ h∈H E fh (x 0 ) -fh (x 0 ) 2 1 Ω c ≤ h∈H c 2 w K 2 ∞ n 2 e -2ns 2 = c 2 w K 2 ∞ n 3 e -2ns 2 < ∞ , as #H ≤ n. Furthermore, E sup h∈H fh (x 0 ) -fh (x 0 ) 2 1 Ω ≤ c 2 w E   sup h∈H 1 n n i=1 | Ĝn (Z i ) -G(Z i )||K h (x 0 -Z i )| 2 1 Ω   ≤ s 2 c 2 w E   sup h∈H 1 n n i=1 |K h (x 0 -Z i )| 2   ≤ 2s 2 c 2 w    E   sup h∈H 1 n n i=1 (|K h (x 0 -Z i )| -E[|K h (x 0 -Z i )|]) 2   + sup h∈H E 1 n n i=1 |K h (x 0 -Z i )| 2    ≤ 2s 2 c 2 w 1 n h∈H Var (|K h (x 0 -Z 1 )|) + sup h∈H [E(|K h (x 0 -Z 1 )|)] 2 .
On the one hand,

1 n h∈H Var (|K h (x 0 -Z 1 )|) ≤ 1 n h∈H E K 2 h (x 0 -Z 1 ) = 1 n h∈H 1 h K 2 2 g ∞ ≤ S K 2 2 d f ∞ ,
where S is defined in (H3) and d in (6). On the other hand,

sup h∈H [E(|K h (x 0 -Z 1 )|)] 2 = sup h∈H |K(z)|g(x 0 -zh)dz 2 ≤ d 2 f 2 ∞ K 2 1 . It follows that E[D 5 ] ≤ µ 1 n 3 e -2ns 2 + µ 2 s 2 , with constants µ 1 = c 2 w K 2 ∞ and µ 2 = 2c 2 w d f ∞ (S K 2 2 + d f ∞ K 2 1 ). Choosing s 2 = 2 log n/n gives (45) E[D 5 ] ≤ µ 1 n + 2µ 2 log n n .
Finally, the study of D 1 follows the same line as the study of D 5 . That is, on the one hand, we have for the same set Ω

E [D 1 1 Ω c ] ≤ c 2 w K 2 ∞ K 2 1 n 3 e -2ns 2 . On the other hand, E [D 1 1 Ω ] ≤ 2s 2 c 2 w 1 n h∈H E (K h ′ * K h (x 0 -Z 1 )) 2 + sup h∈H (E [|K h ′ * K h (x 0 -Z 1 )|]) 2 .
By the generalized Minkowski inequality, we obtain

E (K h ′ * K h (x 0 -Z 1 )) 2 ≤ |K h ′ (u)| K 2 h (x 0 -z -u)g(z)dz 1/2 du 2 ≤ g ∞ K h 2 2 K h ′ 2 1 ≤ d f ∞ K 2 2 K 2 1 /h . Furthermore, sup h∈H (E [|K h ′ * K h (x 0 -Z 1 )|]) 2 ≤ sup h∈H K h ′ (u) |K(v)| g(x 0 -vh -u)dvdu 2 ≤ (d f ∞ K 2 1 ) 2 . It follows with μ1 = µ 1 K 2 1 and μ2 = µ 2 K 2 1 that E[D 1 ] ≤ μ1 n 3 e -2ns 2 + μ2 s 2 . Hence, (46) E[D 1 ] ≤ μ1 n + 2μ 2 log n n ,
with s 2 = 2 log n/n. Now, if we plug (41), ( 42), ( 43), ( 45) and ( 46) into (40), we get

E[A 0 (h, x 0 )] ≤ C1 h 2β + C2 log n n ,
which, associated with Proposition 4.1, can be inserted in (39) to end the proof of Theorem 4.1.

6.4. Proof of Theorem 4.2. In all the proof below, super-indices (2) are omitted. Similar to the pointwise case, we have for any

h ∈ H fĥ -f 2 2 ≤ 6A(h) + 6V (h) + 3 fh -f 2 2 . (47) By the proof of Proposition 4.2, E fh -f 2 2 ≤ 3 f h -f 2 2 + C 4 nh . (48)
Hence, only term E[A(h)] needs to be studied. By analogy to the proof of Theorem 4.1,

A(h) ≤ 5(F 1 + F 2 + F 3 + F 4 + F 5 ) , (49)
where

F 1 = sup h ′ ∈H fh,h ′ -fh,h ′ 2 2 , F 2 = sup h ′ ∈H fh,h ′ -E[ fh,h ′ ] 2 2 - V (h ′ ) 10 + , F 3 = sup h ′ ∈H E[ fh,h ′ ] -E[ fh ′ ] 2 2 , F 4 = sup h ′ ∈H E[ fh ′ ] -fh ′ 2 2 - V (h ′ ) 10 + , F 5 = sup h ′ ∈H fh ′ -fh ′ 2 2 .
First, we study term

F 3 . The inequality u * v 2 ≤ u 1 v 2 yields (50) F 3 = sup h ′ ∈H K h ′ * K h * f -K h ′ * f 2 2 ≤ sup h ′ ∈H K h ′ 2 1 K h * f -f 2 2 = K 2 1 f -f h 2 2 .
To study term F 4 we introduce the centered empirical process ν n,h defined by

ν n,h (ψ) = fh -E[ fh ], ψ = 1 n n i=1 (w • G(Z i )K h (u -Z i ) -E[w • G(Z i )K h (u -Z i )]) ψ(u)du .
As ψ → ν n,h (ψ) is continuous, the supremum can be taken over a countable dense subset of {ψ ∈ L 2 , ψ = 1}, which we denote by B(1). Then, fh -

E[ fh ] 2 2 = sup ψ∈B(1) fh - E[ fh ], ψ 2 = sup ψ∈B(1) ν n,h (ψ). Therefore we obtain E[F 4 ] ≤ h∈H E fh -E[ fh ] 2 2 - V (h) 10 + = h∈H E sup ψ∈B(1) ν 2 n,h (ψ) - V (h) 10 + .
The expectation in the last term can be bounded by Talagrand's inequality (see Subsection 6.7). More precisely, to apply this result, we have to determine the values of the constants H, M and v. Denote

f ψ (z) = w • G(z)K h * ψ(z), so that ν n,h (ψ) = 1 n n i=1 (f ψ (Z i ) -E[f ψ (Z i )]
). First, for any ψ ∈ B(1) the Cauchy-Schwarz inequality gives

f ψ ∞ ≤ 1 a K h * ψ ∞ = 1 a sup z | |K h (• -z)|, |ψ| | ≤ K h 2 ψ 2 a ≤ K 2 a √ h =: M .
Next, we see that 38). Furthermore, let ε 2 = 1/2. To obtain 4H 2 = V (h)/10, we set H = V (h)/40. Lastly, for any ψ ∈ B(1) we show that by ( 7)

E sup ψ∈B(1) |ν n,h (ψ)| 2 ≤ E sup ψ∈B(1) ν 2 n,h (ψ) ≤ E fh -E[ fh ] 2 2 ≤ V (h) κ 1 =: H 2 , by (
Var (f ψ (Z)) ≤ E (w • G(Z)K h * ψ(Z)) 2 ≤ 1 a (K h * ψ(y)) 2 f (y)dy ≤ 1 a f ∞ K h * ψ 2 2 ≤ 1 a f ∞ K h 2 1 ψ 2 2 ≤ 1 a f ∞ K 2 1 =: v .
Finally, Talagrand's inequality yields

E sup ψ∈B(1) ν 2 n,h (ψ) - V (h) 10 + ≤ C1 n e -C2 /h + 1 nh e -C3 √ n ≤ C1 n e -C2 /h + C4 n ,
where Ck > 0, k = 1, . . . , 4 are constants depending on K, f ∞ and a. Consequently, (51)

E[F 4 ] ≤ C1 n h∈H e -C2 /h + C4 n ≤ C5 n ,
as #H ≤ n and h∈H e -C2 /h ≤ Σ(C ′ 2 ) under condition (26). In the same way we obtain for F 2

(52) E[F 2 ] ≤ C n .
Now let us turn to F 5 . We note that fhfh

2 2 ≤ 4 a 2 n 2 n i=1 |K h (u -Z i )| 2 du ≤ 4 a 2 n n i=1 K h (• -Z i ) 2 2 = 4 a 2 h K 2 2 ≤ 4n a 2 K 2 2 .
Therefore,

E sup h∈H fh -fh 2 2 1 Ω c ≤ 4n a 2 K 2 2 P(Ω c ) ,
where Ω = {ω : G -Ĝn ∞ ≤ s} as previously, and we recall that P(Ω c ) ≤ e -2ns 2 . Following the same line as for D 5 in the pointwise case and by choosing s = log n/n, we conclude that

(53) E[F 5 ] ≤ C ′ 1 n + C ′ 2 log n n .
For F 1 , we follow the same line as for F 5 to obtain

(54) E[F 1 ] ≤ K 2 1 C ′ 1 n + C ′ 2 log n n .
Consequently, plugging (50), ( 51), ( 52), ( 53) and ( 54) into (49) gives a bound of E[A(h)].

Combining this with (48), (47) and the definition of V (h) yields Theorem 4.2.

6.5. Proof of Proposition 4.3. Pythagoras formula yields f -fm

2 2 = f -f m 2 2 + f m - fm 2 2
. By definition of the orthogonal projection f m = 2m j=0 a j ϕ j and by using equality (7), we have a j = ϕ j , f = E(ϕ j (Y )) = E(ϕ j (Z 1 )w • G(Z 1 )). This, together with formula (13) implies that f m -fm

2 2 = 2m j=0 (a j -âj ) 2 . If we define (55) ν n (h) = 1 n n i=1 [h(Z i )w • G(Z i ) -E(h(Z i ) w • G(Z i ))], (56) R n (h) = 1 n n i=1 h(Z i )[w • Ĝn (Z i ) -w • G(Z i )] ,
then we get f m -fm 2 2 ≤ 2 2m j=0 (ν n (ϕ j ) 2 + R n (ϕ j ) 2 ). We have, on the one hand,

2m j=0 E(ν 2 n (ϕ j )) = 2m j=0 1 n Var(ϕ j (Z i )w • G(Z i )) ≤ 2m j=0 1 n E ϕ 2 j (Z 1 )(w • G(Z 1 )) 2 ≤ 1 n E   2m j=0 ϕ 2 j ∞ (w • G(Z 1 )) 2   ≤ D m n E[(w • G(Z 1 )) 2 ] ≤ 1 a 2 D m n , (57)
because the basis satisfies 2m j=0 ϕ 2 j = 2m + 1 = D m . On the other hand, we have

2m j=0 E(R 2 n (ϕ j )) ≤ 2m j=0 E   1 n n i=1 ϕ j (Z i )[w • Ĝn (Z i ) -w • G(Z i )] 2   ≤ 1 n n i=1 2m j=0 E ϕ 2 j (Z i )[w • Ĝn (Z i ) -w • G(Z i )] 2 ≤ c 2 w 2m j=0 E G -Ĝn 2 ∞ ϕ 2 j (Z i ) ≤ c 2 w D m E G -Ĝn 2 ∞ ≤ c 2 w D m n , (58)
with ( 5) and because of E G -Ĝn 2 ∞ ≤ 1/n (see e.g. Brunel and Comte, 2005, p. 462). By gathering all terms, we obtain the risk bound stated in Proposition 4.3. 6.6. Sketch of proof of Theorem 4.3. In the following, we omit the super index (2) . It is easy to see that fm = arg min t∈Sm γ n (t) for γ n (t) = t 2 -2n -1 n i=1 w • Ĝn (Z i )t(Z i ). Thus, we can write γ n (t)γ n (s) = tf 2 2sf 2 2 -2ν n (ts) -2R n (ts), where ν n and R n are defined by ( 55) and ( 56). By definition of f m we have for all m ∈ M n , γ n ( f m) + pen( m) ≤ γ n (f m ) + pen(m). This can be rewritten as f mf 2 2 ≤ f mf 2 2 + pen(m) + 2ν n ( f mf m )pen( m) + 2R n ( f mf m ). Using this and and that 2xy ≤ x 2 /θ + θy 2 for all nonnegative x, y, θ, we obtain [R n (t)] 2 . C/n by using Talagrand Inequality in a standard way (see e.g. Brunel et al., 2005). For the last term E sup t∈Sm n , t 2 =1 [R n (t)] 2 , we define Ω G by ( 59)

f -f m 2 2 ≤ f -f m 2 2 + pen(m) + 2ν n ( f m -f m ) -pen( m) + 2R n ( f m -f m ) f -f m 2 2 ≤ f -f m
As f m -f m 2 2 ≤ 2( f m -f 2 2 + f m -f 2 2 ), this yields 1 4 E[ f -f m 2 2 ] ≤ 7 4 f -f m 2 2 + 2pen(m) + 8E sup t∈Sm n , t 2 =1 [R n (t)] 2
Ω G = { √ n Ĝn -G ∞ ≤ log(n)}.
As in (44), we use [START_REF] Massart | The tight constant in the dvoretzky-kiefer-wolfowitz inequality[END_REF] and get ( 60)

P( √ n Ĝn -G ∞ ≥ λ) ≤ 2e -2λ 2 .
This implies that P(Ω c G ) ≤ 2/n 2 . Then we write that E sup t∈Sm n , t 2 =1 [R n (t)] 2 is less than

E sup t∈Sm n , t 2 =1 [R n (t)1 Ω G ] 2 + E sup t∈Sm n , t 2 =1 [R n (t)1 Ω c G ] 2 := R 1 + R 2 .
For the first term, we have

R 1 ≤ c 2 w E Ĝn -G 2 ∞ 1 Ω G E sup t∈Sm n ,
R 2 ≤ j E(R 2 n (ϕ j )1 Ω c ) ≤ c 2 w nE 1/2 ( Ĝn -G 4 ∞ )P 1/2 (Ω c G ) ≤ C n .
This yields E sup t∈Sm n , t 2 =1 [R n (t)] 2 ≤ C log 2 (n)/n. Finally we obtain that, for all

m ∈ M n , E[ f -f m 2 2 ] ≤ 7 f -f m 2
2 + 8pen(m) + Klog 2 (n)/n, which ends the proof.

6.7. The Talagrand inequality. The following result follows from the Talagrand concentration inequality given in [START_REF] Klein | Concentration around the mean for maxima of empirical processes[END_REF] and arguments in [START_REF] Birgé | Minimum contrast estimators on sieves: exponential bounds and rates of convergence[END_REF] (see the proof of their Corollary 2 page 354).

Lemma 6.2. (Talagrand's inequality) Let Y 1 , . . . , Y n be independent random variables, let

ν n,Y (f ) = (1/n) n i=1 [f (Y i ) -E(f (Y i ))
] and let F be a countable class of uniformly bounded measurable functions. Then for ǫ 2 > 0

E sup f ∈F |ν n,Y (f )| 2 -2(1 + 2ǫ 2 )H 2 + ≤ 4 K 1 v n e -K 1 ǫ 2 nH 2 v + 98M 2 K 1 n 2 C 2 (ǫ 2 ) e - 2K 1 C(ǫ 2 )ǫ 7 √ 2 nH M
, with C(ǫ 2 ) = √ 1 + ǫ 2 -1, K 1 = 1/6, and

sup f ∈F f ∞ ≤ M, E sup f ∈F |ν n,Y (f )| ≤ H, sup f ∈F 1 n n k=1
Var(f (Y k )) ≤ v.

By standard denseness arguments, this result can be extended to the case where F is a unit ball of a linear normed space, after checking that f → ν n (f ) is continuous and F contains a countable dense family.

  i = 1, 2. Then select model m(i) given by (21) m(1) = arg min m∈Mn ĝm 2 + pen (1) (m) , m(2) = arg min m∈Mn -f (2)

Figure 1 .

 1 Figure 1. (rescaled) MISE for all estimators as a function of κ for different distributions of f .

Figure 2 .

 2 Figure 2. Each boxplot represents the values of MISE*1000 of the six estimators computed on 1000 datasets for 5 different distributions with n = 500 and for µ = 0.1, 0.8, 2 (from left to right in the figure).

  the projection estimators. Concerning the pointwise kernel estimators, the MISE is rather sensitive to the value of κ (i) 0 . Slight changes may have a strong impact on the result. Furthermore, the estimators have quite the opposite behavior for the Beta B(5, 1) and the Laplace distribution. For the former the minimum of the MISE is attained at κ here). It is not clear which value of κ (i) 0 achieves the best compromise among all distributions. In the following we set κ Comparison of all six estimators. There are several factors potentially influencing the performance of the different estimators. In our simulation study we consider• two different sample sizes (n = 500 and n = 2000),• three levels of the Poisson parameter (µ = 0.1, µ = 0.8 and µ = 2),

Figure 3 .

 3 Figure 3. (a) Noise sample. (b) Fluorescence data and target density f . (c) Estimators for fluorescence data and target density f .

22

  + pen(m) + 2 f mf m 2 sup t∈S m+Sm, t 2 =1 |ν n (t)|pen( m) + 2 f mf m 2 sup t∈S m+Sm, t 2 =1|R n t∈S m +Sm, t 2 =1

.

  Then the term E sup t∈S m+Sm, t 2 =1 [ν n (t)] 2 -(pen(m) + pen( m))/4 + is bounded by

Table 1 .

 1 Mean MISE*1000 values for the six different estimators in 30 different settings.

	µ	Normal N (10, 1) distribution 0.1 0.8 2		Laplace L(10, 3) distribution 0.1 0.8 2
	n	500 2000	500 2000	500 2000	500 2000	500 2000	500 2000
	gl. ker L 3.08 2.76	3.09 2.59	3.23 2.67	1.76 0.940 1.14 0.579 1.22 0.416
	gl. ker P 3.09 2.77	3.37 2.80	5.60 4.66	1.93 1.07	2.00 1.44	2.28 1.42
	p. ker L 3.09 2.68	3.09 2.65	3.26 2.70	5.50 3.33	4.16 2.66	3.81 1.88
	p. ker P 3.10 2.68	3.32 2.83	5.62 4.70	5.70 3.43	5.07 3.51	4.41 2.77
	proj L	1.82 0.421 1.83 0.433 2.26 0.617 3.82 2.72	3.92 2.77	4.24 2.87
	proj P	1.83 0.426 2.25 0.482 4.03 1.18	3.98 2.89	7.21 6.20	26.2 27.4
	µ	0.1	Beta B(5, 1) distribution 0.8	2		0.1	Beta B(5, 2) distribution 0.8	2
	n	500 2000	500 2000	500 2000	500 2000	500 2000	500 2000
	gl. ker L 1431 1350	1444 1367	1589 1363	324 323	323 323	322 323
	gl. ker P 1302 1227	970 971	547 512		318 318	288 288	358 350
	p. ker L 258 232	266 236	279 246		334 332	334 334	334 336
	p. ker P 251 225	236 210	247 268		324 321	252 243	219 201
	proj L	314 272	317 274	339 279		28.7 11.2	28.5 11.5	34.9 13.8
	proj P	292 250	205 166	128 123	27.5 10.9 22.5 10.1 28.9 11.0
			µ	Exponential E(1) distribution 0.1 0.8 2
			n	500 2000 500	2000 500 2000
			gl. ker L 86.1 77.4 89.0 79.3 93.6 84.7
			gl. ker P 88.1 79.0 104.6 92.1 140 124
			p. ker L 89.0 83.7 91.4 86.1 95.7 90.8
			p. ker P 90.7 85.2 106	98.7 142 131
			proj L	62.0 53.8 62.1 53.9 64.9 54.5
			proj P	67.1 59.0 138	136 704 755

2 in (18) are approximated by Riemann-type discretization.

  It is proved inBrunel and Comte (2005) that E sup t∈Sm n , t 2 =1 |ν ′ n (t 2 )| ≤ C log(n) if the density of Z 1 is bounded and N n ≤ O( √ n) for the trigonometric basis. Moreover E(t 2 (Z 1 )) ≤ t 2 2 f ∞ /w 0 . We obtain R 1 ≤ C log 2 (n)/n.On the other hand, we have

								t 2 =1	(	1 n	n i=1	|t(Z i )|) 2
	≤ c 2 w	log(n) n	E	sup t∈Sm n , t 2 =1 (	1 n	n i=1	t 2 (Z i ))
	≤ 2c 2 w	log(n) n	E	sup t∈Sm n , t 2 =1	|ν ′ n (t 2 )| +	t∈Sm n , t 2 =1 sup	E(t 2 (Z 1 ))
	where ν ′ n (t) = 1 n	n i=1 (t(Z

i ) -E(t(Z 1 )).


Wiem Maalel
email: wiem.maalel@gmail.com

Kuang Zhou

Arnaud Martin
email: arnaud.martin@univ-rennes1.fr

Zied Elouedi Larodec

Belief Hierarchical Clustering

Keywords: Clustering, Hierarchical clustering, Belief function, Belief clustering

In the data mining field many clustering methods have been proposed, yet standard versions do not take into account uncertain databases. This paper deals with a new approach to cluster uncertain data by using a hierarchical clustering defined within the belief function framework. The main objective of the belief hierarchical clustering is to allow an object to belong to one or several clusters. To each belonging, a degree of belief is associated, and clusters are combined based on the pignistic properties. Experiments with real uncertain data show that our proposed method can be considered as a propitious tool.

Introduction

Due to the increase of imperfect data, the process of decision making is becoming harder. In order to face this, the data analysis is being applied in various fields.

Clustering is mostly used in data mining and aims at grouping a set of similar objects into clusters. In this context, many clustering algorithms exist and are categorized into two main families: The first family involves the partitioning methods based on density such as kmeans algorithm [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF] that is widely used thanks to its convergence speed. It partitions the data into k clusters represented by their centers. The second family includes the hierarchical clustering methods such as the top-down and the Hierarchical Ascendant Clustering (HAC) [START_REF] Hastie | The elements of statistical learning; data mining, inference and prediction[END_REF]. This latter consists on constructing clusters recursively by partitioning the objects in a bottom-up way. This process leads to good result visualizations. Nevertheless, it has a non-linear complexity.

All these standard methods deal with certain and precise data. Thus, in order to facilitate the decision making, it would be more appropriate to handle uncertain data. Here, we need a soft clustering process that will take into account the possibility that objects belong to more than one cluster.

In such a case, several methods have been established. Among them, the Fuzzy C-Means [START_REF] Bezdek | Fcm: The fuzzy c-means clustering algorithm[END_REF] which consists in assigning a membership to each data point corresponding to the cluster center, and the weights minimizing the total weighted mean-square error. This method constantly converges. Patently, Evidential c-Means (ECM) [START_REF] Denoeux | EVCLUS: Evidential Clustering of Proximity Data[END_REF], [START_REF] Masson | Clustering interval-valued proximity data using belief functions[END_REF] is deemed to be a very fateful method. It enhances the FCM and generates a credal partition from attribute data. This method deals with the clustering of object data. Accordingly, the belief k-Modes method [START_REF] Ben Hariz | Clustering approach using belief function theory[END_REF] is a popular method, which builds K groups characterized by uncertain attribute values and provides a classification of new instances. Schubert has also found a clustering algorithm [START_REF] Schubert | Clustering belief functions based on attracting and conflicting metalevel evidence[END_REF] which uses the mass on the empty set to build a classifier.

Our objective in this paper is to develop a belief hierarchical clustering method, in order to ensure the membership of objects in several clusters, and to handle the uncertainty in data under the belief function framework.

This remainder is organized as follows: in the next section we review the ascendant hierarchical clustering, its concepts and its characteristics. In section 3, we recall some of the basic concepts of belief function theory. Our method is described in section 4 and we evaluate its performance on a real data set in section 5. Finally, Section 6 is a conclusion for the whole paper.

Ascendant hierarchical clustering

This method consists on agglomerating the close clusters in order to have finally one cluster containing all the objects x j (where j = 1, .., N ). Let's consider

P K = {C 1 , ..., C K } the set of clusters. If K = N , C 1 = x 1 , ..., C N = x N .
Thereafter, throughout all the steps of clustering, we will move from a partition P K to a partition P K-1 . The result generated is described by a hierarchical clustering tree (dendrogram), where the nodes represent the successive fusions and the height of the nodes represents the value of the distance between two objects which gives a concrete meaning to the level of nodes conscripted as "indexed hierarchy". This latter is usually indexed by the values of the distances (or dissimilarity) for each aggregation step. The indexed hierarchy can be seen as a set with an ultrametric distance d which satisfies these properties:

i) x = y ⇐⇒ d(x, y) = 0. ii) d(x, y) = d(y, x). iii) d(x, y) ≤ d(x, z) + d(y, z), ∀x, y, z ∈ IR.
The algorithm is as follows:

-Initialisation: the initial clusters are the N-singletons. We compute their dissimilarity matrix. -Iterate these two steps until the aggregation turns into a single cluster:

• Combine the two most similar (closest) elements (clusters) from the selected groups according to some distance rules. • Update the matrix distance by replacing the two grouped elements by the new one and calculate its distance from each of the other classes.

Once all these steps completed, we do not recover a partition of K clusters, but a partition of K -1 clusters. Hence, we had to point out the aggregation criterion (distance rules) between two points and between two clusters. We can use the Euclidian distance between N objects x defined in a space IR. Different distances can be considered between two clusters: we can consider the minimum as follows:

d(C i j , C i j ′ ) = min x k ∈C i j ,x k ′ ∈C i j ′ d(x k , x k ′ ) (1)
with j, j ′ = 1, ..., i. The maximum can also be considered, however, the minimum and maximum distances create compact clusters but sensitive to "outliers". The average can also be used, but the most used method is Ward's method, using Huygens formula to compute this:

∆I inter(C i j ,C i j ′ ) = m Cj m C j ′ m Cj + m C j ′ d 2 (C i j , C i j ′ ) (2)
where m Cj and m C j ′ are numbers of elements of C j and C j ′ respectively and C i j , C i j ′ the centers. Then, we had to find the couple of clusters minimizing the distance:

(C k l , C k l ′ ) = d(C i l , C i l ′ ) = min C i j ,C i l ′ ∈C i d(C i j , C i j ′ ) (3)
3 Basis on the theory of belief functions

In this Section, we briefly review the main concepts that will be used in our method that underlies the theory of belief functions [START_REF] Shafer | Mathematical Theory of evidence[END_REF] as interpreted in the Transferable Belief Model (TBM) [START_REF] Smets | The Transferable Belief Model[END_REF]. Let's suppose that the frame of discernment is Ω = {ω 1 , ω 2 , ..., ω 3 }. Ω is a finite set that reflects a state of partial knowledge that can be represented by a basis belief assignment defined as:

m : 2 Ω → [0, 1] A⊆Ω m(A) = 1 (4)
The value m(A) is named a basic belief mass (bbm) of A. The subset A ∈ 2 Ω is called focal element if m(A) > 0. One of the important rules in the belief theory is the conjunctive rule which consists on combining two basic belief assignments m 1 and m 2 induced from two distinct and reliable information sources defined as:

m 1 ∩ m 2 (C) = A∩B=C m 1 (A) • m 2 (B), ∀C ⊆ Ω (5)
The Dempster rule is the normalized conjunctive rule:

m 1 ⊕ m 2 (C) = m 1 ∩ m 2 (C) 1 -m 1 ∩ m 2 (∅) , ∀C ⊆ Ω (6)
In order to ensure the decision making, beliefs are transformed into probability measures recorded BetP, and defined as follows [START_REF] Smets | The Transferable Belief Model[END_REF]:

BetP(A) = B⊆Ω |A ∩ B | | B | m(B) (1 -m(∅))
, ∀A ∈ Ω (7)

Belief hierarchical clustering

In order to set down a way to develop a belief hierarchical clustering, we choose to work on different levels: on one hand, the object level, on the other hand, the cluster level. At the beginning, for N objects we have, the frame of discernment is Ω = {x 1 , ..., x N } and for each object belonging to one cluster, a degree of belief is assigned. Let P N be the partition of N objects. Hence, we define a mass function for each object x i , inspired from the k-nearest neighbors [START_REF] Denoeux | A k-Nearest Neighbor Classification Rule Based on Dempster-Shafer Theory[END_REF] method which is defined as follows:

m Ωi i (x j ) = αe -γd 2 (xi,xj ) m Ωi i (Ω i ) = 1 -m Ωi i (x j ) (8)
where i = j, α and γ are two parameters we can optimize [START_REF] Zouhal | An Evidence-Theoric k-NN Rule with Parameter Optimization[END_REF], d can be considered as the Euclidean distance, and the frame of discernment is given by

Ω i = {x 1 , ..., x N } \ {x i }.
In order to move from the partition of N objects to a partition of N -1 objects we have to find both nearest objects (x i , x j ) to form a cluster. Eventually, the partition of N -1 clusters will be given by P N -1 = {(x i , x j ), x k } where k = 1, ..., N \ {i, j}. The nearest objects are found considering the pignistic probability, defined on the frame Ω i , of each object x i , where we proceed the comparison by pairs, by computing firstly the pignistic for each object, and then we continue the process using arg max. The nearest objects are chosen using the maximum of the pignistic values between pairs of objects, and we will compute the product pair one by one.

(x i , x j ) = arg max xi,xj ∈P N BetP Ωi i (x j ) * BetP Ωj j (x i ) (9)
Then, this first couple of objects is a cluster. Now consider that we have a partition P K of K clusters {C 1 , . . . , C K }. In order to find the best partition P K-1 of K -1 clusters, we have to find the best couple of clusters to be merged. First, if we consider one of the classical distances d (single link, complete link, average, etc), presented in section 2, between the clusters, we delineate a mass function, defined within the frame Ω i for each cluster C i ∈ P K with C i = C j by:

m Ωi i (C j ) = αe -γd 2 (Ci,Cj ) (10)
m Ωi i (Ω i ) = 1 - m Ωi i (C j ) (11)
where

Ω i = {C 1 , . . . , C K } \ {C i }.
Then, both clusters to merge are given by:

(C i , C j ) = arg max Ci,Cj ∈P K BetP Ωi (C j ) * BetP Ωj (C i ) ( 12
)
and the partition P K-1 is made from the new cluster (C i , C j ) and all the other clusters of P K . The point by doing so is to prove that if we maximize the degree of probability we will have the couple of clusters to combine. Of course, this approach will give exactly the same partitions than the classical ascendant hierarchical clustering, but the dendrogram can be built from BetP and the best partition (i.e. the number of clusters) can be preferred to find. The indexed hierarchy will be indexed by the sum of BetP which will lead to more precise and specific results according to the dissimilarity between objects and therefore will facilitate our process. Hereafter, we define another way to build the partition P K-1 . For each initial object x i to classify, it exists a cluster of P K such as x i ∈ C k . We consider the frame of discernment Ω i = {C 1 , . . . , C K } \ {C k }, m, which describes the degree that the two clusters could be merged, can be noted m Ω and we define the mass function:

m Ωi i (C kj ) = xj ∈C k j αe -γd 2 (xi,xj ) (13)
m Ωi i (Ω i ) = 1 - xj ∈C k j m Ωi i (C kj ) (14)
In order to find a mass function for each cluster C i of P K , we combine all the mass functions given by all objects of C i by a combination rule such as the Dempster rule of combination given by equation [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF]. Then, to merge both clusters we use the equation ( 12) as before. The sum of the pignisitic probabilities will be the index of the dendrogram, called BetP index.

Experimentations

Experiments were first applied on diamond data set composed of twelve objects as describe in Figure 1.a and analyzed in [START_REF] Masson | Clustering interval-valued proximity data using belief functions[END_REF]. The dendrograms for both classical and Belief Hierarchical Clustering (BHC) are represented by Figures 1.b and 1.c. The object 12 is well considered as an outlier with both approaches. With the belief hierarchical clustering, this object is clearly different, thanks to the pignistic probability. For HAC, the distance between object 12 and other objects is small, however, for BHC, there is a big gap between object 12 and others. This points out that our method is better for detecting outliers. If the objects 5 and 6 are associated to 1, 2, 3 and 4 with the classical hierarchical clustering, with BHC these points are more identified as different. This synthetic data set is special because of the equidistance of the points and there is no uncertainty. We continue our experiments with a well-known data set, Iris data set, which is composed of flowers from four types of species of Iris described by sepal length, sepal width, petal length, and petal width. The data set contains three clusters known to have a significant overlap.

In order to reduce the complexity and present distinctly the dendrogram, we first used the k-means method to get initial few clusters for our algorithm. It is not necessary to apply this method if the number of objects is not high.

Several experiments have been used with several number of clusters. We present in Figure 2 the obtained dendrograms for 10 and 13 clusters. We notice different combinations between the nearest clusters for both classical and belief hierarchical clustering. The best situation for BHC is obtained with the pignistic equal to 0.5 because it indicates that the data set is composed of three significant clusters which reflects the real situation. For the classical hierarchical clustering the results are not so obvious. Indeed, for HAC, it is difficult to decide for the optimum cluster number because of the use of the euclidean distance and as seen in Figure 2.c it indistinguishable in terms of the y-value. However, for BHC, it is more easy to do this due to the use of the pignistic probability.

In order to evaluate the performance of our method, we use some of the most popular measures: precision, recall and Rand Index (RI). The results for both BHC and HAC are summarized in Table 1. The first three columns are for BHC, while the others are for HAC. In fact, we suppose that F c represents the final number of clusters and we start with F c = 2 until F c = 6. We fixed the value of k init at 13. We note that for F c = 2 the precision is low while the recall is of high value, and that when we have a high cluster number (F c = 5 or 6), the precision will be high but the recall will be relatively low. Thus, we note that for the same number of final clusters (e.g. F c = 4), our method is better in terms of precision, recall and RI. Tests are also performed to a third data base, Congressional Voting Records Data Set. The results presented in Figure 3 show that the pignistic probability value increased at each level, having thereby, a more homogeneous partition. We notice different combinations, between the nearest clusters, that are not the same within the two methods compared. For example, cluster 9 is associated to cluster 10 and then to 6 with HAC, but, with BHC it is associated to cluster 4 and then to 10. Although, throughout the BHC dendrograms shown in Figure 3.c and Figure 3.d, the best situation indicating the optimum number of clusters can be clearly obtained. This easy way is due to the use of the pignistic probability. For this data set, we notice that for F c = 2 and 3, the precision is low while the recall is high. However, with the increasing of our cluster number, we notice that BHC provides a better results. In fact, for F c = 3, 4, 5 and 6 the precision and RI values relative to BHC are higher then the precision and RI values relative to HAC, which confirmed the efficiency of our approach which is better in terms of precision and RI.

Conclusion

Ultimately, we have introduced a new clustering method using the hierarchical paradigm in order to implement uncertainty in the belief function framework. This method puts the emphasis on the fact that one object may belong to several clusters. It seeks to merge clusters based on its pignistic probability. Our method was proved on data sets and the corresponding results have clearly shown its efficiency. The algorithm complexity has revealed itself as the usual problem of the belief function theory. Our future work will be devoted to focus on this peculiar problem.

Fig. 1 .

 1 Fig. 1. Clustering results for Diamond data set.

Fig. 2 .

 2 Fig. 2. Clustering results on IRIS data set for both hierarchical (HAC) (Fig. a and b) and belief hierarchical (BHC) (Fig. c and d) clustering (Kinit is the cluster number by k-means first).

Fig. 3 .

 3 Fig. 3. Clustering results on Congressional Voting Records Data Set for both hierarchical and belief hierarchical clustering.

Table 1 .

 1 Evaluation results

		BHC		HAC
		Precision Recall	RI Precision Recall	RI
	Fc = 2	0.5951 1.0000 0.7763	0.5951 1.0000 0.7763
	Fc = 3	0.8011 0.8438 0.8797	0.6079 0.9282 0.7795
	Fc = 4	0.9506 0.8275 0.9291	0.8183 0.7230 0.8561
	Fc = 5	0.8523 0.6063 0.8360	0.8523 0.6063 0.8360
	Fc = 6	0.9433 0.5524 0.8419	0.8916 0.5818 0.8392

Table 2 .

 2 Evaluation results for Congressional Voting Records Data Set

		BHC		HAC
		Precision Recall	RI Precision Recall	RI
	Fc = 2	0.3873 0.8177 0.5146	0.5951 1.0000 0.7763
	Fc = 3	0.7313 0.8190 0.8415	0.6288 0.8759 0.7892
	Fc = 4	0.8701 0.6833 0.8623	0.7887 0.7091 0.8419
	Fc = 5	0.8670 0.6103 0.8411	0.7551 0.6729 0.8207
	Fc = 6	0.9731 0.6005 0.8632	0.8526 0.6014 0.8347