M.-Anas Knefati
email: maknefati@hotmail.com

Pierre E Chauvet
email: pierre.chauvet@uco.fr

Sylvie Nguyen
email: sylvie.nguyenthe@gmail.com

Bassam Daya

Sylvie N'guyen

Reference Curves Estimation Using Conditional Quantile and Radial Basis Function Network with Mass Constraint

Keywords: Quantile regression, Radial basis function neural network, Reference curve, Brain maturation 1 Motivation

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

The reference curves are part of the basic tools of the physician practices, in that they can decide on the vulnerability of an individual against a certain disease. Their construction, which is based on measurements made on a given set of individuals, shows a partition between normal and at risk individuals. It was during a project in pediatric neurology that we have developed a general method for construction of reference curves based on conditional quantiles and Radial basis function networks (RBFN). Because cerebral injury in newborns tends to be clinically silent, tools and techniques for neurological evaluation are essential. The electroencephalogram (EEG) is such a bedside, non-invasive and low cost technique. The EEG consists in recording the spontaneous electrical activity of the brain through several electrodes placed on the scalp. In preterm infants, the normal background EEG activity has the unique characteristic of being spontaneously discontinuous with periods of electrical activity alternating with periods of quiescence -called interburst intervals (IBIs). Figure 1 illustrates this phenomenon, which is normal in infancy. The diagnostic and prognostic value of neonatal EEG abnormalities in the preterm infant are well established and the IBI duration have been shown to be related to abnormal brain maturation in preterm infants [START_REF] Watanabe | eonatal EEG: a powerful tool in the assessment of brain damage in preterm infants[END_REF]. In a recent study comparing maturation of cerebral activity and cortical folding, IBI duration was the only parameter significantly linked to morphologic brain maturation [START_REF] Biagioni | Maturation of cerebral electrical activity and development of cortical folding in young very preterm infants[END_REF].

We have developed inside a dedicated web portal a Java application allowing the physician to extract automatically all the IBIs from an EEG. In this area, artificial neural networks are used generally for features extraction and classification, particularly for brain computer Fig. 1 Normal discontinuous EEG tracing in a preterm infant at 26 weeks of gestationnal age consisting in burst of electrical activities separated by interburst intervals (IBI) interfaces (see for example [START_REF] Mateo | Dynamic fuzzy neural network based learning algorithms for ocular artefact reduction in EEG recordings[END_REF][START_REF] Wang | Harmonic mean of Kullback-Leibler divergences for optimizing multi-class EEG spatiotemporal filters[END_REF]). In this work, the features (IBIs) are already detected and the objective is to build reference curves for brain maturation using these data, easy to recalculate from new sets of EEG and robust against artefacts and entry mistakes (e.g. age).

Our paper is structured as follows. First we start by reviewing the state of the art in Sect. 2. Then, after a few reminders about the basics of conditional quantile, we introduce in Sect. 3 the RBFN modeling of the conditional quantile curve; the learning algorithm is based on a constructive approach with a control of the integral of the model over its input domain (its mass), and the coupling with the conditional quantile estimator of Yu and Jones (QYJ), to obtain the so-called QRBFc model. Experiments are carried out in Sect. [START_REF] Chen | Orthogonal least squares learning algorithm for radial basis function networks[END_REF], where we compare on simulated data with known properties the QYJ, the QRBFc, the quantile regression neural network models (QRNN), and the possible improvements of the mass constraint. Finally we present and discuss in Sect. 5 the brain maturation reference curves obtained with the different models, based on the correlation between the mean length of interburst intervals in the newborn EEG and its age. We conclude this work in Sect. 6, providing our future work directions.

State of the Art

