Sahar Slama
email: sahar.slama@essths.u-sousse.tn

Yousri Slaoui
email: yousri.slaoui@math.univ-poitiers.fr

Hamdi Fathallah
email: hamdi.fathallah@essths.u-sousse.tn

Statistical inference for multivariate conditional cumulative distribution function estimation by stochastic approximation method

Keywords: conditional cumulative distribution function (CCDF), kernel estimation method, stochastic approximation algorithm, plug-in principle AMS 2010 subject classications 62G20, 62H12, 62L20, 65D10

This paper handles non-parametric estimation of a conditional cumulative distribution function (CCDF). Using a recursive approach, we set forward a multivariate recursive estimator dened by stochastic approximation algorithm. Our basic objective is to investigate the statistical inference of our estimator and compare it with that of non-recursive Nadaraya-Watson's estimator. From this perspective, we rst derive the asymptotic properties of the proposed estimator which highly depend on the choice of two parameters, the stepsize (γ n ) as well as the bandwidth (h n ). The second generation plug-in method, a method of bandwidth selection minimizing the Mean Weighted Integrated Squared Error (MW ISE) of the estimator in reference, entails the optimal choice of the bandwidth and therefore maintains an appropriate choice of the stepsize parameter. Basically, we demonstrate that, under some conditions, the Mean Squared Error (MSE) of the proposed estimator can be smaller than the one of Nadaraya Watson's estimator. We corroborate our theoretical results through simulation studies and two real dataset applications, namely the Insurance Company Benchmark (COIL 2000) dataset as well as the French Hospital Data of COVID-19 epidemic.

Introduction

Assume that we observe independent identically distributed vectors (X 1 , Y 1 ), ..., (X n , Y n ) of a bivariate random variable (X, Y ) with common cumulative distribution function π(x, y)

where one is interested in modeling the functional dependence of the observation Y on the covariable X by the conditional cumulative distribution function (CCDF) of Y given X = x, denoted by for all real y and x π(y|x)

:= P [Y y|X = x] .
We shall also assume that the bivariate random variable (X, Y ) (resp. the random variable X) has a density function f (X,Y ) (resp. f X ) with respect to the Lebesgue measure. Recall that for all real y and x such that f X (x) = 0, the CCDF of Y given X = x is expressed by

π(y|x) = R 1 {u y} f (X,Y ) (x, u)du f X (x) .
In a variety of non-parametric statistical problems, the estimation of a CCDF is a key aspect of inference. Remember that the CCDF has the merit of characterizing the conditional law of the considered bivariate random variables. Notably, the CCDF is often useful in reliability or survival analysis. More specically, the conditional survival function S (y|x) dened by, for all real y and x, S (y|x) := 1 -π(y|x) is of extreme interest, either by itself, or by its independence with the conditional hazard function h (y|x) indicated by, for all real y and x, h (y|x) := f (y|x) S (y|x) where f (y|x) denotes the conditional density of Y given X = x. Furthermore, conditional quantiles can also be deduced from the CCDF π by (pseudo)-inversion given x of the function y → π(y|x) and the same procedure may be applied to the estimator of CCDF to nd conditional quantile estimators. Several non-parametric estimators have been elaborated to estimate the CCDF. Many of them rely initially on estimating the R 1 {u y} f (X,Y )(x,u) du. The conditional cumulative distribution function was rst extensively explored by [START_REF] Stute | Conditional empirical processes[END_REF] using a nearest-neighbor-type conditional empirical process. Subsequently, [START_REF] Hall | Methods for estimating a conditional distribution function[END_REF], motivated by the problem of setting prediction intervals in time series analysis, developed a new non-parametric method for CCDF estimation resting on an adjusted form of NadarayaWatson estimator. Afterwards, [START_REF] Ferrigno | Certainty bands for the conditional cumulative distribution function and applications[END_REF] established uniform asymptotic certainty bands for the CCDF using the same strategy. For a general non-parametric regression model, [START_REF] Kiwitt | Estimating the conditional error distribution in nonparametric regression[END_REF] set up two estimators using a kernel approch, where the distribution of the error given the covariate is modeled by a CCDF provided by P ( ≤ y|X = x). On a given compact set, [START_REF] Brunel | Minimax estimation of the conditional cumulative distribution function[END_REF] constructed a minimax estimator of the CCDF. Thereafter, [START_REF] Veraverbeke | Preadjusted non-parametric estimation of a conditional distribution function[END_REF] built up a new estimator of CCDF investing a method of pre-adjusting the original observations non-parametrically. Recently, [START_REF] Bouanani | Local linear conditional cumulative distribution function with mixing data[END_REF] introduced a new method to settle CCDF estimation problem based on local polynomial technique. Many functional estimations are grounded on estimating the CCDF see e.g [START_REF] Laksaci | Conditional cumulative distribution estimation and its applications[END_REF], [START_REF] Laksaci | Note on the functional linear estimate of conditional cumulative distribution function[END_REF] and [START_REF] Almanjahie | kNN local linear estimation of the conditional cumulative distribution function: Dependent functional data case[END_REF]. The CCDF is involved in a wide range of applications, for instance, in medicine see [START_REF] Gannoun | Reference curves based on non-parametric quantile regression[END_REF], econometrics see [START_REF] Li | Optimal bandwidth selection for nonparametric conditional distribution and quantile functions[END_REF] or machine learning domain see the recent work of [START_REF] Chilinski | Neural likelihoods via cumulative distribution function. roeedings of the QTth gonferene on nertinty in erti[END_REF]. In a broader context, extensive state of art works including various non-parametric approaches tackled the conditional estimation. We can state for example [START_REF] Yu | Local linear quantile regression[END_REF], [START_REF] Fan | Estimation of conditional densities and sensitivity measures in nonlinear dynamical systems[END_REF], Berlinet et al. (1998a,b), Honda (2000a,b) and [START_REF] Plancade | Adaptive estimation of the conditional cumulative distribution function from current status data[END_REF]. For recent references see [START_REF] Benziadi | Recursive kernel estimate of the conditional quantile for functional ergodic data[END_REF], [START_REF] Chowdhury | Nonparametric estimation of conditional distribution functions with longitudinal data and time-varying parametric models[END_REF], [START_REF] Al-Awadhi | Functional data analysis: local linear estimation of the L 1 -conditional quantiles[END_REF], [START_REF] Chikr-Elmezouar | FDA: strong consistency of the kNN local linear estimation of the functional conditional density and mode[END_REF] and [START_REF] Slaoui | Adaptive recursive kernel conditional density estimators under censoring data[END_REF].

Over the past decade, data streams have become an increasingly important area of research. Some of the most common data streams include Internet packet data, Twitter activity, Facebook newsfeed, credit card transactions and more recently COVID-19 epidemic data. In these situations, the data arrives regularly so that it is impossible to store it in a traditional database. In such a context, it is very interesting to build a recursive multivariate conditional cumulative distribution estimator that does not need to store all the data in memory and that can be easily updated to handle the online data. The basic target of this paper is to provide a non-parametric strategy to recursively estimate the CCDF. The paper is organized as follows. In section 2 we introduce our method for the estimation procedure of CCDF. The main results of our recursive estimator are displayed in section 3. Section 4 exhibits the asymptotic properties of non-recursive Nadaraya-Watson's estimator. A second generation plug-in scheme for constructing data-driven bandwidth selection procedures is elaborated in section 5. Section 6 highlights the performances of our estimator on simulated data as well as a real dataset. Finally, the conclusion is drawn in section 7.

Preliminaries

Let (X, Y ) ∈ R d × R q and (X 1 , Y 1 ), ..., (X n , Y n ) be independent random vectors identically distributed as (X, Y ) with joint density function f (X,Y ) and let f X denote the probability density of X. In this paper, our central focus is upon the problem of estimating the CCDF of Y given X = x provided by

π : R q × R d -→ R (y|x) -→ P[Y y|X = x] = R q 1 {u y} f (X,Y ) (x, u)du f X (x) .
Here, we introduce our recursive estimator π n specied by

π n (y|x) =    a n (x, y) f n (x) if f n (x) = 0 0 otherwise , (1)
with

a n (x, y) = Π n n k=1 Π -1 k γ k χ k (y)h -d k K x -X k h k and f n (x) = Π n n k=1 Π -1 k γ k h -d k K x -X k h k ,
where

K is a multivariate kernel satisfying R d K(t)dt = 1.
χ is a multivariate indicator function identied by χ k : R q -→ R, y -→ 1 {Y k y} .

(h n ) is the bandwidth: a sequence of positive real numbers that tends to zero.

(γ n ) is the stepsize: a positive sequence of real numbers decreasing towards zero and

Π n = n j=1
(1 -γ j ).

Our main purpose is to examine the asymptotic properties of the proposed multivariate estimator of the CCDF and to corroborate its performances.

We shall compare our estimator to the generalized kernel CCDF estimator of Nadaraya-Watson Nadaraya (1964) and [START_REF] Watson | Smooth regression analysis[END_REF] π n expressed by

π n (y|x) =    a n (x, y) f n (x) if f n (x) = 0 0 otherwise , (2)
with

a n (x, y) = 1 nh d n n k=1 χ k (y) K x -X k h n and f n (x) = 1 nh d n n k=1 K x -X k h n .
The recursive estimator was constructed based on dint of stochastic approximation method. Indeed, incorporating stochastic approximation algorithms in the context of non-parametric statistics dates back to the papers of [START_REF] Robbins | A stochastic approximation method[END_REF] and [START_REF] Kiefer | Stochastic estimation of the maximum of a regression function[END_REF]. Their research works have been extended in several directions. We refer the reader to [START_REF] Blum | Multidimensional stochastic approximation methods[END_REF], [START_REF] Nadaraya | On estimating regression[END_REF], [START_REF] Ruppert | Almost sure approximations to the Robbins-Monro and Kiefer-Wolfowitz processes with dependent noise[END_REF][START_REF] Duo | ndom itertive models, volume 34 of epplitions of wthemtis @xew orkA[END_REF]. Subsequently, in the paper of [START_REF] Mokkadem | The stochastic approximation method for the estimation of a multivariate probability density[END_REF], the multidimensional case was investigated in order to estimate of a multivariate probability density using the estimation by condence intervals. Additionally, [START_REF] Slaoui | The stochastic approximation method for estimation of a distribution function[END_REF] reused stochastic approximation methods to build an univariate distribution function estimator and lately [START_REF] Slama | Multivariate distribution function estimation using stochastic approximation method[END_REF] elaborated the multidimensional case in order to enhance the qualities of the multivariate distribution function estimator. Following the same recursive approach, we intend to establish a multivariate conditional cumulative distribution function estimator.

To build up a stochastic algorithm, which approaches the function

a : (x, y) -→ R q 1 {u y} f (X,Y ) (x, u)du
at a given couple of vectors (x, y), we dene an algorithm of search of the zero function φ : z -→ a(x, y) -z and we set:

(i) a 0 (x, y) ∈ R (ii) for all n ≥ 1, a n (x, y) = a n-1 (x, y) + γ n U n (x, y),
where U n (x, y) corresponds to an observation of the function φ at the point a n-1 (x, y).

Note that to dene U n (x, y), we adopt the approach of [START_REF] Révész | How to apply the method of stochastic approximation in the nonparametric estimation of a regression function[END_REF] and [START_REF] Tsybakov | Recurrent estimation of the mode of a multidimensional distribution[END_REF].

By considering U n (x, y) = χ n (y) h -d n K
x -X n h n -a n-1 (x, y), the stochastic approximation algorithm that is devoted to estimate recursively the function a at a couple of vectors (x, y) can be stated as follows :