Conditional quantile regression has gained particular attention during the recent three decades due to their useful applications in various disciplines, such as finance, economics, medicine, and biology. See for exemple Fan and Gijbels [START_REF] Fan | Estimation of conditional densities and sensitivity mesures in nonlinear dynamic systems[END_REF] and Cai and Wang [START_REF] Cai | Nonparametric methods for estimating conditional VaR and expected shortfall[END_REF]. Of particular interest is the conditional median which is more explainable and more robust than the mean regression function for asymmetric conditional distributions. For τ ∈ (0, 1), the quantile regression function gives the τ th quantile q τ (x) in the conditional distribution of a response variable Y given X = x. It measures the effect of covariates not only in the center of the population, but also in the lower and upper tails. For x varying in a given real interval, q τ (x) is a reference curve that predicts vulnerability of an individual with the probability τ or 1τ as its associated measured pair (X = x, Y = y) is below or above the curve. A classical approach to evaluate the conditional quantile function from a sample (X i , Y i ) on a discrete finite set is the local linear quantile regression method from Yu and Jones (denoted QYJ) in [START_REF] Yu | Local Linear quantile rregression[END_REF]. As its name suggests, q τ (x) is approximated by a piecewise linear curve, resulting is a broken line that emphasizes local fidelity to the detriment of regularity, difficult to interpret as a reference chart.

Cannon [START_REF] Cannon | Quantile regression neural networks: implementation in R and application to precipitation downscaling[END_REF] developed in R the QRNN package (for Quantile Regression Neural Network) based on the work from Taylor [START_REF] Taylor | A quantile regression neural network approach to estimating the conditional density of multiperiod returns[END_REF] in the field of time series forecasting. J. W. Taylor used a one hidden-layer feedforward neural network (FNN) to fit a non-local quantile model. FNN where successfully used in several domains like robotics and image processing (see for example [START_REF] Daya | A neural control system of a two joints robot for visual conferencing[END_REF][START_REF] Daya | A multilayer perceptrons model for the stability of a bipedal robot[END_REF]) as non linear parametric models. FNN realizes a global approximation of the data, and by varying the number of neurons it is possible to slide between high/low precision and smoothness of the quantile modelling. However, standard non-linear FNN with sigmoid transfer functions does not allow local settings on some parts of the data, and the learning process can be very slow. RBFN has the advantage to allow a faster constructive approach since each neuron puts emphasis on local data: neurons can be added iteratively, as in Chen et al. [START_REF] Chen | Orthogonal least squares learning algorithm for radial basis function networks[END_REF].

we explain how to calculate them introducing constraint on the weights. Finally, we end this section by explaining our method to improve the estimation of the conditional quantile with a RBFN network, leading to the QRBFc method.

An Overview of the Conditional Quantile

Let (X, Y ) be a bivariate random variable, F(y|x) the conditional distribution function of Y given X = x,andτ ∈ (0, 1).Theτ th conditional quantile, noted by q τ (x),isgivenby

q τ (x) = inf{y ∈ R : F(y|x) ≥ τ }≡F -1 (τ |x) (1)
or equivalently, by

q τ (x) = arg min θ E{ρ τ (Y -θ)|X = x}, ( 2
)
where ρ τ is the "check function" ρ τ (u) = 0.5(|u|+(2τ -1)u).

The building of predictive intervals is an important application of quantile regression. Suppose the observations can be modelled as Y i = m(X i ) + γǫ i where m(.) is an unknown function, γ is unknown reel parameter and the residual ǫ i are uncorrelated random variables with zero mean and one variance. Then the quantile function can be written as

q τ (x)) = m(x) + γ F -1 ǫ (τ ), (3)
where F ǫ (.) denotes the distribution function of ǫ. A predictive interval is an interval that predicts, with certain coverage probability, the future value of the response variable Y for a given covariate X = x. The pairs of extreme conditional quantiles q inf (x) and q sup (x) map out a conditional prediction interval within which one expects the majority of individual points to lie. These "reference curves" are popular in medicine (see, e.g, Cole in [START_REF] Cole | Fitting smoothed centile curves to reference data[END_REF]) and have provided a stimulus for much of the recent statistical work in this area. Parametric techniques for estimating conditional quantile can be efficient if the underlying functions are correctly specified. But a misspecification may cause serious bias, and model constraint may distort the underlying distribution. Therefore, we will concentrate in this article on the nonparametric quantile regression, with the advantage that little or no restrictive prior information on functionals is needed.

Nonparametric estimation of conditional Quantile has been tackled by several authors, with direct and indirect methods. Direct methods use the "check" function ρ τ , taking roots in definition [START_REF] Cannon | Quantile regression neural networks: implementation in R and application to precipitation downscaling[END_REF]. See, for instance, Fan et al. [START_REF] Fan | Robust Non-parametric function estimation[END_REF] and Yu and Jones [START_REF] Yu | Local Linear quantile rregression[END_REF] for a local linear estimate of conditional quantiles. Indirect methods, inspired from (1), are performed in two steps: the estimation of the conditional distribution is performed first, then the inverse of the obtained estimator is used to estimate the desired conditional quantile. Examples of indirect methods are provided by Cai and Wang [START_REF] Cai | Nonparametric methods for estimating conditional VaR and expected shortfall[END_REF]. Our work is based on the direct estimator proposed by Yu and Jones [START_REF] Yu | Local Linear quantile rregression[END_REF], because of its robustness in the y direction. The idea is to approximate the unknown τ th conditional quantile q τ (x) by the linear function

q τ (z) = q τ (x) + q ′ τ (x)(z -x) ≡ a + b(z -x)
for z in the neighborhood of x. This motivated us to define an estimator by setting qτ (x) = â(x), with

â(x), b(x) = arg min (a,b)∈R 2 n i=1 ρ τ (Y i -a -b(X i -x)) × K ( x -X i h ), (4)
where K is a kernel density function and h is the smoothing (also called bandwith) parameter, which is a nonnegative number controlling the size of the local neighborhood. Yu and Jones have proposed a "rule-of-thumb" for selecting h:

h τ = h mean τ(1 -τ) φ(Φ -1 (τ )) 2 1 5 , (5)
where h mean is the optimal choice of h for regression mean estimation, φ and Φ are the standard normal density and distribution functions. They recommended to use the technique proposed by Ruppert et al. [START_REF] Ruppert | An effective bandwidth selector for local least squares regression[END_REF] for selecting h mean .

The solution of the problem (4) can be obtained by using the iteratively reweighted least squares algorithm, see for exemple Yu [START_REF] Yu | Smooth regression quantile estimation[END_REF]. We will denote by QYJ this algorithm (as well as the resulting quantile curves) used in conjunction with the rule-of-thumb from Yu and Jones in what follows.

The RBFN Neural Network Model

A RBFN can be viewed as a FNN with a specific structure that allows an easier constructive approach. It has one hidden layer in which each neuron computes its output using a radial basis function (RBF) receiving the inputs, and an output layer which builds a linear weighted sum of hidden neuron outputs and supplies the network's response.

An RBF function is a function φ :[ 0, ∞) → R that depends only on the distance from some point c, called a center, so that it has the form φ( xc ).I no t h e rw o r d s ,ar a d i a l basis function is radially symmetric with respect to a given norm. We generally choose the

Gaussian function φ( x -c ) = e -x-c 2 2 .
Because our goal is to produce reference curves, i.e. real functions, we use RBFN with only one output neuron. Its model is given by:

f (x) = N j=1 ω j φ( x -c j ), (6)
where x ∈ R d , N is the number of hidden neurons, ω j (1 ≤ j ≤ N ) are the weights of network output layer, c j (1 ≤ j ≤ N ) is the center of the jth hidden neuron, φ(x) is a gaussian RBF function and . denotes the distance function that is taken, in general, to be the Euclidean norm. In this work we use the Mahalanobis distance, because it takes into account the correlations inside the data set and is scale-invariant; in our case, it contributes to a better adjustment of the neurons width with the data. This distance is defined as