a n (x, y) = (1 -γ n )a n-1 (x, y) + γ n χ n (y) h -d n K x -X n h n . (3)
Throughout this paper, we consider that a 0 (x, y) = 0. Therefore, by recurrence, we get

a n (x, y) = Π n n k=1 Π -1 k γ k χ k (y)h -d k K x -X k h k . ( 4
)
Within this framework, we use the recursive multivariate probability density estimator of the density function f X noted f n and dened in [START_REF] Mokkadem | The stochastic approximation method for the estimation of a multivariate probability density[END_REF]. It was constructed with the same tools of stochastic approximation algorithm and under the condition that f 0 (x) = 0, we have:

f n (x) = Π n n k=1 Π -1 k γ k h -d k K x -X k h k . ( 5
)
The following notations are highly useful as they are invested throughout the whole paper. * ∀ i, j ∈ {1, . . . , d}, a

i (., y) := ∂a ∂x i (., y), π

i (y|.) :=

∂π ∂x i (y|.), f (1)
X i (.) := ∂f X ∂x i (.), a (2)
ij (., y) := ∂ 2 a ∂x i ∂x j (., y), π

(2)

ij (y|.) := ∂ 2 π ∂x i ∂x j (y|.), f (2)
X ij (.) := ∂ 2 f X ∂x i ∂x j (.), µ i (K) := R d z 2 i K(z)dz. * ξ := lim n→+∞ (nγ n ) -1 . * R(K) := R d K 2 (z) dz. * I 1 := R d+q d j=1 µ j (K) π (2) jj (y|x)f X (x) + 2π
(1) j (y|x)f

(1)

X j (x) 2 f X,Y (x, y)dxdy. * I 2 := R d+q π(y|x) (1 -π(y|x)) f X (x)f X,Y (x, y)dxdy. * ∀ x ∈ R d , y ∈ R q , Z n (x, y) := h -d n χ n (y) K x -X n h n and W n (x) := h -d n K
x -X n h n .

In the sequel, let us present the following denition of class of regularly varying sequences introduced by [START_REF] Galambos | Regularly varying sequences[END_REF].

Denition 1. vet (v n ) e nonrndom positive sequene nd γ ∈ RF e sy tht

(v n ) ∈ GS(γ) if lim n→+∞ n 1 - v n-1 v n = γ.
In what follows, we introduce a lemma that will be widely invested throughout the study of our estimator π n . It is worth noting that the proof of this lemma was recorded in [START_REF] Mokkadem | The stochastic approximation method for the estimation of a multivariate probability density[END_REF].

Lemma 1. vet

(v n ) ∈ GS(v * ), (γ n ) ∈ GS(-α) nd let m > 0 suh tht m -v * ξ > 0. henD lim n→+∞ v n Π m n n k=1 Π -m k γ k v k = 1 m -v * ξ .
woreoverD for ny positive sequene (α n ) suh tht lim

n→+∞ α n = 0 nd ll C ∈ RD lim n→+∞ v n Π m n n k=1 Π -m k γ k v k α k + C = 0.
In order to introduce our theoretical main results, we need the following technical assumptions.

Assumptions:

(A 1 ) K : R d -→ R is a continuous bounded function satisfying:

R d K(u)du = 1 , ∀j ∈ {1, . . . , d}, R d u j K(u)du = 0 and R d u 2 j K(u)du < ∞. (A 2 ) (i) (γ n ) ∈ GS(-α), with α ∈ 1 2 , 1 . (ii) (h n ) ∈ GS(-a), with a ∈ (0, 1]. (iii) lim n→+∞ (nγ n ) ∈ min{2a, α-ad 2 }, ∞ . (A 3 ) (i) f X is
bounded, twice dierentiable with respect to x and for all i, j ∈ {1, . . . , d}, f

(2) X ij is bounded.

(ii) a is bounded, twice dierentiable with respect to x and for all i, j ∈ {1, . . . , d}, a

(2) ij is bounded.

(iii) f (X,Y ) is bounded and twice continuously dierentiable with respect to x.

Discussion of the assumptions:

All these assumptions are standard and are generally assumed within the context of nonparametric estimation. Classical assumption (A 1 ) provides regularity conditions on the kernel density estimator introduced by [START_REF] Rosenblatt | Remarks on some nonparametric estimates of a density function[END_REF] and [START_REF] Parzen | On estimation of a probability density function and mode[END_REF].It is widely used in the non-parametric framework for the functional estimation. Assumption (A 2 ) on the stepsize and the bandwidth was used in the recursive framework for the estimation of the density function in [START_REF] Mokkadem | The stochastic approximation method for the estimation of a multivariate probability density[END_REF]; Slaoui (2014a) and for the distribution function estimation in [START_REF] Slaoui | The stochastic approximation method for estimation of a distribution function[END_REF]. Furthermore, assumption (A 2 )(iii) on the limit of (nγ n ) as n goes to innity is frequent in the framework of stochastic approximation algorithms. It implies, in particular, that the limit of (nγ n ) -1 is nite. Moreover, assumptions in (A3) are technical conditions imposed in order to ensure the reliability of proofs. Those conditions on the density of the couple (X, Y ) were applied in the non-recursive framework for the estimation of the regression function [START_REF] Nadaraya | On estimating regression[END_REF]; [START_REF] Watson | Smooth regression analysis[END_REF] and in the recursive framework Mokkadem et al. (2009a); [START_REF] Slaoui | Large and moderate deviation principles for averaged stochastic approximation method for the estimation of a regression function[END_REF].

In this section, we need to recall the following proposition which introduces the bias and the variance of f n . The proof of this result was depicted in [START_REF] Mokkadem | The stochastic approximation method for the estimation of a multivariate probability density[END_REF].

Bias and variance of f n : Proposition 1. nder ssumptions (A 1 ) -(A 3 )D nd ssuming thtD for ll i, j ∈ {1, . . . , d}D f

(2)

X ij is ontinuous t xD we otin IF sf a ∈ 0, α d+4 D then E[f n (x)] -f X (x) = 1 2(1 -2aξ) d j=1 µ j (K)f X (2) jj (x) h 2 n + o h 2 n . (6) sf a ∈ α d+4 , 1 D then E[f n (x)] -f X (x) = o γ n h -d n . ( 7
) PF sf a ∈ 0, α d+4 D then V ar[f n (x)] = o h 4 n . (8) sf a ∈ α d+4 , 1 D then V ar[f n (x)] = 1 2 -(α -ad)ξ f X (x)R(K) γ n h d n + o γ n h -d n . ( 9
)
3 Main results

In order to explore the asymptotic properties of our estimator π n , we need rst to introduce the following proposition which provides the bias and the variance of a n .

3.1

Bias and variance of a n : Proposition 2. nder ssumptions

(A 1 ) -(A 3 )D nd ssuming thtD for ll i, j ∈ {1, . . . , d}D a (2) ij is ontinuous t xD we otin IF sf a ∈ 0, α d+4 D then E[a n (x, y)] -a(x, y) = 1 2(1 -2aξ) d j=1 µ j (K)a (2) jj (x, y) h 2 n + o h 2 n . (10) sf a ∈ α d+4 , 1 D then E[a n (x, y)] -a(x, y) = o γ n h -d n . ( 11
) PF sf a ∈ 0, α d+4 D then V ar[a n (x, y)] = o h 4 n . ( 12
) sf a ∈ α d+4 , 1 D then V ar[a n (x, y)] = 1 2 -(α -ad)ξ a(x, y)R(K) γ n h d n + o γ n h -d n . ( 13
)
roofF We have

E[a n (x, y)] -a(x, y) = Π n n k=1 Π -1 k γ k (E[Z k (x, y)] -a(x, y)) + Π n [a 0 (x, y) -a(x, y)].
Relying upon the assumption (A 1 ), we have

R d K(z)dz = 1. Hence, it follows that E[Z k (x, y)] -a(x, y) = R d h -d k K x -t h k E[χ k (y)|X = t]f X (t)dt - R d K (t) a(x, y)dt = R d K (z) [a(x -zh k , y) -a(x, y)] dz.
Moreover, Taylor's expansion with integral remainder ensures that

E[Z k (x, y)] -a(x, y) = R d K(z) d i=1 ∂a ∂x i (x, y)z i h k + 1 0 (1 -t) d i,j=1 ∂ 2 a ∂x i ∂x j (x -tzh k , y)z i z j h 2 k dt dz = h 2 k 2 d j=1 µ j (K)a (2) jj (x, y) + h 2 k η k (x).
where,

η k (x) := d i,j=1 R d 1 0 (1 -t) a (2) ij (x -tzh k , y) -a (2) ij (x, y) z i z j K(z)dtdz. Hence, E[a n (x, y)] -a(x, y) = 1 2 d j=1 µ j (K)a (2) jj (x, y)Π n n k=1 Π -1 k γ k h 2 k + Π n n k=1 Π -1 k γ k h 2 k η k (x) + Π n [a 0 (x, y) -a(x, y)].
• For the case a α/(d + 4), we have lim n→+∞ (nγ n ) > 2a and then 1 -2aξ > 0. The application of lemma 1 enables us to write

E[a n (x, y)] -a(x, y) = 1 2(1 -2aξ) d j=1 µ j (K)a (2) jj (x, y) h 2 n + o h 2 n + o (1) + O (Π n ) .
We thus obtain the following desired result

E[a n (x, y)] -a(x, y) = 1 2(1 -2aξ) d j=1 µ j (K)a (2) jj (x, y) h 2 n + o h 2 n .
• For the case a > α/(d+4), we have lim

n→+∞ (nγ n ) > α-a 2 which yields that h 2 n = o γ n h -d n .
Then, the use of lemma 1 leads to

E[a n (x, y)] -a(x, y) = o γ n h -d n .
Therefore, the claimed result ( 11) is established.

For the variance, and owing to the independence of X i , for i = 1, . . . , n, it's obvious that

V ar[a n (x, y)] = Π 2 n n k=1 Π -2 k γ 2 k V ar[Z k (x, y)] = Π 2 n n k=1 Π -2 k γ 2 k E[Z 2 k (x, y)] -E[Z k (x, y)] 2 = Π 2 n n k=1 Π -2 k γ 2 k R d h -d k K 2 (z) a(x -zh k , y)dz - R d K (z) a(x -zh k , y)dz 2 .
As matter of fact, the Taylor's expansions theorem ensures that

V ar[a n (x, y)] = Π 2 n n k=1 Π -2 k γ 2 k h -d k a(x, y) R d K 2 (z) dz + ν k (x) -h d k η k (x) ,
where,

ν k (x) = R d K 2 (z) [a(x -zh k , y) -a(x, y)]dz and η k (x) = R d K (z) a(x -zh k , y)dz 2 .
• For the case a α/(d + 4), we have lim

n→+∞ (nγ n ) > α-ad 2 and then 1 -2aξ > 0. Since we have lim k→+∞ ν k (x) = 0 and lim k→+∞ h k η k (x) = 0, then the application of lemma 1 ensures that V ar[a n (x, y)] = Π 2 n n k=1 Π -2 k γ 2 k h -d k a(x, y)R(K) + ν k (x) -h d k η k (x) = 1 2 -(α -ad)ξ γ n h d n [a(x, y)R(K) + o (1)].
Thus, this leads to the result displayed in (13).

• For the case a < α/(d + 4), we have lim

n→+∞ (nγ n ) > 2a which provides that γ n h -d n = o (h 4 n ) . By applying lemma 1, we infer that V ar[a n (x, y)] = Π 2 n n k=1 Π -2 k γ k o h 4 k = o h 4 n .
In the following theorem, we introduce our main result which provides the bias and the variance of our CCDF multivariate estimator π n .

3.2

Bias and variance of π n : Theorem 1. vet ssumptions (A 1 ) -(A 3 ) hold nd ssume thtD for ll i, j ∈ {1, . . . , d}D

a (2) ij nd f (2) X ij re ontinuous t xD we otin IF sf a ∈ 0, α d+4 D then E[π n (y|x)] -π(y|x) = 1 2(1 -2aξ) 1 f X (x) d j=1 µ j (K) π (2) jj (y|x)f X (x) (14)
+ 2π

(1) j (y|x)f

(1)

X j (x) h 2 n + o h 2 n . sf a ∈ α d+4 , 1 D then E[π n (y|x)] -π(y|x) = o γ n h -d n . ( 15
) PF sf a ∈ 0, α d+4 D then V ar[π n (y|x)] = o h 4 n . ( 16
) sf a ∈ α d+4 , 1 D then V ar[π n (y|x)] = R(K) 2 -(α -ad)ξ π(y|x)(1 -π(y|x)) f X (x) γ n h d n + o γ n h d n . ( 17
)
It is noteworthy that the bias and the variance of the estimator π n dened by the stochastic approximation algorithm (1) mainly depend on the choice of the stepsize (γ n ).

roofF Our proof rests upon the following decomposition, for

f n = 0, n ≥ 0 π n (y|x) -π(y|x) = A n (x, y) f X (x) f n (x) , (18)
with

A n (x, y) = 1 f X (x) (a n (x, y) -a(x, y)) - π(y|x) f X (x) (f n (x) -f X (x)) .
It follows from ( 18) that the asymptotic behavior of π n (y|x) -π(y|x) can be deduced from the one of A n (x, y).

For the bias of π n , we can state

E[A n (x, y)] = 1 f X (x) (E[a n (x, y)] -a(x, y)) - π(y|x) f X (x) (E[f n (x)] -f X (x)) .
Now, using the rst bias part of proposition 1 and proposition 2 and considering the fact that a(x, y) = π(y|x)f X (x), then by combining the assertions ( 6), ( 7),( 10) and ( 11); we obtain the relations ( 14) and ( 15). In order to conrm the variance of π n statement, we have

V ar[A n (x, y)] = 1 (f X (x)) 2 V ar[a n (x, y)] + (π(y|x)) 2 (f X (x)) 2 V ar[f n (x)] -2 π(y|x) (f X (x)) 2 Cov(a n (x, y), f n (x)). ( 19
)
Given that the X k 's are independent, then for all i = 1, . . . , n, k = 1, . . . , n and i = k, we have Cov(Z k (x, y), W i (x)) = 0. Using Taylor's expansion with integral remainder and lemma 1, classical computations entail

Cov(a n (x, y), f n (x)) = R(K) 2 -(α -ad)ξ π(y|x)f X (x) γ n h d n + o γ n h d n . (20)
Consequently, relations ( 16) and ( 17) follow from the combination of assertions ( 8), ( 9), ( 12), ( 13) and (20).

For the case a ≥ α/(d + 4), we deduce with (19) that

V ar[π n (y|x)] = R(K) 2 -(α -ad)ξ π(y|x)(1 -π(y|x)) f X (x) γ n h d n + o γ n h d n .
Proceeding with the same reasoning applied for the case a < α/(d + 4), we obtain the desired result ( 16).

In the sequel, let us present the following theorem which identies the asymptotic normality of our recursive estimator π n . Throughout this paper, we shall denote convergence in probability, convergence in distribution and the Gaussian distribution by

3.3

Weak pointwise convergence rate of π n Theorem 2. vet ssumptions

(A 1 ) -(A 3 ) holdF IF sf there exists nonEnegtive rel c suh tht γ -1 n h d+4 n -→ n→+∞ cD then γ -1 n h d n (π n (y|x) -π(y|x)) D -→ n→+∞ N √ c m(x, y) , σ 2 (x, y) . (21) with m(x, y) = 1 2(1 -2aξ) 1 f X (x) d j=1 µ j (K) π (2) jj (y|x)f X (x) + 2π (1) j (y|x)f (1) X j (x) nd σ 2 (x, y) = R(K) 2 -(α -ad)ξ π(y|x)(1 -π(y|x)) f X (x) . PF sf γ -1 n h d+4 n -→ n→+∞ ∞D then 1 h 2 n (π n (y|x) -π(y|x)) P -→ n→+∞ m(x, y).
roofF We can write, for all n ≥ 0,

x ∈ R d , y ∈ R q , A n (x, y) -E[A n (x, y)] = 1 f X (x) [a n (x, y) -E[a n (x, y)]] - π(y|x) f X (x) [f n (x) -E[f n (x)]] = 1 f X (x) Π n n k=1 S k (x, y), (22)
where

S k (x, y) := Π -1 k γ k (T k (x, y) -E[T k (x, y)]) and T k (x, y) = Z k (x, y) -π(y|x)W k (x).
This proof falls naturally into the application of Lyapunov's theorem for S k (x, y).

On the one hand, we can write

u 2 n := n k=1 V ar[S k (x, y)] = n k=1 Π -2 k γ 2 k V ar [Z k (x, y) -π(y|x)W k (x)] -2π(y|x) Cov (Z k (x, y), W k (x)) .
Therefore, by applying lemma 1, it can be deduced that

u 2 n = n k=1 Π -2 k γ 2 k h -d k R(K)f X (x)π(y|x) (1 -π(y|x)) + o(1) = f X (x) 2 Π 2 n γ n h d n [σ 2 (x, y) + o (1)]. (23)
On the other hand, we can write

E[|T k (x, y)| 2+p ] = O 1 h d(1+p) k , ∀p > 0.
This yields,

n k=1 E[|S k (x, y)| 2+p ] = O n k=1 Π -(2+p) k γ 2+p k E |T k (x, y)| 2+p = O n k=1 Π -(2+p) k γ 2+p k 1 h d(1+p) k
.

For the application of lemma 1, let us assume that there exists a positive real p such that

lim n→+∞ (nγ n ) > 1 + p 2 + p (α -ad).
Then, we obtain

n k=1 E[|S k (x, y)| 2+p ] = O γ 1+p n Π 2+p n h d(1+p) n . It follows that 1 u 2+p n n k=1 E[|S k (x, y)| 2+p ] = O γ 1+p n u 2+p n Π 2+p n h d(1+p) n
.

As a sequel, with the assertion (23), we can write

1 u 2+p n n k=1 E[|S k (x, y)| 2+p ] = O γ n h d n p/2 = o (1) .
Moreover, since we have

lim n→+∞ 1 u 2+p n n k=1 E |S k (x, y) -E[S k (x, y)]| 2+p = lim n→+∞ 1 u 2+p n n k=1 E[|S k (x, y)| 2+p ] = 0,
then, by applying the Lyapunov theorem, we get

1 u n n k=1 (S k (x, y) -E[S k (x, y)]) D -→ n→+∞ N (0 , 1) . This implies 1 u n n k=1 S k (x, y) D -→ n→+∞ N (0 , 1) .
Additionally, the relations ( 18) and ( 22) ensure that

1 u n Π n f X (x) (π n (y|x) -E[π n (y|x)]) D -→ n→+∞ N (0 , 1) . ( 24
)
Given that

u 2 n = f X (x) 2 Π 2 n γ n h d n [σ 2 (x, y) + o (1)],
and by replacing u n with its value in relation ( 24), we deduce that

γ -1 n h d n (π n (y|x) -E[π n (y|x)]) D -→ n→+∞ N 0 , σ 2 (x, y) . ( 25
)
Since we have γ -1 n h d+4 n -→ n→+∞ √ c, then the convergence (21) follows from the combination of relations ( 14), ( 15) and convergence (25).

In order to assess the asymptotic quality of the CCDF recursive estimator π n , we set up the Mean Weighted Integrated Squared Error (M W ISE).

Asymptotic expressions of M W ISE[π n ]

We rst introduce the M W ISE expression:

M W ISE[π n ] = R d+q (E[π n (y|x)] -π(y|x)) 2 + V ar[π n (y|x)] f X 2 (x)f X,Y (x, y)dxdy. (26) Proposition 3. he M W ISE of the estimtor π n is mintined s followsF sf a ∈ 0, α d+4 D then M W ISE[π n ] = 1 4 I 1 (1 -2aξ) 2 h 4 n + o(h 4 n ). sf a = α d+4 D then M W ISE[π n ] = I 2 2 -(α -ad)ξ R(K)γ n h -d n + 1 4 I 1 (1 -2aξ) 2 h 4 n + o(h 4 n ). sf a ∈ α d+4 , 1 D then M W ISE[π n ] = I 2 2 -(α -ad)ξ R(K)γ n h -d n + o γ n h -d n .
Based on the relation ( 26) and by distinguishing the dierent possible cases according to the expressions of the Bias (( 14) and ( 15)) as well as the Variance (( 16) and ( 17)), one can prove this proposition and nd the required result.

The following corollary ensures that the bandwidth which minimizes the M W ISE of π n depends on the choice of the stepsize (γ n ). As a matter of fact, the corresponding M W ISE depends also on (γ n ).

Corollary 1. vet ssumptions

(A 1 ) -(A 3 ) holdF o minimize the M W ISE of π n , the ndwidth (h n ) must e equl to   d (1 -2aξ) 2 2 -(α -ad)ξ I 2 I 1 R(K) 1 d+4 γ 1 d+4 n   . reneD the orresponding M W ISE is determined y M W ISE[π n ] = d + 4 4d d d+4 I 1 (1 -2aξ) 2 d d+4 I 2 2 -(α -ad)ξ 4 d+4 R(K) 4 d+4 γ 4 d+4 n + o γ 4 d+4 n .
The following corollary holds in the special case where (γ n ) is chosen as