x -c = (x -c) t Σ -1 (x -c), ( 7
)
where Σ is the covariance matrix of (x 1 ,...,

x n ),andx i ∈ R d (i = 1,…, n) are the training data input.
We use in this work an iterative design of the network: neurons are added one at a time until the mean sum-squared error (MSE) falls beneath an error goal or a maximum number of neurons has been reached. The choice of the number N of neurons on the hidden layer of FNN and RBFN networks is crucial. A low number means a very poor performance, or fidelity, of the network. Instead, a large number of neurons will allow the network to fit exactly all the data (resulting in a very low MSE), including noise and biased observations: its regularity will be low. Because N is an integer defining the structure of the network, it cannot be adjusted like the synaptic weights. It exists for FNN some meta-heuristics (like genetic algorithms) to adjust N , but this is a slow process. Unlike FNN, RBFN authorizes a faster constructive approach because of the locality of the radial basis transfer function in the hidden layer.

Let x i ∈ R d be the training data inputs and y i ∈ R be the training data outputs for i = 1,…, n. We now explain our RBFN learning algorithm.

Computation of the Centers

The algorithm for computation of the centers can be stated as follows:

-Provide the error threshold (err ) and the maximum number of neurons Nmax (Nmax ≤ n). -Initialize N to 0: the initial network does not have any radial functions. -Store the output vector (y 1 , ..., y n ) in an other vector (y * 1 , ..., y * n ). -While the total error ( 1

n n i=1 ( f (x i ) -y i ) 2
) is greater than err and N < Nmax do:

calculate the network output f (x i ) using ( 6) for each input and the error r i =| f (x i )y * i | for i =1, ..., n; find the input x l that causes the greatest error, with l = arg max i∈{1,...,n} r i ; add a radial basis function whose center is this entry (c N = x l , N = N + 1); recalculate the weight vector W = (ω 1 , ..., ω n ); recalculate f (x i ), i =1,...,n; initialize to 0 the output value that causes the greatest error, i.e. y * l = 0and f (x l ) = 0, so it will not be chosen again.

We determined here the number of hidden neurons N and the centers c j . We explain in the next section the computation of the weights.

Computation of the Weights

Due to the fact that the mapping from hidden layer to output layer is linear, the weights computation becomes a linear problem. Minimization of the MSE yields to the well-known least square solution:

W = (Φ t Φ) -1 Φ t Y, (8)
where Y = (y 1 , ..., y n ) t and

Φ = ⎛ ⎜ ⎝ φ(x 1 -c 1 ) ••• φ(x 1 -c N ) . . . . . . ••• . . . φ(x n -c 1 ) ••• φ(x n -c N ) ⎞ ⎟ ⎠
A first possibility for a robust calculation of the weights is to use Tikhonov regularization (also called ridge regression). Another suggestion to avoid possible numerical difficulties due to a singular or near-singular matrix, is to use the pseudo-inverse of a matrix which generalizes the inverse of a squared matrix. Any matrix can be factored using its singular value decomposition (SVD) from which its Moore-Penrose or generalized inverse can be obtained.

After testing these two methods we finally prefer the pseudo-inverse method, that provides the same results with far less computation time.

Computation of the Weights using a Mass Constraint

The objective in this section is to modify the previous algorithm to control the mass c of the RBFN model (its integral over the input domain). This is easier than controlling the energy of the RBFN (integral of the squared function), since we will be able to bring a linear constraint on the weights. In the case where the model is positive for all x, the mass of the model is equivalent to its norm L 1 ,andc appears as a global smoothing parameter. We calculate the weights of the model [START_REF] Daya | A neural control system of a two joints robot for visual conferencing[END_REF]

such that it satisfies R d f (x)dx = c.W ehave R d f (x)dx = N j=1 α j ω j , where α j = R d φ( x-c x j σ j
)dx for j = 1, ..., N . In our case, φ( x -c x j ) = exp(-

x-c x j 2 2
) is a Gaussian function, then α j = α = ( √ π) d for all j = 1, ..., N . Therefore, for calculating the weights under the assumption that R d f (x)dx = c we obtain the linear constraint:

α N j=1 ω j = c, with α = ( √ π) d .
To solve this least squares constrained problem, let introduce the Lagrange function given by

L(W,λ) = 1 2 n i=1 y i -f (x i ) 2 + λ ⎛ ⎝ α N j=1 ω j -c ⎞ ⎠ = 1 2 n i=1 ⎛ ⎝ y i - N j=0 ω j φ ij ⎞ ⎠ 2 + λ ⎛ ⎝ α N j=1 ω j -c ⎞ ⎠ ,
where λ is the Lagrange multiplier, φ ij = φ(

x i -c x j 2
) and W is the vector of weights. The necessary conditions for optimality are

∂ L(W,λ) ∂ W = 0 ⇒- n i=1 ⎛ ⎝ y i - N j=1 φ ij ω j ⎞ ⎠ φ ik + d 0 λ = 0 and ∂ L(W,λ) ∂λ = 0 ⇒ α N j=1 ω j -c = 0.
Thus, we have to solve the N + 1 following equations with N + 1 variables:

N j=1 n i=1 φ ij φ ik ω j + αλ = n i=1 φ ik y i , k = 1, ..., N α N j=1 ω j = c
Like in paragraph III-B2, we use the SVD method to solve these equations.

Improvement of Conditional Quantile Estimation with a RBFN Network

We explain below our algorithm to improve the QYJ estimator using the RBFN model. Let (X 1 , Y 1 ), ..., (X n , Y n ) be a random sample from (X, Y ), and assume the homoscedastic model:

Y i = m(X i ) + γǫ i , i = 1, ..., n, (9)
where m(.) is an unknown function, γ unknown reel parameter and the residuals ǫ i are uncorrelated random variables of known distribution with zero mean and one variance. Therefore the conditional quantile function is defined by:

q τ (x) = m(x) + γ F -1 ǫ (τ ). (10)
The proposed algorithm is explained as follow :

-Estimate first the conditional quantile function q τ (x) using the QYJ estimator with the Yu and Jone method to select its smoothing parameter; we get q QY J τ (x). -F r o m( 10), we have

m QY J (x k ) =q QY J τ (x k ) -γ F -1 ǫ (τ )
, where γ 2 is the variance of the data outputs, and x k , k = 1, ..., n 0 are either regularly spaced values along x-axis or the values (X 1 , ..., X n ) taken in the random sample itself.

-Use the RBFN model with entries x k and y k =m QY J (x k ) to approximate the unknown function m(x) to obtain m QRBF (x). -Plug-in m QRBF (x) and γ in [START_REF] Ruppert | An effective bandwidth selector for local least squares regression[END_REF] to get the improved RBFN conditional quantile estimator:

q QRBF τ (x) =m QRBF (x) +γ F -1 ǫ (τ ) (11)
for all x in the input space.

We denote by QRBF (when we do not use the mass constraint) and QRBFc (when using the mass constraint) this algorithm and the resulting quantile curves in what follows.

Simulations with a Known Underlying Model

The following applications have been made using the R software, with the package "quantreg" from Koenker [START_REF] Koenker | Quantreg R package[END_REF] for our implementation of the QYJ algorithm and the QRNN R-package implemented by Cannon [START_REF] Cannon | QRNN R package[END_REF].

Consider the model

Y i = exp(-X i ) + exp -4(X i -1) 2 + ǫ i
for i = 1, ..., n,where{ǫ i } are independent and identically N (0, 1) random variables, and {X i } are exponential random variables with mean 1 independent from {ǫ i }. Then the true conditional quantile function is given by

q τ (x) = exp(-x) + exp -4(x -1) 2 + F -1 ǫ (τ )
, where F ǫ is the cumulative N (0, 1) distribution function.