(γ n ) = (γ 0 n -1 ) in order to minimize the M W ISE[π n ]. Corollary 2. vet ssumptions (A 1 ) -(A 3 ) holdF o minimize the M W ISE of π n D we need to opt for the stepsize (γ n ) in GS(-1) suh tht lim n→∞ (nγ n ) = γ 0 F hen the ndwidth (h n ) must e equl to d (γ 0 (d + 4) -2) 2(d + 4) I 2 I 1 R(K) 1 d+4 n -1 d+4 . gonsequentlyD the orresponding M W ISE is identi(ed y M W ISE[π n ] = (d + 4) 3d+8 d+4 4 d + 6 d+4 d d d+4 γ 2 0 (γ 0 (d + 4) -2) -(2d+4) d+4 I d d+4 1 I 4 d+4 2 R(K) 4 d+4 n -4 d+4 + o n -4 d+4 .
In order to get the optimal choice of (γ n ), we deduce that the minimum of M W ISE[π n ] is achieved at γ 0 = 1. Hence, we introduce the following corollary.

Corollary 3. vet ssumptions

(A 1 ) -(A 3 ) holdF o minimize the M W ISE of π n D we must selet the stepsize (γ n ) in GS(-1) suh tht lim n→∞ (nγ n ) = 1F hereforeD the optiml ndwidth (h n ) must equl d(d + 2) 2(d + 4) I 2 I 1 R(K) 1 d+4 n -1 d+4 . ( 27
)
es resultD the orresponding M W ISE is expressed y

M W ISE[π n ] = (d + 4) 3d+8 d+4 4 d+6 d+4 d d d+4 (d + 2) 2d+4 d+4 I d d+4 1 I 4 d+4 2 R(K) 4 d+4 n -4 d+4 + o n -4 d+4 .
Remark 1. xote thtD for the prtiulr se where the stepsize

(γ n ) is in GS(-1) suh tht lim n→∞ (nγ n ) = 1 nd the ndwidth (h n ) is hosen suh tht lim n→∞ nh d+4 n
= 0 @whih orreE sponds to undersmoothingAD the symptoti normlity of our proposed estimtor is indited s follows

nh d n (π n (y|x) -π(y|x)) D -→ n→+∞ N 0 , d + 4 2(d + 2) R(K) π(y|x)(1 -π(y|x)) f X (x) . ( 28
)
The statistical inference of the CCDF multivariate non-recursive estimator π n is addressed in our next section.

The following results can be handled in nearly the same way as π n . The unique dierence lies in the fact that it pertains to a non-recursive case. (See [START_REF] Hall | Methods for estimating a conditional distribution function[END_REF] for more details of the univariate case.)

Asymptotic properties of π n

In order to tackle the asymptotic properties of our estimator π n , we need rst to introduce the following proposition which provides the bias and the variance of π n .

4.1

Bias and variance of π n Proposition 4. vet ssumptions (A 1 ) nd (A 3 ) holdF hen the is nd vrine of xdryE tson9s estimtor re displyed s followsF

IF he is of π n X E[ π n (y|x)]-π(y|x) = 1 2f X (x) d j=1 µ j (K) π (2) jj (y|x)f X (x) + 2π (1) j (y|x)f X (1) j (x) h 2 n +o h 2 n . PF he vrine of π n X V ar[ π n (x)] = R(K) π(y|x)(1 -π(y|x)) f X (x) 1 nh d n + o 1 nh d n .
The following proposition yields the distribution convergence rate of the non-recursive estimator.

4.2

Asymptotic normality of π n Theorem 3. vet ssumptions

(A 1 ) -(A 3 ) hold nd suppose tht nh d+4 n -→ n→+∞ 0F henD nh d n ( π n (y|x) -π(y|x)) D -→ n→+∞ N 0 , R(K) π(y|x)(1 -π(y|x)) f X (x) . ( 29
)
Remark 2. st is ovious to infer from the expressions (28) nd (29) tht our gghp proposed estimtor is etter thn nonEreursive one in terms of vrineF

In the next subsection, we exhibit the expression of the Mean Weighted Integrated Squared Error of Nadaraya-Watson's estimator.

M W ISE[ π n ] = 1 4 I 1 h 4 n + I 2 R(K) 1 nh d n + o h 4 n + 1 nh d n .
Proposition 5. vet ssumptions

(A 1 ) nd (A 3 ) holdF o minimize the M W ISE of π n D the ndwidth (h n ) must e equl to d I 2 I 1 R(K) 1 d+4 n -1 d+4 . ( 30
)
reneD the orresponding M W ISE is determined y

M W ISE[ π n ] = d + 4 4d d d+4 I 4 d+4 2 I d d+4 1 R(K) 4 d+4 n -4 d+4 + o n -4 d+4 .

Bandwidth selection

Although theoretical asymptotic study yields the optimal bandwidth, the fact that we do not know the density function makes it hard to interpret it in practice. Hence, kernel smoothing in non-parametric statistics requires the choice of a bandwidth parameter. This choice is crucial to obtain a good rate of consistency of the kernel estimators. It has a signicant inuence on the size of the bias. One has to nd an appropriate bandwidth that produces an estimator which has a good balance between the bias and the variance of the estimator of the function a(•, •) as well as f (•). It is worth noticing that the bandwidth selection methods reported in the literature can be divided into three broad classes: the cross-validation techniques, the plug-in ideas and the bootstrap procedure. In this investigation, we are basically interested in the plug-in method. [START_REF] Altman | Bandwidth selection for kernel distribution function estimation[END_REF] developed an ecient method of bandwidth selection, which minimizes an estimate of the mean weighted integrated squared error, using the density function as a weight function. For this reason, we followed the work of Slaoui (2014a) see appendix 5.1 for additional details.

Plug-in bandwidth selection

As a result to the plug-in procedure, based on the expression of the M W ISE, we estimate the unknown quantities I 1 and I 2 by elaborating asymptotic unbiased estimators. This process is known as a plug-in estimate. Basically, we introduce (b n ) ∈ GS(-δ), δ ∈ (0, 1) . In practice, [START_REF] Altman | Bandwidth selection for kernel distribution function estimation[END_REF]

set b n = n -δ min s, Q 3 -Q 1 1.349
, with s being the sample standard deviation and Q 1 , Q 3 being the rst and third quartiles.

In the following and for the sake of simplicity, the kernel K we shall use is considered as a product of univariate kernels K satisfying R K(x)dx = 1.

In addition, we note:

I 1 = µ 2 (K) (J 1 -2J 2 + J 3 ) ,
where

J 1 = R d+q d j=1 a (2) jj (x, y) 2 f X,Y (x, y)dxdy, J 3 = R d+q d j=1 f X (2) jj (x) 2 π 2 (y|x)f X,Y (x, y)dxdy, J 2 = R d+q d j=1 a (2) jj (x, y) d j=1 f X (2) jj (x) π(y|x)f X,Y (x, y)dxdy, µ(K) = R z 2 K(z)dz.

5.2

Recursive estimator π n :

To estimate the optimal bandwidth (27), we need to estimate I 1 and I 2 . Here we can write

a n (x, y) = Π n n k=1 Π -1 k γ k h -d k χ k (y) K x -X k h k = Π n n k=1 Π -1 k γ k h -d k d i=1 K x -X k i h k χ k i (y)
and

f n (x) = Π n n k=1 Π -1 k γ k h -d k K x -X k h k = Π n n k=1 Π -1 k γ k h -d k d i=1 K x i -X k i h k .
Estimation of I 1 :

J 1 = Π 2 n n n i,j,k=1 i =j =k Π -1 j Π -1 k γ j γ k b -(d+2) j b -(d+2) k    d v=1 K (2) b X iv -X jv b j d l=1 l =v K b X i l -X j l b j    ×    d v=1 K (2) b X iv -X kv b k d l=1 l =v K b X i l -X k l b k    q s=1 χ (j+1) s (Y is )χ (k+1) s (Y is ), J 2 = Π 2 n n n i,j,k=1 i =j =k Π -1 j Π -1 k γ j γ k b -(d+2) j b -(d+2) k    d v=1 K (2) b X iv -X jv b j d l=1 l =v K b X i l -X j l b j    ×    d v=1 K (2) b X iv -X kv b k d l=1 l =v K b X i l -X k l b k    q s=1 χ (i+1) s (Y is )χ (j+1) s (Y is ), J 3 = Π 2 n n n i,j,k,m=1 i =j =k =m Π -1 j Π -1 k γ j γ k b -(d+2) j b -(d+2) k    d v=1 K (2) b X iv -X jv b j d l=1 l =v K b X i l -X j l b j    ×    d v=1 K (2) b X iv -X kv b k d l=1 l =v K b X i l -X k l b k    q s=1 χ (i+1) s (Y is )χ (m+1) s (Y is ),
where, K b stands for a kernel with bandwidth b n such that δ = -2 5 , and K

(2) b corresponds to the second derivative of a kernel K b with the associated bandwidth b n such that δ = -3 14 . Note that our choice of the parameter δ is based on the work of Slaoui (2014a). At this stage, we obtain

I 1 = µ 2 (K) J 1 -2 J 2 + J 3 .
Estimation of I 2 :

I 2 = Π n n n i,k=1 i =k Π -1 k γ k b -1 k d l=1 K b X i l -X k l b k q s=1 χ (i+1) s (Y is ) 1 -χ (k+1) s (Y is ) ,
where, K b is a kernel with bandwidth b n such that δ = -2 5 .

As a result, the plug-in estimator of ( 27) is determined by

(h n ) =   d(d + 2) 2(d + 4) 1 d+4 I 2 I 1 1 d+4 R(K) 1 d+4 n -1 d+4   , (31)
Eventually, an estimator of M W ISE[π n ] is specied by

M W ISE[π n ] = (d + 4) 3d+8 d+4 4 d+6 d+4 d d d+4 (d + 2) d+6 d+4 I 1 d d+4 I 2 4 d+4 R(K) 1 d+4 n -4 d+4 + o n -4 d+4 .

5.3

Non-Recursive estimator π n :

To estimate the optimal bandwidth (30), we need to estimate I 1 and I 2 .Therefore, we can state

a n (x, y) = 1 nh d n n k=1 χ k (y) K x -X k h n = 1 nh d n n k=1 d i=1 K x -X k i h k χ k i (y)
and

f n (x) = 1 nh d n n k=1 K x -X k h n = 1 nh d n n k=1 d i=1 K x i -X k i h k .
Estimation of I 1 :

J 1 = 1 n 3 b 6 n n i,j,k=1 i =j =k    d v=1 K (2) b X iv -X jv b n d l=1 l =v K b X i l -X j l b n    ×    d v=1 K (2) b X iv -X kv b n d l=1 l =v K b X i l -X k l b n    q s=1 χ (j+1) s (Y is )χ (k+1) s (Y is ), J 2 = 1 n 3 b 6 n n i,j,k=1 i =j =k    d v=1 K (2) b X iv -X jv b n d l=1 l =v K b X i l -X j l b n    ×    d v=1 K (2) b X iv -X kv b n d l=1 l =v K b X i l -X k l b n    q s=1 χ (i+1) s (Y is )χ (j+1) s (Y is ), J 3 = 1 n 4 b 6 n n i,j,k,m=1 i =j =k =m    d v=1 K (2) b X iv -X jv b n d l=1 l =v K b X i l -X j l b n    ×    d v=1 K (2) b X iv -X kv b n d l=1 l =v K b X i l -X k l b n    q s=1 χ (i+1) s (Y is )χ (m+1) s (Y is ),
where, K b is a kernel with bandwidth b n such that δ = -2 5 , and K

(2) b

is the second derivative of a kernel K b with the associated bandwidth b n such that δ = -3 14 . Then, we obtain

I 1 = µ 2 (K) J 1 -2 J 2 + J 3 .
Estimation of I 2 :

I 2 = 1 n 2 b n n i,k=1 i =k d l=1 K b X i l -X k l b n q s=1 χ (i+1) s (Y is ) 1 -χ (k+1) s (Y is ) ,
where, K b is a kernel with bandwidth b n such that δ = -2 5 .

As a result, the plug-in estimator of ( 30) is denoted by

(h n ) =   I 2 I 1 1 d+4 R(K) 1 d+4 n -1 d+4   , (32)
Finally, a non-recursive estimator of M W ISE[π n ] is provided by

M W ISE[ π n ] = 5 4 I 2 4 d+4 I 1 1 d+4 R(K) 1 d+4 n -4 d+4 + o n -4 d+4 .
The major aim of our next section lies in comparing the performance of our recursive estimator (1) with that of non-recursive Nadaraya-Watson one (2).

Numerical applications

Let's start our numerical studies with some simulations with dierent dimensions Models.

Simulation studies

In order to compare the proposed recursive estimator with the Nadaraya-Watson non-recursive one, we consider three sample sizes: n = 100, 200 and 500, a xed number of simulations: N = 500 and four distribution models:

• Model 1: (X, Y ) ∈ R × R: Y = 2 sin(πX) +
, where X follows the binomial distribution B (2, 1/3) and follows the normal distribution N (0, 1) .

• Model 2: (X, Y ) ∈ R 2 × R: Y = exp(-X/2) +
, where X follows the exponential distribution P (1/2) and follows the normal distribution N (0, 1/2) .

• Model 3:

(X, Y ) ∈ R 3 × R 2 : Y = AX + with A = 1 1 1 1 -2 1 , X = 0 × 1 Z<=0.5 + 1 Z>0.5 where Z follows the uniform distribution U 0 0 0 , 1 1 1
and follows the normal distribution N (0, 1/2) .

• Model 4: (X, Y ) ∈ R 10 × R 10 : Y = exp(X) + , with X = Z where Z follows the 10-dimensional normal distribution N (0 10 , I 10 ) and follows the normal distribution N (0, 1) .

We denote by π * i the reference CCDF and by π i the test CCDF. Then, we calculate the following two measures:

• Mean squared error: M SE = 1 n n i=1 (π i -π * i ) 2 .
• The linear correlation:

Cor = Cov(π i , π * i ) σ(π i )σ(π * i )
.

In what follows, we portray the dierent steps of the simulation algorithm in the multivariate case.

Simulation Algorithm

Algorithm 1 K is the Gaussian kernel, d the dimension size, n the simple size, N p the number of observations and N the number of iterations.

Input: K, d, n, N p and N .

1: A random initialization of Π (0) (resp. Π (0) .)

2: for l = 1, . . . , N do 3:

A random sample vectors X 1 , . . . , X d and Y of length n.

4:

A choice value for the recursive bandwidth vectors h 1 , . . . , h n . (resp. the non-recursive bandwidth values h n ) using the plug-in approach given in (31) (resp. ( 32)).

5:

The choice of the stepsize (γ n ) = (n -1 ).

6:

We x x 1 , . . . , x d and consider an arbitrary sampling vector T of y of length N p .

7:

π l (y|x) = n k=1 kγ k 1 {Y k y} h -d k d i=1 K x i -X k i h k n k=1 kγ k h -d k d i=1 K x i -X k i h k
for the multivariate recursive CCDF estimator. (resp. π l (y|x) = n k=1

1 {Y k y} d i=1 K x i -X k i h k n k=1 d i=1 K x i -X k i h k
for the multivariate nonrecursive CCDF estimator). Π (l) = π l (T |x). (resp. Π (l) = π l (T |x).)

8: end for l) .)

9: π = N -1 N l=1 Π (l) (resp. π = N -1 N l=1 Π (
output: The vectors π and π. Table 3: Quantitative comparison between the recursive estimator and the non-recursive one with stepsize (γ n ) = (n -1 ) through a plug-in method for Model 3.

x = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0) Departing from Tables 1, 2, 3 and 4, we conclude that:

1. The MSE of the proposed recursive estimator with stepsize (γ n ) = (n -1 ) through a plug-in method is smaller than that of Nadaraya-Watson's non-recursive estimator.

2. The estimators get closer to the true CCDF function as sample size increases, i.e., the M SE decreases as the simple size increases and therefore the Cor increases as the sample size increases.

6.3

Real Datasets:

In this section, our focal point is to examine two real datasets Models, namely the Insurance Company Benchmark (COIL 2000) dataset as well as the French Hospital Data of COVID-19.

6.3.1 Application 1: The Insurance Company Benchmark (COIL 2000) dataset

The (COIL 2000) dataset is found in data.world website https://data.world/uci/insurance-company-benchmark-coil-2000). As far as our application is concerned, we shall consider the following two models:

• Model 1: X corresponds to the sociodemographic data attribute number 16 and Y stands for the whole 5822 observations of customer records.

• Model 2: X corresponds to the sociodemographic 5-dimensional data attributes number 6,8,11,12 and 13 and Y stands for the whole 5822 observations of customer records.

x=0 x=1 (0, 0, 0, 0, 0).

Application 2: French Hospital Data of COVID19

The French Hospital data of the COVID-19 epidemic are extracted from https://www.data.gouv.fr/fr/datasets/donnees-hospitalieres-relatives-a-lepidemie-de-covid-19/.

The Santé publique France's mission is to improve and protect the health of populations.

During the health crisis linked to the COVID-19 epidemic, Santé publique France has taken in charge monitoring and understanding the dynamics of the epidemic, anticipating the dierent scenarios and implementing actions to prevent and limit the transmission of this virus on the national territory.

Description of the dataset

This dataset provides information on the hospital situation regarding the COVID-19 epidemic. We have opted for the rst proposed le: Hospital data related to the COVID-19 epidemic by department (dep), sex of the patient (sex), number of hospitalized patients (hosp), number of persons currently in intensive care or resuscitation (rea), number of persons currently in follow-up and rehabilitation care (SSR) or long-term care units (USLD), number of persons currently in conventional hospitalization (HospConv), number of persons currently hospitalized in another type of service (autres), cumulative number of persons returning home (rad) or cumulative number of dead persons (dc).

The data have been daily updated. For the current application, we considered the data of 28/07/2021, with a total of 150894 observations. For simplicity reason, we have chosen to just study the department of 'Vienne' database. Therefore, for our application, we served of a dataframe of 1494 observations and 6 variables. Hence, we shall consider the following three models:

• Model 1: X = dc and Y = hosp.

• Model 2: X1 = sex, X2 = rea, X3 = dc and Y = hosp.

• Model 3: X1 = sex, X2 = rea, X3 = dc, Y 1 = hosp and Y 2 = rad.

x=119 x=17

Data interpretation:

Referring to Tables 5, 6 and Figures 11, 12, 13 and 14 , we conclude that:

1. For all considered Models, the proposed recursive estimator with stepsize (γ n ) = (n -1 ) through a plug-in method outperformed the non-recursive one in terms of estimation error M SE and Cor.

2. The proposed recursive estimator is closer to the true CCDF function, compared with Nadaraya-Watson's non-recursive estimator.

Concerning the COVID-19 epidemic Model 1, we can infer that, for a xed number of deaths x = 17, the proportion of hospitalized cases less then 20 is 50% and the proportion of hospitalized cases less then 50 is 85%. Moreover, 99% of the population have a number of hospitalized cases less than 100.

Likewise, the COVID-19 epidemic Model 2 yielded the same results as model 1. Indeed, we recorded for the women gender a xed value of sex x 1 = 2, a xed number of REA persons x 2 = 0 and deaths x 3 = 17.

Conclusion and Perspectives

In this work, we elaborated a multivariate recursive CCDF estimator. We tackled the asymptotic properties of the proposed estimator by providing the bias as well as the variance in order to demonstrate that our estimator asymptotically follows a normal distribution. Subsequently, we revealed that the use of our recursive estimator with an appropriate choice of the bandwidth and the stepsize enables us to get closer to the true conditional cumulative distribution function rather than non-recursive one. The basic merit of recursive estimators resides in the fact that one can update the estimation with each additional new observation. Therefore, instead of re-running the data each time, it is possible to rewrite our considered estimator as a combination of two (or more) estimators, where each estimator is based on separate datasets, which can be very interesting to keep the computational cost reasonably low. It is noteworthy that all computation and simulation have been done using the R statistical software.

A future research direction would be to extend our ndings to the setting of serially dependent observations, α-mixing framework like in [START_REF] Khardani | Recursive kernel density estimation and optimal bandwidth selection under alpha-mixing data[END_REF]. Another outstanding direction lies in recursively estimating the modal regression, which requires non trivial mathematics. This would go well beyond the scope of the present paper.

  M W ISE[ π n ] Corollary 4. he M W ISE expression of the nonEreursive gghp estimtor is given y

Figure 1 :

 1 Figure 1: The reference CCDF for Model 1 for one simple simulation.

Figure 2 :

 2 Figure 2: The recursive CCDF estimator for Model 1 for one simple simulation.

Figure 3 :

 3 Figure 3: The non-recursive CCDF estimator for Model 1 for one simple simulation.

Figure 4 :

 4 Figure 4: Qualitative comparison between the recursive estimator and the non-recursive one for Model 1 with n = 200, N = 500 and x = 0.

Figure 5 :

 5 Figure 5: Qualitative comparison between the recursive estimator and the non-recursive one for Model 1 with n = 500, N = 500 and x = 0.

Figure 6 :

 6 Figure 6: Qualitative comparison between the recursive estimator and the non-recursive one for Model 2 with n = 100, N = 500 and x = (0, 0).

Figure 7 :

 7 Figure 7: Qualitative comparison between the recursive estimator and the non-recursive one for Model 2 with n = 500, N = 500 and x = (0, 0).

Figure 8 :

 8 Figure 8:The reference CCDF for Model 3 for one simple simulation with n = 500 and x = (1, 1, 1).

Figure 9 :

 9 Figure 9: The recursive CCDF estimator for Model 3 for one simple simulation with n = 500 and x = (1, 1, 1).

Figure 10 :

 10 Figure 10: The non-recursive CCDF estimator for Model 3 for one simple simulation with n = 500 and x = (1, 1, 1).

Figure 11 :

 11 Figure 11: Qualitative comparison between the recursive estimator and Nadara Watson's non-recursive one for the dataset Model 1 with x = 1.

Figure 12 :

 12 Figure 12: Qualitative comparison between the recursive estimator and the non-recursive one for the dataset Model 2 with x =

Figure 13 :

 13 Figure 13: Qualitative comparison between the recursive estimator and the non-recursive one for the COVID-19 epidemic dataset Model 1 with x = 17.

Figure 14 :

 14 Figure 14: Qualitative comparison between the recursive estimator and the non-recursive one for the COVID-19 dataset Model 2 with x = (2, 0, 17).

Table 1 :

 1 Quantitative comparison between the recursive estimator and the non-recursive one with stepsize (γ n ) = (n -1 ) through a plug-in method for Model 1.

	x = 0	x = 2

Table 2 :

 2 Quantitative comparison between the recursive estimator and the non-recursive one with stepsize (γ n ) = (n -1 ) through a plug-in method for Model 2.

	x = (0, 0, 0)	x = (1, 1, 1)

Table 4 :

 4 Quantitative comparison between the recursive estimator and the non-recursive one with stepsize (γ n ) = (n -1 ) through a plug-in method for Model 4.

  Information about customers consists of 86 variables and includes product usage data and sociodemographic data derived from zip area codes. The data are supplied by the Dutch data mining company Sentient Machine Research and rest on a real world business problem. The training set involves over 5000 descriptions of customers, including the information of whether or not they have a caravan insurance policy. A test set includes 4000 customers whom only the organizers know if they have a caravan insurance policy. This corresponds to a Dataset to train and validate prediction models and build up a description (5822 customer records). Each record consists of 86 attributes, incorporating sociodemographic data (attribute 1-43) and product ownership (attributes 44-86). The sociodemographic data are derived from zip codes. All customers living in areas with the same zip code have the same sociodemographic attributes.

Table 5 :

 5 Quantitative comparison between Nadaraya-Watson estimator and the proposed estimator with stepsizes (γ n ) = (n -1 ) through a plug-in method for the Insurance Company Benchmark (COIL 2000) dataset case.

		Nadaraya's estimator	Recursive estimator	Nadaraya's estimator	Recursive estimator
	Model 1 M SE	0.006979541	0.004748803	0.001999743	0.001390595
	Cor	0.987112335	0.990653245	0.995770041	0.996574025
		x=(0,0,0,0,0)	x=(1,1,1,2,2)
		Nadaraya's estimator	Recursive estimator	Nadaraya's estimator	Recursive estimator
	Model 2 M SE	0.000708417	0.0005969373	0.00414241	0.003705582
	Cor	0.996664739	0.9979940979	0.97540614	0.976614363

Table 6 :

 6 Quantitative comparison between Nadaraya-Watson estimator and the proposed estimator with stepsizes (γ n ) = (n -1 ) through a plug-in method for the COVID-19 epidemic dataset case.

		Nadaraya's estimator	Recursive estimator	Nadaraya's estimator	Recursive estimator
	Model 1 M SE	0.02538461	0.01900265	0.008079359	0.007509241
	Cor	0.86514927	0.88894736	0.781752320	0.789593904
		x = (2, 0, 17)	x = (1, 0, 0)
		Nadaraya's estimator	Recursive estimator	Nadaraya's estimator	Recursive estimator
	Model 2 M SE	0.02176777	0.01416798	0.01026319	0.005681196
	Cor	0.55644221	0.63197054	0.74621610	0.819493783
	Model 3 M SE	0.02107916	0.009991914	0.04130757	0.03504694
	Cor	0.85362436	0.920385555	0.61896604	0.64659332


Lisheng Kuang

Maud Marchal

Marco Aggravi

Paolo Robuffo Giordano

Claudio Pacchierotti

Design of a 2-DoF Haptic Device for Motion Guidance

Keywords:

We present a 2-degrees-of-freedom (2-DoF) haptic device, which can be either used as a grounded or a hand-held device. It is composed of two platforms moving with respect to each other, actuated by two servomotors housed in one of structures. The device implements a rigid coupling mechanism between the two platforms, based on a threelegged 3-4R constrained parallel linkage, with the two servomotors actuating two of these legs. The device can apply position/kinesthetic haptic feedback to the user hand(s). This paper presents the device and its kinematics, together with a human subjects experiment where we evaluate its capabilities to provide meaningful directional information.

Introduction

Kinesthetic haptic interfaces have been very popular in the past, for applications ranging from industrial to surgical robotics. Researchers have designed many different types of such interfaces, focusing on improving their, e.g., peak force, bandwidth, workspace, and/or price, according to the target field of application. In this respect, we can identify two main categories of kinesthetic interfaces: grounded and ungrounded. Grounded devices have their base placed on an external support, such as a table, while ungrounded devices have their base on the user's body [START_REF] Pacchierotti | Wearable haptic systems for the fingertip and the hand: taxonomy, review, and perspectives[END_REF]. Grounded kinesthetic devices include popular commercial kinesthetic systems such as the Virtuose (Haption, FR), Omega.x (Force Dimension, CH), Falcon (Novint Tech., USA) and the Phantom (Geomagic, 3D Systems, USA) series. More recently, in research, Jang et al. [START_REF] Jang | Evaluation of grounded isometric interface for whole-body navigation in virtual environments[END_REF] presented a grounded isometric interaction device to induce whole-body interaction; Okui et al. [START_REF] Okui | Delta-type four-dof forcefeedback device composed of pneumatic artificial muscles and magnetorheological clutch and its application to lid opening[END_REF] designed a delta-type 4-degrees-of-freedom (4-DoF) grounded haptic device actuated with a magnetorheological clutch, thus able to adjust its stiffness and viscosity; and Satler et al. [START_REF] Satler | Control of a desktop mobile haptic interface[END_REF] devised a portable interface composed of controlled wheel torques to render forces to a user handle placed on the top of the device. On the other hand, ungrounded kinesthetic devices come in very different forms, spanning from hand-held devices for gaming and VR interaction [START_REF] Kovacs | Haptic pivot: On-demand handhelds in vr[END_REF][START_REF] De Tinguy | Weatavix: wearable actuated tangibles for virtual reality experiences[END_REF] or guidance [START_REF] Spiers | Design and evaluation of shape-changing haptic interfaces for pedestrian navigation assistance[END_REF][START_REF] Walker | Holdable haptic device for 4-dof motion guidance[END_REF] to body-worn exoskeletons [START_REF] Chinello | A modular wearable finger interface for cutaneous and kinesthetic interaction: control and evaluation[END_REF][START_REF] Sarac | Design requirements of generic hand exoskeletons and survey of hand exoskeletons for rehabilitation, assistive, or haptic use[END_REF]. This paper presents a 2-degrees-of-freedom (2-DoF) kinesthetic device that can be either attached onto an external support or held between the two hands. Its design is inspired by the anti-parallelogram mechanism [START_REF] Okada | Development of a cybernetic shoulder-a 3-dof mechanism that imitates biological shoulder motion[END_REF] featuring a quaternion joint [START_REF] Kim | Quaternion joint: Dexterous 3-dof joint representing quaternion motion for high-speed safe interaction[END_REF]. It is composed of two structures, connected by three articulated legs which are in turn actuated by two servo motors. By moving with respect to each other, they can provide the user with directional information.

Device design and actuation

The proposed device is shown in Fig. 1. A video is available at https://youtu. be/vc6B-OOj590.

Mechanism and Structure

The design of the proposed device is inspired by the principle of the antiparallelogram mechanism [START_REF] Okada | Development of a cybernetic shoulder-a 3-dof mechanism that imitates biological shoulder motion[END_REF], which consists of three identical supporting linkages forming an interlaced structure with no interference between them [START_REF] Kim | Quaternion joint: Dexterous 3-dof joint representing quaternion motion for high-speed safe interaction[END_REF]. Compared to standard serial mechanisms, parallel mechanisms enable fast dynamics and high payload with relatively small size and low weight. Moreover, a large range of motion and uniform manipulability can be obtained by choosing appropriate dimensional parameters and actuation.

As shown in Fig. 1, the device consists of two platforms, a lower and an upper platforms, connected by three legs each having a spiral curved link. The device has dimensions of 15×15×23 cm and weighs 150 g. The two ends of each link are connected with two serial revolute joints to the platforms, forming a 3-4R coupling parallel mechanism which can move freely with 2 degrees of freedom, according to the orthogonal rolling motion on each leg [START_REF] Kim | Quaternion joint: Dexterous 3-dof joint representing quaternion motion for high-speed safe interaction[END_REF]. If the lower platform is fixed on an external support (as in Fig. 1(a)), the motion of the upper platform is confined on the surface of a sphere centered in the center of the lower platform (see grey dotted lines in Fig. 1(b)). Two Hitech-625MG servo motors are attached to two legs on the lower platform. One revolute joint on the short linkage is mounted to the motor shaft, the other joint connects to the leg, as shown in Fig. 1. In this configuration, the two legs equipped with the servomotors have an active rotation on the lower joints and a passive rotation on the other three joints. By changing the actuation of the two motors, the upper platform can reach any point within its workspace (the surface of a sphere, as mentioned above). The graphical method using reciprocal screw system theory can be used [START_REF] Yu | Mobility and singularity analysis of a class of two degrees of freedom rotational parallel mechanisms using a visual graphic approach[END_REF] to analyze the mobility of the proposed parallel mechanism. We study the kinematics of the device so as to evaluate the relationship between the motion of the servomotors and that of the upper platform with respect to the lower one. As illustrated in Figs. 1 and2, three legs are regarded as individual open-chain manipulators. The connection points between the upper platform and legs 1, 2, 3 are defined as P 1 , P 2 , P 3 , respectively, while the connection points between the lower platform and the legs is denoted by points B 1 , B 2 , B 3 . Let us consider the device in its grounded configuration, with the base attached to an external support. The coordinate frames of the upper (mobile) and lower (static) platforms, C p (O p -x p y p z p ) and C b (O b -x b y b z b ), are fixed on their geometric centers, O p and O b , respectively. The frames are defined as indicated in Figs. 1 and2.

Kinematics Analysis

The plane where the z b and z p axes coexist is defined as the bending plane, which is always perpendicular to the surfaces of the two platforms. Here, the orientation angle θ and bending angle γ denote the angle between z b and z p on the bending plane and the angle from the x b axis to the bending plane, respectively. Considering the bending and orientation angles, the homogeneous transformation from C b to C p is as follows:

b p T = R(Z, γ)R(Y, θ 2 )T (h)R(Y, θ 2 )R(Z, -γ) =     1 -2C 2 (γ)S 2 ( θ 2 ) -S(2γ)S 2 ( θ 2 ) C(γ)S(θ) hC(γ)S( θ 2 ) -S(2γ)S 2 ( θ 2 ) 1 -2S 2 (γ)S 2 ( θ 2 ) S(γ)S(θ) hS(γ)S( θ 2 ) -C(γ)S(θ) -S(γ)S(θ) C(θ) hC( θ 2 ) 0 0 0 1     (1)
where S(x) = sin(x) and C(x) = cos(x), R(., .) denotes a rotation, T an homogeneous transformation around the z axis, and h the distance between O b and O p , which is also the diameter of the sphere onto which the upper platform moves. Geometrically, we can also derive the relationship (forward and inverse kinematics) between γ and θ with respect to the motor's inputs α 1 and α 2 . As shown in Fig. 2b

(α 1 ) = tan( O p EK) = O p K EK = O p K O b K sin(γ + π 3 ) = h cos θ 2 h sin θ 2 sin(γ + π 3 ) = 1 tan θ 2 sin(γ + π 3 ) (2) tan(α 2 ) = tan( O p F K) = O p K F K = 1 tan θ 2 sin(γ -π 3 ) (3)

Experimental evaluation

We carried out an experiment evaluating the capabilities of the device in providing directional/motion information. The device renders a set of shapes with its end-effector, i.e., the moving upper platform, that users are asked to recognize.

Setup

The experimental setup is shown in Fig. 1. The device is used in its grounded configuration, with the lower platform attached to a table. Users are seated in front of the device and are asked to place their dominant hand on the device upper platform. A cardboard prevents the user from seeing his or her hand on the device during the experiment. A computer screen is also placed in front of the user, from which he or she can receive information about the experiment and answer the related questions.

Participants

Fourteen subjects (1 female and 13 males, aged from 23 to 33 years) participated in the experiment. Six are left-handed, eight are right-handed. Participants received an information sheet with the experiment details and signed a consent form. The study has been approved by Inria's ethics committee (Saisine 513).

Procedure

The experiment is divided in two blocks, carried out one after the other.

In the first one, the device moves along eight linear patterns, shown in Figs. 3a and3c. We recall that the center of the upper platform O p moves across the surface of a sphere having diameter h and centered at the center of the lower platform O b (see Sec. 2.2). Each linear pattern starts from the resting position of O p , then moves towards its designated direction and back, along a circular arc with chord of 20 cm (see also the video). For example, the linear pattern referred to as "South" in Fig. 3, starts from the resting position of O p , then moves towards the S (South) direction (see Fig. 3a), then back to the resting position of O p , then moves towards the N (North) direction, and finishes in the resting position of O p . All linear patterns start and finish in the resting position of O p . The user sits on the South side of the device.

In the second block, the device moves along three shape patterns, shown in Figs. 4a and4b. Similarly as before, these shapes are actuated over the surface of the sphere that is the workspace of our device. As indicated in Fig. 4a, they all start from a point located at the North of the resting position of O p (named P in the Figure ). Before starting the experiment, the experimenter explained the procedure to the user and spent about two minutes adjusting the chair armrest. Subjects placed their palm on the upper platform of the device and lean on the chair armrest to keep the forearm at the same level of the upper platform, ensuring maximum comfort. After each pattern was rendered, users were asked to select which pattern, in their opinion, the device just rendered. This choice was made through a GUI on the computer screen in front of the user. The device actuating all the patterns is shown at https://youtu.be/vc6B-OOj590.

Each pattern was provided five times, yielding (8 linear patterns + 3 shape patterns)×5 = 55 repetitions of this pattern recognition task.

Immediately after the experiment, participants were asked to fill in a questionnaire where we asked which pattern was the harder/easier to recognize and if they had any further comment about the rendering and the experiment.

Results

Figs. 3b and 4b report the results of the experiment, in the form of two confusion matrices showing the percentage of recognition of the rendered vs. the recognized patterns. For the linear patterns the chance level is 1/8 (12.5%), while for the shape patterns the chance level is 1/3 (33.3%). The questionnaire showed that "East" and "West" linear patterns were considered as the easiest to identify, while the square shape pattern resulted to be the hardest. Diagonal linear patterns were reported to be harder to identify than other linear patterns.

Discussion and Conclusions

This paper presented a kinesthetic haptic device composed of two platforms connected by three interleaved legs. Two servo motors, actuating two of the three legs, can move the upper platform on the surface of a sphere centered in the lower platform. While the device can be used both in a grounded or handheld configuration, this paper focused on the former. In this case, the lower base is firmly attached to an external support, and the upper platform moves so as to provide the user with haptic feedback. We carried out an experimental evaluation of our device. 14 human subjects were asked to place their hand on the upper platform of the device and recognize the pattern being displayed. The device was commanded to render 8 linear and 3 shape patterns, which can be used to provide, e.g., navigational information. Results show a very high recognition rate for all patterns, with the highest and lowest recognition rate being 97.1% and 65.8%, respectively, with a chance level of 12.5%. For the linear patterns recognition, it was easier to recognize horizontal patterns (lateral-median motions with respect to the user, e.g., West and East patterns), while diagonal patterns were the hardest to recognize. The sense of motion along a certain pattern did not seem to affect much performance. The most common mistake was to identify diagonal patterns as anterior-posterior ones (North and South patterns). For the shape patterns recognition, recognition rates are very high in all conditions. In this case, subjects reported to be mostly relying on the number of edges being rendered by the device.

The design of this device shows also some limitations. First, in the current configuration, as we are using servo motors, the device is not backdrivable and cannot be used as a standard impedance-type haptic interface. Second, the motors are both placed on the same platform, which is a good idea if the device is used in a grounded configuration, but not ideal if held between the two hands. In a hand-held configuration, it is better to distribute the weight equally between the two platform, which should be straightforward to do.

Given such promising results, in the next future we plan to test a wider range of patterns, evaluating, e.g., the smallest change in orientation between two linear patterns and the smallest change in the number of edges between two shape patterns users can discriminate. We will also test the device when held between two hands, as well as the performance of the device in terms of resolution, precision, repeatability, and maximum output force. This device was developed within a larger effort to design a multi-modal haptic handle for various

Fig. 1 .

 1 Fig. 1. The device is composed of two servo motors actuating an anti-parallelogram mechanism moving two platforms with respect to each other. (a) Device in its grounded configuration, with the lower platform secured to an external support and the user hand posed on the upper one. The device can be either used as a grounded interface, like in this figure and in our experiment of Sec. 3, or it can be held between the two hands (see the video for this configuration). (b) CAD of the device. The grey dotted lines show the surface of the sphere on which the upper platform moves. The spherical cap posed on the upper platform makes the interaction with the user hand more comfortable. Of course, it can be changed with other shapes according to the task at hand.

  Kinematic model on the two actuating legs.

Fig. 2 .

 2 Fig. 2. Schematic of the proposed device.

  , motor 1 is mounted at B 1 with the motor shaft rotating along B 1 O b , while motor 2 is mounted at B 2 , with its shaft rotating along B 2 O b . The motor's rotating angle is limited within (-π/2, π/2). For any bending pose of the upper platform, O p has a projection point K in the plane B 1 B 2 B 3 . Obviously, O p K is perpendicular to B 1 O b and B 2 O b , and it is in the bending plane. Moving perpendicularly from K to B 1 O b and B 2 O b , we obtain two pedal points, E, F , which lie in the shaft axis of motor 1 and motor 2, respectively. Right-angled triangle O p KE and O p KF share the same edge O p K. As O p E and O p B 1 lie in the same plane O p O b B 1 , then the rotation angle of motor 1, α 1 , is equal to O p EK, and that of motor 2, α 2 , is equal to O p F K. From a simple geometrical derivation, we obtain tan

  View of the patterns from the top of the device.

  The linear patterns considered.

Fig. 3 .

 3 Fig. 3. Linear patterns. (a) The patterns as rendered by the device, view from the top of the upper platform; (b) recognition rates of linear patterns; (c) the linear patterns. All linear patterns start and finish in the resting position of Op, moving along a circular arc having a chord of 20 cm.

  View of the patterns from the top of the device.

Fig. 4 .

 4 Fig. 4. Shape patterns. (a) The patterns as rendered by the device, view from the top of the upper platform. All shape patterns start and finish north of the resting position of Op, in the point named P. (b) Recognition rates of shape patterns;

This work has received funding from the Inria Défi project "DORNELL" and the China Scholarship Council No. 201908440309.

mobility aids, e.g., power wheelchairs, walkers, prewalkers. For this reason, we also plan to test this device effectiveness in rendering directional information when mounted on one of these mobility aids, similarly to [START_REF] Devigne | Power wheelchair navigation assistance using wearable vibrotactile haptics[END_REF].


P-A Cabaret

T Howard

C Pacchierotti

M Babel

M Marchal

Perception of spatialized vibrotactile impacts in a hand-held tangible for virtual reality

Informative and realistic haptic feedback significantly enhances virtual reality (VR) manipulation. In particular, vibrotactile feedback (VF) can deliver diverse haptic sensations while remaining relatively simple. This has made it a go-to solution for haptics within hand-held controllers and tangible props for VR. However, VF in hand-helds has solely focused on monolithic vibration of the entire hand-held device. Thus, it is not clear to what extent such hand-held devices could support the delivery of spatialized information within the hand. In this paper, we consider a tangible cylindrical handle that allows interaction with virtual objects extending beyond it. This handle is fitted with a pair of vibrotactile actuators with the objective of providing in-hand spatialized cues indicating direction and distance of impacts. We evaluated its capability for rendering spatialized impacts with external virtual objects. Results show that it performs very well for conveying an impact's direction and moderately well for conveying an impact's distance to the user.

Introduction and related work

Vibrotactile feedback (VF) is a popular haptic feedback modality for virtual reality (VR) interaction because it combines relatively low technological complexity in its implementation with a wide variety of achievable haptic effects [START_REF] Choi | Vibrotactile display: Perception, technology, and applications[END_REF]. There is evidence for a positive impact of VF on many success metrics for interactions with virtual environments (VEs), such as improved task performances [START_REF] Brasen | Effects of vibrotactile feedback in commercial virtual reality systems[END_REF], improved user immersion [START_REF] Brasen | Effects of vibrotactile feedback in commercial virtual reality systems[END_REF], increased perceived realism [START_REF] Peng | Walkingvibe: Reducing virtual reality sickness and improving realism while walking in vr using unobtrusive head-mounted vibrotactile feedback[END_REF] and increased presence [START_REF] Kaul | Increasing presence in virtual reality with a vibrotactile grid around the head[END_REF][START_REF] Cooper | The effects of substitute multisensory feedback on task performance and the sense of presence in a virtual reality environment[END_REF]. VF can be used to communicate both physical cues relating to the VE (e.g. vibrating objects [START_REF] Passalenti | No strings attached: Force and vibrotactile feedback in a virtual guitar simulation[END_REF], contacts [START_REF] Chinello | A three revolute-revolute-spherical wearable fingertip cutaneous device for stiffness rendering[END_REF], impacts [START_REF] García-Valle | Evaluation of presence in virtual environments: haptic vest and user's haptic skills[END_REF], interaction forces [START_REF] Cheng | Vibrotactile feedback in delicate virtual reality operations[END_REF], texture roughness [START_REF] Culbertson | The Penn Haptic Texture Toolkit for Modeling, Rendering, and Evaluating Haptic Virtual Textures[END_REF]) as well as abstract cues (e.g. for guidance, notification or communication [START_REF] Choi | Vibrotactile display: Perception, technology, and applications[END_REF]). Many technologies can deliver VF in VR, such as wearable [START_REF] Choi | Vibrotactile display: Perception, technology, and applications[END_REF], grounded [START_REF] Sreng | Using vibration patterns to provide impact position information in haptic manipulation of virtual objects[END_REF], hand-held [START_REF] Adilkhanov | Vibero: Vibrotactile stiffness perception interface for virtual reality[END_REF], and even mid-air haptic devices [START_REF] Howard | Pumah: Pan-tilt ultrasound mid-air haptics for larger interaction workspace in virtual reality[END_REF].

Despite the wide use of VF delivered through ungrounded hand-held devices in VR, the approach is mostly restricted to monolithic VF [START_REF] Choi | Vibrotactile display: Perception, technology, and applications[END_REF]. This has the advantage of simplicity as it requires only a single actuator. However it remains inadequate for providing spatial information. That is, cues originating from different directions relative to the user are identical and thus indiscriminable for the user, unless a mapping is created between direction and waveform parameters. Conversely, localized VF through multiple actuators is widely used to convey spatial information to the user, in particular in wearables (e.g. [START_REF] Kaul | Increasing presence in virtual reality with a vibrotactile grid around the head[END_REF][START_REF] Choi | Vibrotactile display: Perception, technology, and applications[END_REF][START_REF] De | Designing a vibrotactile head-mounted display for spatial awareness in 3d spaces[END_REF]) and surface haptics (e.g. [START_REF] Pantera | Multitouch vibrotactile feedback on a tactile screen by the inverse filter technique: Vibration amplitude and spatial resolution[END_REF]). In this paper, we begin to explore the possibilities offered by localized VF within handheld tangible objects. In particular, we focus on rendering spatialized impacts happening on a virtual hand-held object larger than the tangible held by the user. To render the impacts, the tangible houses two vibrotactile motors at its extremities. We hypothesize that by using two actuators, we can provide localized vibrotactile feedback which can inform the user about where the impact occurred on the larger virtual object they are manipulating (see Fig. 1).

In early work on impact rendering in interactions within VEs, Wellman et al. [START_REF] Wellman | Towards realistic vibrotactile display in virtual environments[END_REF] used a data-driven approach to play back recorded impact vibrations during virtual contacts on a voice-coil actuator embedded into a grounded forcefeedback device handle. Okamura et al. expanded on this, compiling a vibration waveform library for impacts generated by fitting a simplified vibration model based on an exponentially decaying sinusoid to recorded impact data [START_REF] Okamura | Vibration feedback models for virtual environments[END_REF]. Because this model (see Sec. 2.2) provided an interesting compromise between perceived realism, impact property discrimination, and computing requirements, it has since been widely adopted in interactions with VEs [START_REF] Kuchenbecker | Improving contact realism through event-based haptic feedback[END_REF][START_REF] Sreng | Using vibration patterns to provide impact position information in haptic manipulation of virtual objects[END_REF].Some work on spatialization in VR was performed by Gongora et al. [START_REF] Gongora | Experiments on two-handed localization of impact vibrations[END_REF]. They studied vibrotactile impacts delivered in a bimanual task using a pair of monolithic handheld vibrotactile devices, with the aim of rendering localized vibrotactile impacts along a virtual bar connecting both hands. There have also been a few research attempts at systems spatializing vibrotactile cues inside hand-held devices using multiple vibrotactors [START_REF] Ryu | T-hive: Vibrotactile interface presenting spatial information on handle surface[END_REF] or asymmetric vibrations [START_REF] Tappeiner | Good vibrations: Asymmetric vibrations for directional haptic cues[END_REF] but to our knowledge none have been leveraged in VR interactions.

Experimental design

Research questions and hypotheses

We seek to provide VF to render impacts between one manipulated virtual object and other virtual objects in a VE. Our question concerns the extent to which spatializing impact cues by distributing them between two actuators embedded in a cylindrical tangible handle (see Fig. 1-B) is effective in providing users with Values were determined based on literature and a pilot study. We do not consider any impact occurring within the hand, hence the null values between x th and xp.

information on impact direction. We also seek to understand how this approach affects perceived realism and impact properties, and whether it is compatible with existing approaches to rendering impact distance in a setup using a single actuator (e.g. [START_REF] Sreng | Using vibration patterns to provide impact position information in haptic manipulation of virtual objects[END_REF][START_REF] Gongora | Experiments on two-handed localization of impact vibrations[END_REF]). To investigate this, we compare distance and direction discrimination performances, as well as perceived realism and virtual object material properties in VR, using different impact vibration models (see Table 1). We formulate the following hypotheses: H1 Spatialization of impacts in hand by assigning impact waveforms to distinct vibrotactors will allow discrimination of impact direction, regardless of the chosen impact vibration model. H2 Impact models coding distance with more redundant parameters (see Sec. 2.2 for the details of the models) will yield better distance discrimination performance.

Rendering impacts distance and direction

We use the simplified impact vibration model introduced by Okamura et al. [START_REF] Okamura | Vibration feedback models for virtual environments[END_REF], where α(x, t) denotes the waveform amplitude at the instant t for an impact at a distance x from the hand (see Fig. 2): α(x, t) = A(x)e -β(x)t sin(ω(x)t). In realistic impacts, the peak amplitude A, decay β and angular frequency ω would all be functions of impact distance as well as impact dynamics and properties of the materials involved. However, such impact models can sometimes be less effective at communicating usable information on impact distance [START_REF] Sreng | Using vibration patterns to provide impact position information in haptic manipulation of virtual objects[END_REF]. An alternative is to select a subset of model parameters (A, β, ω) to encode impact distance, possibly leaving the remainder free for encoding other impact properties. Given these three parameters, there are seven different possibilities (see Table 1) for encoding impact distance (see Fig. 2).

Materials and methods

To investigate the formulated hypotheses, we designed a pair of experiments assessing impact direction and distance perception in VR. Hardware. Subjects sat at a table, wearing an HTC Vive Pro head-mounted display (HMD). They held the vibrotactile handle in their dominant hand which was tracked using an HTC Vive Tracker attached using an adhesive fixture to keep the palm and inside of the fingers unobstructed. They used an HTC Vive Controller held in their non-dominant hand to answer experimental questions (see Fig. 3-A). The handle was equipped with a pair of symmetrically mounted Actronika HapCoil One voice-coil actuators (see Fig. 1-A).

Experimental task. The common experimental task for both experiments was inspired from Sreng et al. [START_REF] Sreng | Using vibration patterns to provide impact position information in haptic manipulation of virtual objects[END_REF]. Subjects were asked to hold the tangible handle in their dominant hand. They observed the VE showing their virtual hand holding a virtual rod with the same diameter as the tangible handle but extending symmetrically outward 0.5m beyond the edges of the tangible handle. By moving this virtual rod up and down, it could impact a lightweight and unconstrained object at one of four distances d i = 0.05, 0.2, 0.35, 0.50m from either the thumb or the pinkie side of the hand (see Fig. 2-B,C). These impacts were rendered according to one of the impact models summarized in Table 1.

During the experiment, the impacted object was obstructed from view so as to provide no visual feedback of the impact location (see Fig. 3-C). Subjects placed the rod at the starting location, then were prompted to move it downward. On the way down, the stick impacted a first virtual object which appeared randomly on the left or right at one of the distances d i . Upon reaching the target location, subjects were prompted to return the stick to the start location and repeat the process. A second object appeared on the same side as the first, at one of the four possible distances, generating a second impact, after which subjects answered a pair of experimental questions:

Q1 Which side did the impacts occur on? (Left/Right) Q2 Was the second impact further away from the hand than the first? (Y/N)

Experimental design and protocol For achieving a shorter experiment, we split our investigation into two identical experiments containing 4 blocks each. Impacts were rendered respectively using Amp, Dec, AmpDec, AmpDecFreq for experiment 1 and Freq, AmpFreq, DecFreq, AmpDecFreq for experiment 2 (see Table 1). 24 subjects (19 m., 5 f., ages 21-30 (Mean:24.9y), 20 right-handed) took part in the study after providing written informed consent. Subjects were randomly assigned to one of the two experiments.

In each experiment, subjects first performed a familiarisation task where the VE was fully visible, showing the hand-held virtual stick and the impacted virtual objects (see Fig. 3-B). During this task we ensured that subjects moved at a similar speed, though the vibration did not depend on it. They were informed that the rod and impacted object properties might vary during the course of the subsequent experiment. Subjects filled out an initial questionnaire indicating personal data and prior experience with haptics, VR and perception studies.

Each experiment contained one block per impact model, whose order was counterbalanced between subjects. Within each block, subjects performed 3 repetitions of the task for each of the 16 combinations of impact distances occurring on either side, totalling 96 trials presented in a fully random order. Post-block questionnaires assessed perception of the stick and impacted object's material and geometric properties, their variability, and perceived impact realism.

Results

Impact directions were consistently correctly identified between 94% and 97% of the time across all impact models. Most errors occurred for pairs of low amplitude and duration stimuli.

To test for H2, subjects were assigned, for each impact model, to one of two groups based on whether they interpreted the impact model as intended (increased impact distance perceived as an increased impact distance) or in an inverted manner (increased impact distance perceived as a decreased impact distance). Inversion rates (percentage of subjects interpreting an increase in impact distance as a decrease) were around 50% for all models not involving Freq, and varied between 92% and 100% for all models involving Freq.

We then computed the 75%-just-noticeable-difference (JND) for distance discrimination as a Weber fraction for each subject by fitting cumulative Gaussians to the data. Finally, we compared the distribution of JNDs across impact vibration models (see Fig. 4). Data from the experiment 1 (Fig. 4-B) were not normally distributed, and a Friedman test showed no significant differences between conditions. Data from the experiment 2 were normally distributed, and a 2-way ANOVA showed a significant effect of impact model (F (3) = 4.132, p = 0.021) but no significant differences between participants. A post-hoc Tukey HSD test revealed the only significant difference to lie between the JNDs for the Freq and AmpDecFreq conditions (p = 0.016).

Rod properties were rated most consistent (median 2 of 7) in all conditions but Amp (median 3 of 7) and AmpDecFreq (median 4 of 7), however none of these differences were significant. The properties reported as changing between trials were rod material (Freq, AmpDec, AmpDecFreq), stiffness (all models except Dec), length (Dec, AmpDec, AmpFreq), fill (Dec, AmpDec), weight (Freq, Am-pDecFreq). Subjectively reported rod materials were dominated by "metal" and "plastic" for all models involving Amp, as well as the Dec model, with qualifiers such as "resonating" and "tube". Models involving Freq but not Amp yielded more "wood" and "plastic" responses, with qualifiers such as "soft", "damped" and "warm". AmpDecFreq yielded an almost even mix of all three material categories. Realism was consistently rated as average across all models (median 3 of 7) and was considered slightly variable across all models (median 3 of 7).

Discussion

The impact direction identification rates between 94% and 97% indicate that regardless of the chosen impact model, spatializing the impacts between two actuators allowed subjects to correctly and intuitively identify the side on which the impact occurred with a high degree of accuracy. Hypothesis H1 is therefore verified. Looking at inversion rates, it is interesting to note that all models involving Freq tended to be systematically inverted (92% to 100% of subjects perceived an increase in distance as a decrease) which would indicate the evolution of ω may be the cause for this.

Weber fractions for distance discrimination were consistently high across all impact models except AmpDecFreq (m=0.17), ranging from 0.6 (DecFreq, experiment 2) to 1.32 (Amp, experiment 1). This indicates that while distance discrimination was mostly possible, it was far from an easy task. The only statistically significant difference observed (Freq-AmpDecFreq, experiment 2) is in favor of hypothesis H2, and the mean JNDs seem to also support the hypothesis. However, given the poor performance of AmpDecFreq in experiment 1 and the fact that none but one of the differences are statistically significant, we cannot conclude that H2 is supported. This may be due to H2 being wrong, or to flaws in the stimulus or experimental task design. If H2 is not verified, there may be a lot of headroom to encode various impact properties by distributing them across different parameters without adversely impacting performance.

The high inversion rates due to using ω as a parameter led us to hypothesize that models combining ω with A, β or both may have been confusing to half the subjects that did not invert their interpretation of Amp and Dec. This hypothesis cannot be easily tested because subjects that performed Amp and Dec did not perform AmpFreq and DecFreq. Yet, analysing the results from experiment 1 revealed that 6 out of 12 subjects had inverted both Amp and Dec while 5 of 12 had not (the remaining subject inverted only one of both models). By looking at the JNDs for each of these groups of subjects in the AmpDecFreq condition, we note that the group that inverted both Amp and Dec performed better at AmpDecFreq (JNDs: 0.07 to 1.96, mean=0.83) than the group that did not invert Amp and Dec (JNDs: 2.15 to 9.33, mean=4.63). This would tend to support our interpretation and argue for the need to redesign the function ω(x) in our rendering approach. Furthermore, it may be necessary to consider the frequency dependence of vibration amplitude perception in such a redesign. However, given the very small sample size, this conclusion must be seen as tentative.

The spread in JNDs between subjects indicates a large inter-subject variability in the ability to perform the task. During the experiment, several subjects noted that the task was really difficult until they "chose" a way to understand the mapping of the stimuli to impact distance. Thus, we believe this variability shows that subjects displayed different capacities for adapting to the difficulty of the experimental tasks and choosing an effective response strategy. This means that the haptic representation of impact distance is far from intuitive or natural with the chosen models, although AmpDecFreq shows some promise in experiment 2. This may indicate poor model design, or the fact distance discrimination is really hard without any context such as e.g. visual feedback of impacts.

All models were perceived as equally (un)realistic, indicating that either the impact model used is unrealistic, that spatialization impacted realism, or both.

Conclusion and perspectives

We presented an investigation into the use of spatialized in-hand vibrotactile feedback for VR interactions. Our study focused on the ability of a handle equipped with two vibrotactors to deliver realistic, discriminable and understandable sensations of impacts through which users could determine the location (direction and distance) of impacts on a virtual manipulated object.

We determined direction identification scores, JNDs for impact distance, and perceived impact realism for 7 impact vibration models.Results showed excellent direction identification performances, but distance discrimination performances were mediocre. Impacts were perceived as only moderately realistic, which may be due to the impact models studied as well as our spatialization technique.

In future work, we plan to investigate whether differently distributing the vibrations between both actuators can improve perceived realism and consistency between impacts while preserving distance and direction discrimination performance. This study also highlighted certain avenues for improving the perception of vibration impact models which we intend to investigate. Finally, we plan to extend the approach to 2D and 3D spatialization using more actuators.

Because of the good direction discrimination performance observed in this initial study, we believe there is potential for using multiple actuators in manipulated tangible objects or controllers. This seems particularly promising for VR applications which could benefit from the use of in-hand directional cues.

Fig. 1 .

 1 Fig. 1. (A) Close up CAD of the tangible handle being held by the user's avatar; (B) User manipulating the handle in VR; (C,D) VR interactions causing impacts at different distances and in different directions from the hand.

Fig. 2 .

 2 Fig. 2. (A) Virtual rod manipulated in the experiment with 4 possible impact distances extending symmetrically around the virtual hand. x th and xp respectively denote the thumb and pinkie side actuator positions. (B) Possible evolution of vibration amplitude A, decay β and frequency f = 2πω as a function of impact distance for both actuators.Values were determined based on literature and a pilot study. We do not consider any impact occurring within the hand, hence the null values between x th and xp.

Fig. 3 .

 3 Fig. 3. (A) A subject performing the experiment. (B) VR view of the familiarization task, where the impacted objects are visible. (C) VR view of the actual task, where the impacted objects are hidden and only haptic feedback of impacts is provided.

Table 1 .

 1 Impact vibration models for encoding impact distance x from the hand studied in our experiments. The model names indicate the vibration parameters that vary as a function of impact distance, with Amp referring to amplitude A, Dec referring to the decay β and Freq referring to the frequency ω. (Left) Models used in experiment 1; (Right) Models used in experiment 2; AmpDecFreq was common to both experiments.

Fig. 4 .

 4 Fig. 4. JNDs for impact distance, expressed as Weber fractions, for both experiments.

This project has received funding from the EU Horizon 2020 program, grant agreement No 801413, project "H-Reality"; and from the Inria Défi project "DORNELL".