First, we compare for different sample sizes n and different quantile values τ the performance of the QRBF estimator without constraint with the QYJ estimator and the QRNN estimator with 5 neurons on the hidden layer (after several tries, 5 neurons appear to be a good compromise between the error goal and the calculation speed). For this purpose, we ran 100 replications per experiment for several combination of the parameters values n and τ .We also demonstrate the performance of the estimators in terms of the mean absolute deviation error (MADE):

MADE = 1 n 0 n 0 j=1 |q τ (x j ) -q τ (x j )|,
where x j , j = 1, ..., n 0 are regularly spaced values between 0 and 3.5. τ=0.9

Fig. 2 Mean absolute deviation error boxplots obtained for sample sizes n = 300 and n = 500, with the two quantile values τ = 0.1andτ = 0.9. QY J i denotes the Yu and Jones QYJ estimator, QRBF i its QRBF improvement, and QRNN i denotes the QRNN estimator, where i equals 1 for size 300 and 2 for size 500 We can see on Fig. 2 the boxplots of the 100 MADE values for the QYJ, QRBF and QRNN estimators, with τ = 0.1andτ = 0.9. The QRBF estimator is better than its Yu and Jones counterpart and the QRNN estimator : (i) the MADE median of QRBF is lower than QYJ and QRNN; (ii) the MADE spread of QYJ and QRBF is smaller than QRNN. The computation time for the QRNN model is an average of 3 times that of QRBF and QYJ.

We then studied the impact of the mass constraint on the MADE median and on the number N of neurons in the hidden layer. We denote N c the number of hidden neurons obtained by the constructive learning method under the constraint R f (x)dx = c.W eran 100 replications for each simulation with different values of c. We found that the median of N c values is smaller than the median of N values for all values of τ and some values of c.

In the Table 1,weseethemedianofMADE, N and N c obtained for several values of c,in the case of τ = 0.9 and sample sizes equal to 300 and 500, after 100 simulations. It appears clearly that controlling the RBFN mass (i.e calculating weights under constraint) uses, for some values of c, less neurons and the MADE obtained is nearly the same one obtained from RBFN without constraint. In our example, we can see from the Table 1 that c = 10 or c = 20 could be a suitable choice, but we think, according to this table, that an optimal choice should exist. This is the center of our actuel research. The Fig. 3 illustrates, also, this fact for five quantile values, two sample sizes and three values of c. We think that the mass constraint act like a smoothing or regularization parameter and it would reduce the variance of the RBFN estimator by keeping, nearly, the same MADE error obtained by RBFN without constraint.

In Table 2, we put the median of MADE and computation time for every estimator, with the sample size equals to 500 and the quantile order equals to 0.9. According to this table we can see that the QRBF estimator (with or without mass constraint) reduces in this example the relative error by about 16 % for QRNN and for QYJ. Moreover, introducing the mass constraint in the QRBF estimator improves its computation time by reducing significantly its number of neurons. For this reason, we recommend to calculate the RBFN weights using the mass constraint as explained in paragraph (3.2.3), as we reduce the error against the QYJ estimator while we improve the speed against other neural network estimators.

Application to the Building of Brain Maturation Reference Curves

The next results were obtained from the automatic analysis of 416 EEG recorded between years 2003 and 2004 on newborn infants, with relative informations such that the infant outcome, age (from birthday and gestational), other neurological evaluations, etc. Moreover, each EEG in this database was visually analyzed and categorized as "artefact", "normal", "doubtful" or "pathological". A specific algorithm was built during this research to detect IBIs using an approach mimicking visual analysis. The algorithm is designed to study on each channel the variation of the signal estimated variance between contiguous short timewindows. An IBI on one channel is detected if this variation is lower than a given threshold and EEG (global) IBIs are finally computed as the intersection of channel's IBIs. More precisely, our algorithm is described by the following sequence of operations:

-each channel is filtered at 50Hz (that corresponds to frequency of electricity supply within European Union) with a second order Butterworth designed IIR filter; -each channel is smoothed with a moving window using the simple average; -each channel is processed to produce an estimate of the standard deviation on overlapped windows, a standard deviation series. -this is a two steps operation; (i) for each resulting standard deviation series, if the difference between two successive values is lower than a given V T threshold (in μV )the corresponding time intervals are aggregated and an IBI is defined by an aggregated time interval if it lasts at least m 1 seconds; (iii) finally, IBIs separated by less than m 2 seconds are regrouped. -the intersection of IBIs between all the EEG channels is computed, and only the IBIs of a length greater than m 3 seconds are retained.

By comparisons between the IBIs marked by the specialist and the IBIs detected by this algorithm, we set the parameter values to V T = 15 μV , m 1 = 1s,m 2 = 0.5 s and m 3 = 1s.

For each EEG the mean IBI length is computed and plotted versus the gestational age of the infant. From the couple (age , mean IBI length) we have traced the reference curves using the QRNN, the QYJ and the QRBFc estimators. Figures 4 and5 shows the results obtained with two levels of the quantile (0.1 and 0.9). It clearly appears that the improved QRBFc estimator (Fig. 4) produces smoother curves, easier to use, than the QYJ estimator.

The QRNN estimator produces very different cuves from one try to another as the learning algorithm is stochastic. The QRNN curves are generally less smooth, with sometimes a totally flat curve for τ = 0.1, and other times an inversion between the lower and the higher curves. The two curves define three areas allowing one to build three disjoint EEG sets denoted S A , S IN and S U , by taking EEG whose coordinates (age, mean IBI length) are respectively above the higher curve, between the two curves inclusively and under the lower curve. The union of these three sets is equal the initial total set, and the population at risk is defined by the set S A . To compare these sets we have first calculated for each one the mean of the EEG mean IBI lengths and the mean of the ages (from birthday): the results are provided in Table 3. One can see that the sets S A defined by QYJ and QRBFc have a longer mean IBI length by about three times that of the total set. This is perfectly consistent with the following clinical outcome: the higher the mean IBI length, the greater the risk of abnormal maturation is important. The mean age of the population in S A is about 203 days for QYJ and 216 days for QRBFc. Clinical studies have shown that the older the age, the less the IBIs must be long: in an individual with normal cerebral development, the IBIs should disappear. Therefore the QRBFc best defines the population at risk: its mean age is slightly higher than the mean age of the total population, when QYJ provides a younger population for which it is less unusual to have longer IBIs. We have computed too the percentages of EEG for each set in each 4). The pathological individuals represent 12 % of the total population, 37 % in the S A set defined by QYJ and 54 % in the S A set defined by QRBFc. This confirms the improved ability of the QRBFc model to build the population at risk.

Conclusion and Perspectives

Our approach to build reference curves is based on the nonparametric linear quantile regression from Yu and Jones. Its advantages are its robustness to outliers and measurement errors. Moreover, the fact that it is non-parametric allows us to construct the reference curve without a priori assumption. We improved ease of use and performance by reshaping the QYJ estimator with a RBFN network. Indeed, the network is defined for all values of x in the input space: whatever the input x presented by the user, the network produces a prediction without going through all the steps of QYJ algorithm. Somehow, we constructed a parametric model based on a non-parametric approach.

We have shown with simulations that the median (and to a lesser extent the spread) of the mean absolute deviation error obtained for 100 replications with the QRBF estimator is less than that obtained with the QYJ and QRNN estimators for different values of the quantile level τ . We have also find in several simulations that the QRNN estimator is not stable specially when n is small. We introduced in paragraph II-B3 a constraint on the "estimator's mass", i.e. its integral value over the input space, that can act like a smoothing or regularization parameter if the estimator is a positive function. This mass constraint becomes a linear constraint on the weights that is easily integrated in the algorithm of the RBFN construction, with very low additional computational cost. This constraint significantly reduces the number of hidden layer's neurons for a given range of the values of the constraint c: the computation time for plug-in a RBFN neural network in the classical quantile regression approach becomes marginal when the efficiency in term of error gains ground. Further work will be carried out to improve the qualities of the QRBFc estimator by automating the settings of c.

Concerning the IBI curves reflecting the brain maturation, further studies are in progress in order to establish the relationships between these values and the neurological outcome of the infants. We ought to propose reference curves that could help the medical physician in the assessment of the EEG and the neurological status of the preterm infants.

Fig. 3

 3 Fig.[START_REF] Cai | Nonparametric methods for estimating conditional VaR and expected shortfall[END_REF] Boxplot of the number of hidden neurons in the QRBF estimator for both cases: computing the weights without constraint (N ) and with constraint (N c with c = 10). First row: sample size equals to 300; second row: sample size equals to 500

Fig. 4 Fig. 5

 45 Fig.[START_REF] Chen | Orthogonal least squares learning algorithm for radial basis function networks[END_REF] Brain maturation reference curve for age in days (quantile levels τ = 0.1and τ = 0.9) obtained with the improved QRBF 10 estimator using a regular grid with 100 points and the automatic IBI analysis of 416 EEG recorded in[2003][2004]

Table 1

 1 Median for 100 values of N , N c and MADE with mass constraint c

	Sample size	τ = 0.9	N	N1	N5	N10	N20	N50	N100	N500	N1000	N5000
	n = 300	Neurons 82 MADE 0.242 0.238 0.24 80 78	58 0.24	56 0.24	59 0.24	76 0.24	100 0.284 0.329 2.962 100 100
	n = 500	Neurons 102 MADE 0.197 0.196 0.199 0.198 0.198 0.198 0.199 0.218 0.254 0.757 87 92 62 51 59 88 166 166 166

N denotes the median number of neurons on the RBFN hidden layer with no mass constraint and Nc this same median number using a mass constraint equals to c

Table 2

 2 cfor the different estimators; QYJ denotes the classical Yu and Jones linear quantile estimator, QRBF its improvement using the RBFN without mass constraint, QRBFc the RBFN with mass constraint, and QRNN denotes the quantile regression feedforward neural network from A. J. Cannon (τ = 0.9 and sample size:

	n = 500)
		QRNN	QYJ	QRBF	QRBF 10
	Error	0.237	0.237	0.198	0.198
	Time	2.69	0.91	1.385	1.01

Table 3

 3 Mean IBI's length and age in each EEG set; the Total row is the initial set of 416 EEG

	EEG set	Mean IBI length (s)	Mean age (days)
	Total	13.3369	216.0832
	QYJ
	S A	39.4740	202.9307
	S IN	9.1188	219.2450
	S U	2.5125	218.0835
	QRBF
	S A	40.4922	216.1256
	S IN	11.9003	215.1669
	S U	1.1202	222.3292

Table 4

 4 Percentages of EEG in categories "artefact", "normal", "doubtful" or "pathological" for each set; the Total row is the initial set of 416 EEG

	EEG set	Artefact	Normal	Doubtful	Pathological
	Total	1.9231	65.8654	20.1923	12.0192
	QYJ
	S A	03 3 .3333	29.3333	37.3333
	S IN	6.3291	84.8101	5.0633	3.7975
	S U	1.1450	69.4656	22.1374	7.2519
	QRBF
	S A	02 5 .6410	20.5128	53.8462
	S IN	10.4167	75.0000	6.2500	8.3333
	S U	0.9119	69.3009	22.1884	7.5988
	category normal, doubtful and pathological (results in Table

The RBFN Modeling of the Conditional QuantileIn this section, we first start with an overview of the conditional quantile. Next, the RBFN model is presented by mentioning its parameters which are the centers and the weights, and

Acknowledgments

The authors wish to express their appreciation to the three referees for their constructive and helpful comments and suggestions. The writing of this manuscript was supported by the Agence Nationale de la Recherche grant "BB EEG Platform".