L Zhang

C Cavaro-Ménard

P Le Callet

D Ge

A MULTI-SLICE MODEL OBSERVER FOR MEDICAL IMAGE QUALITY ASSESSMENT

Model observers (MOs) have been developed for the medical image quality assessment. Nowadays, numerous modern medical instruments are capable of producing 3D images, while few researchers have conducted MO studies on 3D data. In this paper, we propose a multi-slice MO when considering a relatively more realistic diagnostic task: the detectionlocalization of simulated multiple-sclerosis (MS) lesions on 3D magnetic resonance (MR) images. The jackknife freeresponse receiver operating characteristic (JAFROC) method was used to quantitatively analyse its performances and compare them with those of human observers. Our preliminary results showed that the proposed framework has the potential to approach human detection-localization task performance.

INTRODUCTION

Recently, the task-based approach [START_REF] Barrett | Foundations of Image Science[END_REF] has been widely accepted for medical image quality assessment. In task-based evaluation, one or more observers perform one or more diagnostic tasks (e.g. detection task) using a set of images obtained from the imaging systems being evaluated. The imaging system that enables observers to get a better diagnostic task performance is said to be better. Since the evaluation using human observers is time-consuming and expensive, the evaluation using model observer (MO) has become an attractive alternative [START_REF] Leng | Correlation between model observer and human observer performance in ct imaging when lesion location is uncertain[END_REF]. While numerous MOs have been proposed to perform a diagnostic task on solitary 2D images, those for 3D images have not been extensively studied yet. Compared to 2D images, 3D images provide more complete information which could allow for a better distinction between true abnormalities and noises/background structures. For some digital imaging modalities (e.g. MRI, CT), the 3D visualization of the image data is even becoming the standard for routine patient diagnostic care in different specialities (e.g. neurology, cardiology). Thus, more researchers have started to investigate 3D MOs which seems more clinically relevant.

Most existing 3D MOs are signal-known-exactly (SKE) MOs [START_REF] Liang | Image browsing in slow medical liquid crystal displays[END_REF][START_REF] Kim | A comparison of planar versus volumetric numerical observers for detection task performance in whole-body PET imaging[END_REF][START_REF] Lartizien | Volumetric model and human observer comparisons of tumor detection for whole-body positron emission tomography[END_REF][START_REF] Chen | Using the hotelling observer on multislice and multi-view simulated SPECT myocardial images[END_REF][START_REF] Gifford | A comparison of human and model observers in multislice LROC studies[END_REF][START_REF] Young | Estimating breast tomosynthesis performance in detection tasks with variable-background phantoms[END_REF] that can only perform the detection of one signal with known parameters (e.g. size and location) [START_REF] Ljiljana | Channelized hotelling observers for the assessment of volumetric imaging data sets[END_REF].

Although SKE MOs have successfully predicted human performance in SKE detection tasks [START_REF] Craig | Human-and modelobserver performance in ramp-spectrum noise: effects of regularization and object variability[END_REF], SKE is not clinically realistic. The truth is that human observers do not have such a priori knowledge in clinical practice. That is why there are increased demands of signal-known-statistically (SKS) MOs for which signal parameters are specified only statistically.

For the detection of a signal with variable location, scanning MOs have been proposed. They scan an image exhaustively, then choose the location that gives the largest test statistic (detection decision) as the tentative location and take the largest test statistic as the final test statistic. The idea has been adopted by Whitaker et al. [START_REF] Whitaker | Estimating random signal parameters from noisy images with nuisance parameters: linear and scanning-linear methods[END_REF] for 2D cases; and by Gifford et al. [START_REF] Gifford | A visual-search model observer for multislice-multiview spect images[END_REF] for 3D cases, in which a two-phase visual search paradigm was proposed to improve the scanning process as a description of human search.

For the detection of a non-symmetrical signal with variable orientation, MOs could be generalized by using efficient channels which do not rely on signal symmetry, such as those proposed in [START_REF] Witten | Partial least squares: a method to estimate efficient channels for the ideal observers[END_REF][START_REF] Goossens | The use of steerable channels for detecting asymmetrical signals with random orientations[END_REF][START_REF] Zhang | Channelized model observer for the detection and estimation of signals with unknown amplitude, orientation, and size[END_REF]. These efficient channels are 2D, but they are not limited by the dimensionality of the problem by definition. An alternative (a simpler way) is to assume that the observer knows the mean reconstructed 3D signal profile to allow some signal variabilities in the test images [START_REF] Gifford | A comparison of human and model observers in multislice LROC studies[END_REF].

In this paper, we propose an alternative to address the SKS task for 3D images. We tested the proposed MO for the detection-localization task of multiple non-symmetrical signals on MR images. Section 2 details the proposed MO. Section 3 presents the human observer study and shows the performance results of both observers. Section 4 concludes the paper and discusses possible future works.

PROPOSED MULTI-SLICE MODEL OBSERVER

3D MS lesions simulation

The lesion is modelled by a 3D Gaussian function here:

[x] p = a exp - 1 2 (p -q) t B t D -1 B (p -q) , (1)
where [x] p denotes the intensity value of the added signal at the 3-D coordinate p. The signal intensity attenuation is hereby modeled by a Gaussian function of peak amplitude a, centerd at q. Without loss of generality, the diagonal matrix D specifies the equatorial radii σ and bσ (along the x and y axes) and the polar radius cσ (along the z-axis). We also used the x-y-z convention with independent Euler angles φ, θ, and ψ to construct the 3D rotation matrix B: the rotations are applied firstly around the x-axis by φ, secondly around the y-axis by θ and thirdly around the z-axis by ψ; each basic rotation appears counter-clockwise and the axis about which it occurs points toward the observer.

An example of a set of original images and those with simulated 3D lesions is shown in Fig. 1. Note that a 3D signal appears brighter and possibly bigger on its central slice compared to its appearances on the adjacent slices.

msPCJO

The proposed MO is illustrated in Fig. 2, called multi-slice perceptually relevant channelized joint observer (msPCJO) hereafter. The msPCJO regards the detection-localization diagnostic task as two steps: (1) the first step is a global search to locate some 3D lesion candidates; (2) the second step is a cognitive analysis and interpretation of each candidate.

First step -Candidates localization

For this step, we take the assumption of conventional multislice MOs that radiologists interpret 3D images by preprocessing the image slice-after-slice (in x-y plane) and then integrate the information across the slices along z-axis [START_REF] Ljiljana | Channelized hotelling observers for the assessment of volumetric imaging data sets[END_REF].

We firstly apply a human visual system (HVS) model, visible differences predictor (VDP) [START_REF] Daly | The visible differences predictor: an algorithm for the assessment of image fidelity[END_REF], on each slice of an input volumetric image. The VDP produces a map of probabilities of detecting the differences between a reference image and a distorted image (as a function of their location in the input images), through several functions modelling various HVS properties (e.g. contrast sensitivity, orientation/frequency selectivity...). On the assumption that radiologists know what a healthy image should be like and refer to it in their minds during the diagnostic process, we took an image of a healthy person as reference image (such as images of the 1st row in Fig. 1) and the same image with simulated signals as distorted image (such as images of the 2nd row in Fig. 1). Then the VDP actually predicted the probability that the signal was visible in the radiologists' perceptual domain.

With detection probability maps for all the slices, 3D coordinates (z-coordinate is the slice number) with detection probability superior to the empirical threshold T p (T p = 0.9 in this paper, since a detection probability often quoted is 90% [START_REF] Georgiou | Probability of detection (pod) curves: derivation, applications and limitations[END_REF] are preselected. These selected coordinates {y n } 1,...,N are then classified based on a Gaussian mixture model (GMM) since the signal intensity is itself a Gaussian function (cf. Eq.( 1)) and the problem here is to inverse the mixture (de-mix) of a linear superposition of an unknown number of Gaussian distributions. We here used one of the well-known methods of classification based on the GMM: the expectation-maximization (EM) algorithm [START_REF] Dempster | Maximum Likelihood from Incomplete Data via the EM Algorithm[END_REF]. For a given class number I, the EM algorithm draws confidence ellipsoids (centers and covariances). The Bayesian Information Criterion (BIC) is then computed to assess the number of classes. The minimization of the BIC solves the over-fitting problem (likelihood tends to infinity if one class is created for each data x n ) by using a penalization term on the number of free parameters k (directly related to the number of classes):

BIC = -2 N n=1 log I i=1 π i N (y n | μ i , Σ i ) + k log(N )
The estimated number of classes is the estimated number of 3D signals and the corresponding centers of ellipsoids are then output as candidates centers.

The outputs of the first step are 3D coordinates of candidate centers and 2D blocks with a lesion candidate at the center extracted from the central slice of candidates. The block size M × M is pre-defined to cover biggest lesions.

Second step -CJO detection decision

In this step, a test statistic is calculated for each output 2D block using the channelized joint detection-estimation observer (CJO). The CJO was originally proposed for the detection of an SKS signal on a 2D image in our earlier study [START_REF] Zhang | Channelized model observer for the detection and estimation of signals with unknown amplitude, orientation, and size[END_REF], that can be summarized into a training stage and a test stage:

Training stage: The inputs are two sets of image blocks (with-signal and without-signal) whose ground truth is known for the CJO, including the SKS signal profile x α (α is the signal parameter vector). Note that the signal is always located at the block center. The outputs of this stage are the estimated channelized reference signal x 0 and a template w:

x 0 = A α x α = A α (U t 0 x α ) , ( 2
)
w = Σ b -1 x 0 , (3)
Σ b = 1 2 (g -g | H 0 ) (g -g | H 0 ) t | H 0 + 1 2 [(g -x α ) -(g -x α ) | H 1 ] • [(g -x α ) -(g -x α ) | H 1 ] t | H 1 ; (4) g = A α U t 0 g (5)
where the channel matrix U 0 reduces the data dimensionality, and A α maps the channelized parametric signal x α onto the channelized reference signal x 0 : A α x α = x 0 .

To search the optimal parameters in the channel domain without loss of accuracy, one possible channel design (U 0 and A α ) for the amplitude-orientation-scale-unknown case was proposed in [START_REF] Zhang | Channelized model observer for the detection and estimation of signals with unknown amplitude, orientation, and size[END_REF] and adopted in this study. Note that K orientation-steerable channels and J scale-shiftable channels were used to construct U 0 .

Test stage: w and x 0 are used to estimate the signal parameters, as well as to calculate the test statistic λ for each input test image block whose ground truth is unknown to the CJO. For a fixed number of equally spaced signal parameters, the maximum of the test statistics yielded by all the possible parameter combinations is chosen as the final test statistic:

λ = max α (λ α ) = max α ⎛ ⎜ ⎝ w t U (A α ) t 2 F A α g - 1 2 x 0 ⎞ ⎟ ⎠ . ( 6
)
Because the optimization in ( 6) takes place in the channel domain, the maximization of ( 6) is still relatively fast, even when done using a brute-force search.

RESULTS AND DISCUSSIONS

Experimental images

We collected 20 healthy subjects' MR brain images (3D axial stacks, T2 FLAIR sequence). Half of them were used for the CJO training stage, in which 1000 different 2D white matter (WM) regions were extracted (500 without signal and 500 with a simulated signal). The other 10 images were used as test images (for both human observers and the msPCJO), on which we added simulated 3D lesions as introduced in Section 2.2. The lesion locations were indicated manually to make sure that the lesions were added to the WM (quasirandomly); the number of lesions in each 3D image was different, but fewer than 15 (to shorten the experiment time).

Fig. 3. Human observer study GUI.

Human observer study

We conducted a multi-slice free-response experiment, in which 6 radiologists were asked to detect and locate simulated MS lesions by scrolling back and forth the test images. When detecting a suspicious lesion, they selected its central slice and clicked the lesion center on the central slice. Then a dialog box popped up to ask the confidence rate (0-100).

Radiologists 1 and 2 were MS experts, with respectively 21 and 10 years of experience; Radiologists 3-6 were not MS experts, with respectively 6, 3, 8 and 5 years of experience.

Results and discussions

We used the software "JAFROC 4.1" [START_REF] Chakraborty | New developments in observer performance methodology in medical imaging[END_REF] to evaluate the task performance, and conducted the DBM-MRMC significance test (integrated in the software) to examine JAFROC1 (figure of merite) values' difference between each pair of observers.

Radiologists' JAFROC1 values are shown in Table 1, pvalues in Table 2. Table 3 shows the msPCJO's JAFROC1 values w.r.t. different number of steerable channels K and number of scale-shiftable channels J (cf. Section 2.2.2). Table 4 shows the p-values for the msPCJO with different (K, J) settings and 6 radiologists. Note that p-values< 0.05 are highlighted in red, which means that there is a significant difference between two corresponding JAFROC1 values. [START_REF] Liang | Image browsing in slow medical liquid crystal displays[END_REF][START_REF] Kim | A comparison of planar versus volumetric numerical observers for detection task performance in whole-body PET imaging[END_REF] 0.7645 0.1057 [START_REF] Liang | Image browsing in slow medical liquid crystal displays[END_REF][START_REF] Lartizien | Volumetric model and human observer comparisons of tumor detection for whole-body positron emission tomography[END_REF] 0.7820 0.0991 [START_REF] Lartizien | Volumetric model and human observer comparisons of tumor detection for whole-body positron emission tomography[END_REF][START_REF] Liang | Image browsing in slow medical liquid crystal displays[END_REF] 0.7029 0.1336 [START_REF] Lartizien | Volumetric model and human observer comparisons of tumor detection for whole-body positron emission tomography[END_REF][START_REF] Kim | A comparison of planar versus volumetric numerical observers for detection task performance in whole-body PET imaging[END_REF] 0.7622 0.1094 [START_REF] Lartizien | Volumetric model and human observer comparisons of tumor detection for whole-body positron emission tomography[END_REF][START_REF] Lartizien | Volumetric model and human observer comparisons of tumor detection for whole-body positron emission tomography[END_REF] 0.7785 0.1033 [START_REF] Gifford | A comparison of human and model observers in multislice LROC studies[END_REF][START_REF] Liang | Image browsing in slow medical liquid crystal displays[END_REF] 0.7099 0.1324 [START_REF] Gifford | A comparison of human and model observers in multislice LROC studies[END_REF][START_REF] Kim | A comparison of planar versus volumetric numerical observers for detection task performance in whole-body PET imaging[END_REF] 0.7622 0.1102 [START_REF] Gifford | A comparison of human and model observers in multislice LROC studies[END_REF][START_REF] Lartizien | Volumetric model and human observer comparisons of tumor detection for whole-body positron emission tomography[END_REF] 0.7797 0.1019

We see that from Table 3 that the number of steerable channels K hardly influences the msPCJO's performances, while the number of scale-shiftable channels J does. Table 4 shows that there is no significant difference between Note that when (K, J) = (3, 5) (i.e. 3 × 5 = 15 channels are used), the msPCJO attains its maximum detectionlocalization task performance (around 0.78), which is nearest to the 6 radiologists' average level 0.8045. That means when we use the msPCJO in practice, we can approach average human performance with not much calculation burden (less channels, less calculation burden).

CONCLUSION AND PERSPECTIVES

In this paper, we have proposed a novel 3D MO called the msPCJO in order to predict human performance in a study that has some relevance for diagnostic tasks which better describe clinical practice (multi-slice, multiple-signal...). To our knowledge, we are the first ones to propose such an MO that can perform the joint detection-localization of multiple 3D signals with random amplitude, orientation, size, location and number of lesions. By integrating an HVS model, the msPCJO also considers human visual characteristics. The inherent training stage in the 2nd step of the msPCJO allows modelling the training and the decision making process of radiologists to some extent.

We also conducted a human observer study close to the clinical paradigm to get the detection-localization task performance of six radiologists when volumetric information was available. The results are promising for our studied modality and pathology: the msPCJO is close to the average human performance with reasonable calculation amount.

One limitation of the study is the small number of 3D test images (mainly because it is difficult to get healthy patients). More test images will be needed to enrich the experimental validation in the further. Another interesting aspect to investigate consist in the extension of the 2D CJO to 3D CJO. Then in the 2nd step of the msPCJO, the CJO will be applied on 3D cubes with a lesion candidate at the center, instead of 2D blocks only on the central slice.

Fig. 1 .

 1 Fig. 1. The 1st row shows contiguous slices extracted from a healthy 3D image. The 2nd row shows the corresponding slices plus simulated 3D signals. One signal is highlighted to illustrate its appearances on 2D slices: marked by a red square on its central slice and by red circles on the adjacent slices. The highlighted regions are zoomed in the 3rd row.

Fig. 2 .

 2 Fig. 2. Structure diagram of the multi-slice PCJO (msPCJO): the outputs of each step are highlighted in red.

Table 1 .

 1 JAFROC1 values of six radiologists.

		JAFROC1 FOM Standard Error
	Radiologist 1	0.9267	0.0161
	Radiologist 2	0.9106	0.0274
	Radiologist 3	0.7837	0.0636
	Radiologist 4	0.9192	0.0420
	Radiologist 5	0.6547	0.0852
	Radiologist 6	0.6320	0.0995

Table 2 .

 2 The p-values for each pair of the radiologists.

		Rad1 Rad2 Rad3 Rad4 Rad5 Rad6
	Rad 1	-	-	-	-	-	-
	Rad 2 0.93	-	-	-	-	-
	Rad 3 0.15	0.27	-	-	-	-
	Rad 4 0.95	0.93	0.26	-	-	-
	Rad 5 0.03	0.04	0.29	0.03	-	-
	Rad 6 0.02	0.02	0.21	0.02	0.85	-

Table 3 .

 3 The JAFROC1 values of the msPCJO.

	(K,J) JAFROC1 FOM Standard Error
	(3,3)	0.7076	0.1320

Table 4 .

 4 The p-values for the msPCJO vs. 6 radiologists.

	(K,J) Rad1 Rad2 Rad3 Rad4 Rad5 Rad6
	(3,3) 0.02	0.02	0.40	0.02	0.56	0.40
	(3,4) 0.07	0.09	0.83	0.09	0.22	0.14
	(3,5) 0.11	0.14	0.98	0.13	0.16	0.10
	(5,3) 0.01	0.02	0.37	0.02	0.59	0.43
	(5,4) 0.07	0.10	0.81	0.08	0.23	0.15
	(5,5) 0.10	0.13	0.95	0.12	0.17	0.11
	(7,3) 0.02	0.04	0.41	0.02	0.54	0.39
	(7,4) 0.07	0.10	0.81	0.08	0.23	0.15
	(7,5) 0.10	0.15	0.96	0.12	0.17	0.10
	the msPCJO's JAFROC1 values and those of 6 radiologists
	(p > 0.05) when J > 3.


Anthony Fiche
email: anthony.fiche@ensta-bretagne.fr

Jean-Christophe Cexus
email: jean-christophe.cexus@ensta-bretagne.fr

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

Ali Khenchaf
email: ali.khenchaf@ensta-bretagne.fr

Features modeling with an α-stable distribution: application to pattern recognition based on continuous belief functions

Keywords: Gaussian and α-stable model, unimodal features, continuous belief functions, supervised pattern recognition, Kolmogorov-Smirnov test, Bayesian approach

published or not. The documents may come

Features modeling with an α-stable distribution: Application to pattern recognition based on continuous belief functions

Introduction

The choice of a model has an important role in the problem of estimation. For example, the Gaussian model is a very efficient model which fits data in many applications as it is very simple to use and saves computation time. However, as is the case for all distribution models, Gaussian laws have some weaknesses and results can end-up being skewed. Indeed, as the Gaussian probability density function (pdf) is symmetrical, it is not valid when the pdf is not symmetrical. It is therefore difficult to choose the right model which fits data for each application. The Gaussian distribution belongs to a family of distributions called stable distributions. This family of distributions allows the representation of heavy tails and skewness. A distribution is said to have a heavy tail if the tail decays slower than the tail of the Gaussian distribution. Therefore, the property of skewness means that it is impossible to find a mode where probability density function is symmetrical. The main property of stable laws introduced by Lévy [START_REF] Lévy | Théorie des erreurs : La loi de Gauss et les lois exceptionnelles[END_REF] is that the sum of two independent stable random variables gives a stable random variable. α-stable distributions are used in different fields of research such as radar [START_REF] Achim | SAR image denoising via Bayesian wavelet shrinkage based on heavy-tailed modeling[END_REF][START_REF] Banerjee | Adaptive target detection in foliage-penetrating SAR images using alpha-stable models[END_REF], image processing [START_REF] Kuruoglu | Skewed α-stable distributions for modelling textures[END_REF] or finance [START_REF] Mcculloch | Financial applications of stable distributions, Handbook of statistics[END_REF][START_REF] Fama | The behavior of stock-market prices[END_REF], ... The aim of this study is to show the interest in fitting data with an α-stable distribution. In [START_REF] Williams | Bayesian Data Fusion of Multi-View Synthetic Aperture Sonar Imagery for Seabed Classification[END_REF], the author proposes to characterize the sea floor from a vector of features modeled by Gaussian mixture models (GMMs) using an Autonomous Underwater Vehicle (AUV). However, the problem with the Bayesian approach is the difficulty in considering the uncertainty of data. We therefore favored the use of an approach based on the theory of belief functions [START_REF] Dempster | Upper and Lower probabilities induced by a multivalued mapping[END_REF][START_REF] Shafer | A mathematical theory of evidence[END_REF]. In [START_REF] Fiche | Bayesian approach and continuous belief functions for classification[END_REF], the authors compared a Bayesian and belief classifier where data from sensors are modeled using GMMs estimated via an Expectation-Maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. This work has been extended to data modeled by α-stable mixture models [START_REF] Fiche | Influence de l'estimation des paramètres de texture pour la classification de données complexes[END_REF]. However, it is difficult to choose between these two models because the results are roughly the same. This paper raises two problems. Firstly, it is necessary to work on a data set where features are modeled by α-stable pdfs. The second problem is to know how to apply the theory of belief functions when data from sensors are modeled by α-stable distributions. This point has been dealt with in [START_REF] Fiche | Continuous belief functions and α-stable distributions[END_REF].

The recent characterization of uncertain environments has taken an important place in several fields of application such as the SOund Navigation And Ranging (SONAR). SONAR has therefore been used for the detection of underwater mines [START_REF] Quidu | Mine classification based on raw sonar data: an approach combining Fourier descriptors, statistical models and genetic algorithms[END_REF]. In [START_REF] Laanaya | Multi-view fusion based on belief functions for seabed recognition[END_REF], the author developed techniques to perform automatic classification of sediments on the seabed from sonar images. In [START_REF] Fiche | Bayesian approach and continuous belief functions for classification[END_REF], the authors classified sonar images by extracting features and modeling them with GMMs. In [START_REF] Pouliquen | Sea-bed identification using echo-sounder signals[END_REF], the authors characterized the sea floor using data from a mono-beam echo-sounder. A data set represents an echo signal amplitude according to time. They compare the time envelope of echo signal amplitude with a set of theoretical reference curves. In this paper, we build a classifier based on the theory of belief functions where a vector of features extracted from a mono-beam echo-sounder is modeled by a Gaussian and an α-stable distribution. The feature pdfs have only one mode. We finally show that it would be interesting to fit features with an α-stable model.

This paper is organized as follows: we first introduce the definition of stability, the method to construct α-stable probability density functions and the methods of estimation in Section 2. The definitions in the multivariate case are also presented. We introduce the notion of belief functions in discrete and continuous cases in Section 3. Belief functions are then calculated in the particular case of α-stable distributions and a belief classifier is constructed in Section 4. We finally classify data by modeling data with Gaussian and α-stable distributions before the generated data are classified by comparing the theory of belief functions with Gaussian and α-stable distributions in Section 5. We compare the results obtained with the theory of belief functions using a Bayesian approach.

The α-stable distributions

Gauss introduced the probability density function called "Gaussian distribution" (in honor of Gauss) in his study of astronomy. Laplace and Poisson developed the theory of characteristic function by calculating the analytical expression of Fourier transform of a probability density function. Laplace found that the Fourier transform of a Gaussian law is also a Gaussian law. Cauchy tried to calculate the Fourier transform of a "generalized Gaussian" function with the expression f n (x) = 1 π +∞ 0 exp (-ct n ) cos (tx)dt, but he did not solve the problem. When the integer n is a real α, we define the family of α-stable distributions. However, Cauchy did not know at that time if he had defined a probability density function. With the results of Polya and Bernstein, Janicki and Weron [START_REF] Janicki | Simulation and chaotic behavior of α-stable stochastic processes[END_REF] demonstrated that the family of α-stable laws are probability density functions. The mathematician Paul Lévy studied the central limit theorem and showed, with the constraint of infinite variance, that the limit law is a stable law [START_REF] Lévy | Théorie des erreurs : La loi de Gauss et les lois exceptionnelles[END_REF]. Motivated by this property, Lévy calculated the Fourier transform for all α-stable distributions. These stable laws which satisfy the generalized central limit theorem are interesting as they allow the modeling of data as impulsive noise when the Gaussian model is not valid. Another property of these laws is the ability to model heavy tails and skewness. In the following, we define the notions of stability, characteristic function, the way to build a pdf from characteristic function, the different estimators of α-stable distributions and the extension of α-stable distribution to the multivariate case.

Definition of stability

Stability is when the sum of two independent random variables which follow α laws, gives an α law. Mathematically, this definition means the following: a random variable X is stable, denoted X ∼ S α (β, γ, δ), if for all (a,b) ∈ R + × R + , there are c ∈ R + and d ∈ R such that:

aX 1 + bX 2 = cX + d, (1)
with X 1 and X 2 two independent α-stable random variables which follow the same distribution as X. If Equation (1) defines the notion of stability, it does not give any indication as to how to parameterize an α-stable distribution. We therefore prefer to use the definition given by characteristic function to refer to an α-stable distribution.

Characteristic function

Several equivalent definitions have been suggested in the literature to parameterize an α-stable distribution from its characteristic function [START_REF] Taqqu | Stable non-gaussian random processes[END_REF][START_REF] Zolotarev | One-dimensional stable distributions[END_REF]. Zolotarev [START_REF] Zolotarev | One-dimensional stable distributions[END_REF] proposed the following:

φ S α (β,γ,δ) (t) =            exp( jtδ -|γt| α [1 + jβ tan( πα 2 )sign(t)(|t| 1-α -1)]) if α 1, exp( jtδ -|γt|[1 + jβ 2 π sign(t) log |t|]) if α = 1, (2)
where each feature has specific values: [START_REF] Lévy | Théorie des erreurs : La loi de Gauss et les lois exceptionnelles[END_REF] is the skewness parameter.

• α ∈]0, 2] is the characteristic exponent. • β ∈ [-1,
• γ ∈ R + * represents the scale parameter.

• δ ∈ R is the location parameter.

The advantage of this parameterization compared to [START_REF] Taqqu | Stable non-gaussian random processes[END_REF] is that the values of the characterization and probability density functions are continuous for all parameters. In fact, the parameterization defined by [START_REF] Taqqu | Stable non-gaussian random processes[END_REF] is discontinuous when α = 1 and β = 0.

The probability density function

The representation of an α-stable pdf, denoted f S α (β,γ,δ) , is obtained by calculating the Fourier transform of its characteristic function:

f S α (β,γ,δ) (x) = ∞ -∞ φ S α (β,γ,δ) (t) exp(-jtx)dt.
(

However, this definition is problematic for two reasons: while the integrand function is complex, its bounds are infinite. Nolan [START_REF] Nolan | Numerical calculation of stable densities and distribution functions[END_REF] therefore proposed a way to represent normalized α-stable distributions (i.e. γ = 1 and δ = 0). The main idea of Nolan [START_REF] Nolan | Numerical calculation of stable densities and distribution functions[END_REF] is to use variable modifications so that the integral has finite bounds. Each parameter has an influence on the shape of the f S α (β,γ,δ) . The curve has a large peak when α is near 0 and a Gaussian shape when α = 2 (Figure 1(a)). The shape of the distribution skews to the left if β = 1, to the right if β = -1 (Figure 1(b)) while the distribution is symmetrical when β = 0. Finally, the scale parameter enlarges or compresses the shape of the distribution (Figure 1(c)) and the location parameter leads to the translation of the mode of the f S α (β,γ,δ) (Figure 1(d)).

An overview of α-stable estimators

The estimators of α-stable distributions are decomposed into three families:

• the sample quantile methods [START_REF] Fama | Parameter estimates for symmetric stable distributions[END_REF][START_REF] Mcculloch | Simple consistent estimators of the parameters of stable laws[END_REF] • the sample characteristic function methods [START_REF] Press | Applied multivariate analysis[END_REF][START_REF] Press | Estimation in univariate and multivariate stable distributions[END_REF][START_REF] Koutrouvelis | Regression-type estimation of the parameters of stable laws[END_REF][START_REF] Koutrouvelis | An iterative procedure for the estimation of the parameters of stable laws[END_REF] • the Maximum Likelihood Estimation [START_REF] Brorsen | Maximum likelihood estimates of symmetric stable distribution parameters[END_REF][START_REF] Mcculloch | Linear regression with stable disturbances[END_REF][START_REF] Dumouchel | Stable distributions in statistical inference[END_REF][START_REF] Nolan | Maximum likelihood estimation and diagnostics for stable distributions[END_REF].

Fama and Roll [START_REF] Fama | Parameter estimates for symmetric stable distributions[END_REF] developed a method based on quantiles. However, the algorithm proposed by Fama and Roll suffers from a small asymptotic bias in α and γ and restrictions on α ∈]0.6, 2] and β = 0. McCulloch [START_REF] Mcculloch | Simple consistent estimators of the parameters of stable laws[END_REF] extended the quantile method to the asymmetric case (i.e. β = 0). The McCulloch estimator is valid for α ∈]0.6, 2].

Press [START_REF] Press | Applied multivariate analysis[END_REF][START_REF] Press | Estimation in univariate and multivariate stable distributions[END_REF] proposed a method based on transformations of characteristic function. In [START_REF] Weron | Performance of the estimators of stable law parameters[END_REF], the author compared the performances of several estimators. For example, the Press method is efficient for specific values. Koutrouvelis [START_REF] Koutrouvelis | Regression-type estimation of the parameters of stable laws[END_REF] extended the Press method and proposed a regression method to estimate the parameters of α-stable distributions. He proposed a second version of his algorithm which is distinct that it is iterative [START_REF] Koutrouvelis | An iterative procedure for the estimation of the parameters of stable laws[END_REF]. In [START_REF] Akgiray | Estimation of stable-law parameters: A comparative study[END_REF], the authors proved that the method proposed by Koutrouvelis is better than both the quantile method and the Press method because it gives consistent and asymptotically unbiased estimates.

The Maximum Likelihood Estimation (MLE) was first studied in the symmetric case [START_REF] Brorsen | Maximum likelihood estimates of symmetric stable distribution parameters[END_REF][START_REF] Mcculloch | Linear regression with stable disturbances[END_REF]. Dumouchel [START_REF] Dumouchel | Stable distributions in statistical inference[END_REF] developed an approximate maximum likelihood method. The MLE was also developed in the asymmetric case.

Nolan [START_REF] Nolan | Maximum likelihood estimation and diagnostics for stable distributions[END_REF] extended the MLE in general case. The problem with the MLE is the calculation of the α-stable probability density function because there is no closed-form expression. Moreover, the computational algorithm is timeconsuming.

Consequently, we estimate a univariate α-stable distribution using the Koutrouvelis method [START_REF] Koutrouvelis | An iterative procedure for the estimation of the parameters of stable laws[END_REF].

Multivariate stable distributions

It is possible to extend the α-stable distributions to the multivariate case. A random vector X ∈ R is stable if for all a, b ∈ R + , there are c ∈ R + and D ∈ R d such that:

aX 1 + bX 2 = cX + D, (4)
where X 1 and X 2 are two independent and identically distributed random vectors which follow the same distribution as X.

The characteristic function of an multivariate α-stable distribution, denoted X ∼ S α,d (σ, δ), has the form:

φ S α,d (σ,δ) (t) =                exp - S d | < t, s > | α (1 -jsgn(< t, s >) tan( πα 2 )σ(ds) + j < δ, t > if α 1, exp - S d | < t, s > |(1 + j π 2 sgn(< t, s >) ln(< t, s >))σ(ds) + j < δ, t > if α = 1, (5)
with

• S d = {x ∈ R d |||x|| = 1} the d-dimensional unit sphere.
• σ(.) a finite Borel measure on S d .

• δ, t ∈ R d .
The expression for the characteristic function involves an integration over the unit sphere S d . The measure σ is called the spectral measure and δ is called the location parameter. The problem with a multivariate α-stable distribution is that characteristic function forms a non-parametric set. To avoid this problem, it is possible to consider a discrete spectral measure [START_REF] Modarres | A method for simulating stable random vectors[END_REF] which has the form:

σ(.) = K i=1 γ i δ s i (.), (6)
with γ i corresponding to weight and δ s i the Dirac measure in s i . For an α-stable distribution in R 2 with K mass points, the quantity s i = {cos(θ i ), sin(θ i )} with θ i = 2π(i-1) K . There are two methods to estimate an α-stable random vector:

• the PROJection method [START_REF] Mcculloch | Estimation of the bivariate stable spectral representation by the projection method[END_REF] (PROJ)

• the Empirical Characteristic Function method [START_REF] Nolan | Estimation of stable spectral measures[END_REF] (ECF) These two algorithms have the same performances in terms of estimation and computation time.

At this point, it is important to underline that the features extracted from mono-beam echo-sounders have the properties of heavy-tails and skewness. Consequently, we decided to estimate the data with an α-stable model. These data are however imprecise and uncertain: the imprecision and uncertainty of data can be linked to poor quality estimation of these data. To take these constraints into consideration, we used an uncertain theory called the theory of belief functions. Consequently, the goal of the next section is to present the theory of belief functions.

The belief functions

The final objective of this paper is to classify synthetic and real data using the theory of belief functions. Data obtained from sensors are generally imprecise and uncertain as noises can disrupt their acquisition. It is possible to consider these constraints by using the theory of belief functions. We will first develop the theory of belief functions within a discrete framework before characterizing it in real numbers.

Discrete belief functions

This section outlines basic tools in relation to the theory of belief functions.

Definitions

Discrete belief functions were introduced by Dempster [START_REF] Dempster | Upper and Lower probabilities induced by a multivalued mapping[END_REF], and formalized by Shafer [START_REF] Shafer | A mathematical theory of evidence[END_REF] where he considers a discrete set of n exclusive events C i called the frame of discernment:

Θ = {C 1 , . . . , C n }. (7)
Θ can be interpreted as all the assumptions to a problem. Belief functions are defined as 2 Θ onto [0,1]. The objective of discrete belief functions is to attribute a weight of belief to each element A ∈ 2 Θ . Belief functions must follow the normalization:

A⊆Θ m Θ (A) = 1, (8)
where m Θ is called the basic belief assignment (bba).

A focal element is a subset of A where m Θ (A) > 0 and several functions in one-to-one correspondence are built from bba:

bel Θ (A) = B⊆A,B ∅ m Θ (B), ( 9
)
pl Θ (A) = A∩B ∅ m Θ (B), (10)
q Θ (A) = B⊂Θ,B⊇A m Θ (B). (11)
The credibility function of A called bel Θ (A) is all the elements B ⊆ A which believe partially in A. This function can be interpreted as a minimum of belief in A. On the contrary, a plausibility function pl Θ (A) illustrates the maximum belief in A while commonality function q Θ (A) represents the sum of bba allocated to the superset of A. This function is very useful as shown below.

Combination rule

We consider M different experts who give mass m Θ i (i = 1, . . . , M) on each element A ⊆ Θ. Alternatively, it is possible to combine them using combination rules. There are several combination rules [START_REF] Sentz | Combination of evidence in Dempster-Shafer theory[END_REF] in the literature which differently address conflicts between sources. The most common rule is the conjunctive combination [START_REF] Ph | The combination of evidence in the transferable belief model[END_REF] where the resultant mass of A is obtained by:

m Θ (A) = B 1 ∩...B n =A ∅ M i=1 m Θ i (B i ) , ∀A ∈ 2 Θ . ( 12
)
The mass of the empty set is given by:

m Θ (∅) = B 1 ∩...B n =∅ M i=1 m Θ i (B i ) , ∀A ∈ 2 Θ . ( 13
)
This rule allows us to stay in the open world. However, this is not practical as calculations are difficult. It is possible to calculate resultant mass with communality functions where each mass m Θ i (i = 1, . . . , M) must be converted into its communality function to calculate the resultant communality function as follows:

q Θ (A) = M i=1 q Θ i (A). ( 14
)
The final mass is obtained by carrying out the inverse operation of Equation ( 11) [START_REF] Shafer | A mathematical theory of evidence[END_REF].

Pignistic probability

To make a decision on Θ, several operators exist such that maximum credibility or maximum plausibility with pignistic probability being the most commonly used operator [START_REF] Ph | Constructing the pignistic probability function in a context of uncertainty[END_REF]. This name comes from pignus, a bet, in Latin. This operator approaches the pair (bel,pl) by uniformly sharing a mass of focal elements on each singleton C i . This operator is defined by:

betP(C i ) = A⊂Θ,C i ∈A m Θ (A) |A|(1 -m Θ (∅)) , (15)
where |A| represents the cardinality of A. We choose the decision C i by evaluating max 1≤k≤n betP(C k ).

Continuous belief functions

The basic description of continuous belief functions was accomplished by Shafer [START_REF] Shafer | A mathematical theory of evidence[END_REF], then by Nguyen [START_REF] Nguyen | On random sets and belief functions[END_REF] and Strat [START_REF] Strat | Continuous belief functions for evidential reasoning[END_REF]. Recently, Smets [START_REF] Ph | Belief functions on real numbers[END_REF] extended the definition of belief functions to the set of reals R = R ∪ {-∞, +∞} and masses are only attributed to intervals of R.

Definitions

Let us consider I = {[x, y], (x, y], [x, y), (x, y); x, y ∈ R} as a set of closed, half-opened and opened intervals of R. Focal elements are closed intervals of R. The quantity m I (x, y) are basic belief densities linked to a specific pdf. If x > y, then m I (x, y) = 0. With these definitions, it is possible to define the same functions as in the discrete case. The interval [a, b] being a set of R with a ≤ b, the previous functions can be defined as follows:

bel R ([a, b]) = x=b x=a y=b y=x m I (x, y)dydx, ( 16
) pl R ([a, b]) = x=b x=-∞ y=+∞ y=max(a,x) m I (x, y)dydx, ( 17
) q R ([a, b]) = x=a x=-∞ y=+∞ y=b m I (x, y)dydx. ( 18
)

Pignistic probability

The definition of pignistic probability for a < b is:

Bet f ([a, b]) = x=+∞ x=-∞ y=+∞ y=x |[a, b] ∩ [x, y]| |[x, y]| m I (x, y)dxdy. ( 19
)
It is possible to calculate pignistic probabilities to have basic belief densities. However, many basic belief densities exist for the same pignistic probability. To resolve this issue, we can use the consonant basic belief density. A basic belief density is said to be "consonant" when focal elements are nested. Focal elements I u can be labeled as an index u such that I u ⊆ I ′ u with u ′ > u. This definition is used to apply the least commitment principle, which consists in choosing the least informative belief function when a belief function is not totally defined and is only known to belong to a family of functions. The least commitment principle relies on an order relation between belief functions in order to determine if a belief function is more or less committed than another. For example, it is possible to define on order based on the commonality function:

(∀A ⊆ Θ, q Θ 1 (A) ≤ q Θ 2 (A)) ⇔ (m Θ 1 ⊆ q m Θ 2 ). ( 20
)
The mass function m Θ 2 is less committed than m Θ 1 according to the commonality function. The function Bet f can be induced by a set of isopignistic belief functions Biso(Bet f ). Many papers [START_REF] Ph | Belief functions on real numbers[END_REF][START_REF] Ristic | Target classification approach based on the belief function theory[END_REF][START_REF] Caron | Least Committed basic belief density induced by a multivariate Gaussian pdf[END_REF] deal with the particular case of continuous belief functions with nested focal elements. For example, Smets [START_REF] Ph | Belief functions on real numbers[END_REF] proved that the least committed basic belief assignment m R for the commonality ordering attributed to an interval I = [x, y] with y > x related to a bell-shaped1 pignistic probability function is determined by2 :

m R ([x, y]) = θ(y)δ d (x -ζ(y)), (21)
θ(y) = (ζ(y) -y) dBet f (y) dy . ( 22
)
The resultant basic belief assignment m R is consonant and belongs to the set Biso(Bet f ). However, it is difficult to build belief functions in the particular case of multimodal pdfs because the frame of discernment has connected sets.

In [START_REF] Doré | Constructing of a consonant belief function induced by a multimodal probability density function[END_REF][START_REF] Vannobel | Continuous belief functions: singletons plausibility function in conjunctive and disjunctive combination operations of consonant bbds[END_REF], the authors propose a way to build belief functions with connected sets by using a credal measure and an index function.

Credal measure and index function

In [START_REF] Doré | Constructing of a consonant belief function induced by a multimodal probability density function[END_REF], the authors propose a way to calculate belief functions from any probability density function. They use an index function f and a specific index space I to scan the set of focal elements F :

f I : I → F , (23)
y -→ f I (y). ( 24
)
The authors introduce a positive measure µ Ω such that I dµ Ω (y) ≤ 1 describes unconnected sets. The pair ( f I , µ Ω ) defines a belief function. For all A ∈ Ω, they define subsets which belong to the Borel set:

F ⊆A = {y ∈ I| f I (y) ⊆ A}, ( 25
)
F ∩A = {y ∈ I| f I (y) ∩ A ∅}, (26)
F ⊇A = {y ∈ I| f I (y) ⊇ A}. ( 27
)
From these definitions, they compute belief functions:

bel

Ω (A) = F ⊆A dµ Ω (y), ( 28
)
pl

Ω (A) = F ∩A dµ Ω (y), (29)
q Ω (A) =

F ⊇A dµ Ω (y). ( 30
)
They continue their study by considering consonant belief functions. The set of focal elements F must be ordered from the operator ⊆. They define an index function f from R + to F such that:

y ≥ x =⇒ f (y) ⊆ f (x). ( 31
)
They generate consonant sets by using a continuous function g from R d to I = [0, α max [. The α-cuts are the set:

f I cs = {x ∈ R d |g(x) ≥ α}. ( 32
)
They finally define the index function:

f I cs : I = [0, α max ] → { f I cs (α)|α ∈ I}, (33) α -→ f I cs (α). ( 34
)
The information available is the conditional pignistic density Bet f . However, many basic belief densities exist for the same pignistic probability Bet f . The least commitment principle allows the least informative basic belief density to be chosen, where the focal elements are the α-cuts of Bet f such that3 :

dµ Ω (y)(α) = λ( f I cs (α))dλ(α). ( 35
)

Continuous belief functions and α-stable distribution

In this section, we model data distributions as a single α-stable distribution. We must however introduce the notion of plausibility function for an α-stable distribution. We first describe how to calculate the plausibility function knowing the pdf in R before we extend to in R d . We finally explain how we construct our belief classifier.

Link between pignistic probability function and plausibility function in R

The information available is the conditional pignistic density Bet f

[C i ] with C i ∈ θ. The function Bet f [C i
] is supposed to be bell-shaped for all α-stable distributions (proved by Yamazato [START_REF] Yamazato | Unimodality of infinitely divisible distribution functions of class L[END_REF]).

The plausibility function from a mass m R is obtained by an integral of Equation ( 22) between [x, +∞[:

pl R [C i ](I) = +∞ x (ζ(t) -t) dBet f (t) dt dt. ( 36
)
By assuming that Bet f is symmetrical, an integration by parts can simplify Equation ( 36) :

pl R [C i ](I) = 2(x -µ)Bet f (x) + 2 +∞ x Bet f (t)dt. (37)
Now let us consider a particular case where symmetrical Bet f is an α-stable distribution (the parameter β = 0). We already know that Bet f (x) = f S α (β,γ,δ) (x). We can calculate

+∞ x

Bet f (t)dt by using the Chasles' theorem:

+∞ -∞ f S α (β,γ,δ) (t)dt = x -∞ f S α (β,γ,δ) (t)dt + +∞ x f S α (β,γ,δ) (t)dt. ( 38
)
By definition, a pd f has the quantity

+∞ -∞ f S α (β,γ,δ) (t)dt = 1 and
x -∞ f S α (β,γ,δ) (t)dt represents the definition of the α-stable cumulative density function F S α (β,γ,δ) .

Consequently, Equation (37) can be simplified:

pl R [C i ](I) = 2(x -µ) f S α (β,γ,δ) (x) + 2(1 -F S α (β,γ,δ) (x)). ( 39
)
The plausibility function related to an interval I = [x, y] can also be seen as the area defined under the α cut -cut such that

α cut = Bet f (x). The notation pl[C i ](x) is equivalent to pl[C i ](I).
Now us let consider an asymmetric α-stable probability density function. We must proceed numerically to calculate plausibility function at point x 1 > µ, with µ the mode of the probability density function (Figure 2). The plausibility function related to an interval I 1 = [x 1 , y 1 ] is defined by the area defined under the α cut -cut such that

α cut = Bet f (x 1 ): pl R [C i ](I 1 ) = x 1 -∞ Bet f (t)dt + (y 1 -x 1 )Bet f (x 1 ) + +∞ y 1 Bet f (t)dt. ( 40
)
By definition, a pd f has the quantity

+∞ y 1 f S α (β,γ,δ) (t)dt = 1 -F S α (β,γ,δ) (y 1 ) and x 1 -∞ f S α (β,γ,δ) (t)dt = F S α (β,γ,δ) (x 1 )
. In general, we know only one point y 1 . We estimate numerically x 1 such that f S α (β,γ,δ) (y 1 ) = f S α (β,γ,δ) (x 1 ). Finally, the plausibility function related to the interval I 1 is:

pl R [C i ](I 1 ) = 1 + F S α (β,γ,δ) (x 1 ) -F(y 1 ) + (y 1 -x 1 ) f (x 1 ). ( 41
)
In practice, we base the classification on several features defined for different classes Θ = {C 1 , .., C n }. For example, the Haralick parameters "contrast" and "homogeneity" have different values for the classes rock, silt and sand. The features are modeled by a probability density function because the features have continuous values. We can calculate a plausibility function related to its probability density functions by using the least commitment principle. Several plausibility functions associated to the same feature can be combined by using the generalized Bayes theorem [START_REF] Ph | Belief functions: The disjunctive rule of combination and the generalized Bayesian theorem[END_REF][START_REF] Delmotte | Target identification based on the transferable belief model interpretation of Dempster-Shafer model[END_REF] to calculate mass functions allocated to A of an interval I:

m R [x](A) = C j ∈A pl j (x) C j ∈A c (1 -pl j (x)). ( 42
)
For several features, it is possible to combine the mass functions with combination rules.

To validate our approach, we classify planes using kinematic data as in [START_REF] Ristic | Target classification approach based on the belief function theory[END_REF][START_REF] Caron | Least Committed basic belief density induced by a multivariate Gaussian pdf[END_REF] and compare the decision with the approach of Caron et al. [START_REF] Caron | Least Committed basic belief density induced by a multivariate Gaussian pdf[END_REF]. We associate a Gaussian probability density function for each speed's target (Figure 3):

• Commercial defined by the probability density function of speed S 2 (0, 8, 722.5).

• Bomber defined by the probability density function of speed S 2 (0, 7, 690).

• Fighter defined by the probability density function of speed S 2 (0, 10, 730).

We can observe that the decision is the same with the both approaches in the particular case of Gaussian probability density functions.

Link between pignistic probability function and plausibility function to

R d It is possible to extend plausibility function in R d .
In [START_REF] Caron | Least Committed basic belief density induced by a multivariate Gaussian pdf[END_REF], the authors calculate plausibility function in the Gaussian pdf situation of mode δ and matrix of covariance Σ. Mass function is built in such a way that isoprobability surfaces S i with 1 ≤ i ≤ n are focal elements. In R 2 , isoprobability points of a multivariate Gaussian pdf of mean δ and covariance matrix Σ are ellipses. In dimension d, the authors defined focal elements as the nested sets HV α enclosed by the isoprobability hyperconics HC α = {x ∈ R d |(xδ) t Σ -1 (xδ)}. They obtain bbd by applying the least commitment principle:

m R d (HV α ) = α d+2 2 -1 2 d+2 2 Γ( d+2 2 ) exp (- 1 2 α) with α ≥ 0. ( 43
)
Equation ( 43) defines a χ 2 distribution with d + 2 degrees of freedom. The plausibility function at point x belonging to a surface S i corresponds to the volume delimited by S i (we give an example in Figure 4 for an α-stable pdf). Consequently, the plausibility function at point x is defined by:

pl R d (x ∈ R d ) = α=+∞ α=(x-δ) t Σ -1 (x-δ) α d+2 2 -1 2 d+2 2 Γ( d+2 2 ) exp (- 1 2 α)dα. (44)
Equation ( 44) can be simplified by:

pl R d (x ∈ R d ) = 1 -F d+2 ((x -δ)(Σ) -1 (x -δ)). (45)
The function F d+2 is a cumulative density function of the χ 2 distribution with d + 2 degrees of liberty (d is the dimension of vector x).

The authors also calculate plausibility functions in the particular case of GMMs with their pdfs denoted by p GMM :

p GMM (x) = k=n k=1 w k N(x, δ k , Σ k ), (46)
where N(x, δ k , Σ k ) is the normal distribution of the k components with mean δ k and matrix of covariance Σ k and w k is the weight of each mixture. They assign belief to nested sets belonging to a component. Consequently, the plausibility function of the GMMs can be seen as a weighted sum of plausibility functions defined by each component of the mixture:

pl R d (x ∈ R d ) = 1 - k=n k=1 w k F d+2 ((x -δ k )(Σ k ) -1 (x -δ k )). ( 47
)
We want to extend the calculation of plausibility function for any α-stable pdf. However, there is no closed-form expression for α-stable pdf. We use the approach developed by Doré et al. [START_REF] Doré | Constructing of a consonant belief function induced by a multimodal probability density function[END_REF] to build belief functions.

Belief classifier

Let us consider a data set with N samples from d sensors. For example, each feature of a vector x ∈ R d can be seen as a piece of information from a sensor. The classification is divided into two steps. N × p (with p ∈]0, 1[) samples are first picked out randomly from the data set for the learning base, noted X, such that:

X =              x 1 1 . . . x N×p 1 . . . . . . . . . x 1 d . . . x N×p d             
All columns of X belong to a class C i with 1 ≤ i ≤ n. Probability density functions (Gaussian or α-stable models) are estimated from samples belonging to classes. The rest of the samples is used for the test base (N × (1p) vectors). We use a validation test to determine if samples belong to an estimated model. If the test is not valid, we stop the classification, otherwise we continue the classification (Figure 5). The classification step diagram in one and d dimensions is shown in Figure 6.

Plausibility functions knowing classes for x are calculated either from their probability density functions by using Equation ( 41) for an unsymmetric α-stable probability density function or by using Equation ( 45) for a Gaussian probability density function. We obtain d mass functions (one for each feature) at point x with the generalized Bayes theorem [START_REF] Ph | Belief functions: The disjunctive rule of combination and the generalized Bayesian theorem[END_REF][START_REF] Delmotte | Target identification based on the transferable belief model interpretation of Dempster-Shafer model[END_REF] (Equation ( 42)). These d mass functions are combined by a combination rule to obtain a single mass function (section 3.1.2). Finally, the mass function is transformed into pignistic probability (Equation ( 19)) to make the decision.

It is also possible to work with a vector x of d dimensions. For an α-stable probability density function, plausibility functions are calculated for each feature by using the approach of Doré et al [START_REF] Doré | Constructing of a consonant belief function induced by a multimodal probability density function[END_REF] (section 3.3). For a Gaussian probability density function, we use Equation [START_REF] Vannobel | Continuous belief functions: singletons plausibility function in conjunctive and disjunctive combination operations of consonant bbds[END_REF] to calculate plausibilty functions. We calculate one mass function at point x by using the generalized Bayes theorem (Equation ( 42)). There is no combination step because we are calculating a single mass function. We use Equation [START_REF] Zolotarev | One-dimensional stable distributions[END_REF] to transform the mass function into pignistic probabilities. The decision is chosen by using the maximum number of pignistic probabilities.

Application to pattern recognition

The aim is to build a belief classifier and to perform a classification of synthetic and real data by estimating features using a Gaussian and α-stable pdf. In [START_REF] Fiche | Influence de l'estimation des paramètres de texture pour la classification de données complexes[END_REF], synthetic data are classified by modeling features using Gaussian and α-stable mixture models. By observing confidence intervals, the hypothesis of α-stable mixture models is significantly better than the hypothesis of Gaussian mixture models. However, when the number of Gaussian distributions increases, classification accuracies are significantly the same. Images from a side-scan sonar are automatically classified by extracting Haralick features [START_REF] Haralick | Textural features for image classification[END_REF][START_REF] Haralick | Statistical and structural approaches to texture[END_REF]. However, classification accuracies are roughly the same.

In this section, we limit algorithms with a vector of features in dimension d ≤ 2. Indeed, the generalization of α-stable pdf in dimension d > 2 increases CPU time because there is no closed-form expression. Consequently, we distinguish two cases during the classification step:

• The one dimension case: each dimension is considered as a feature.

• The two case: we considered a vector of two features.

We also use a statistical test called Kolmogorov-Smirnov test (K-S test) to evaluate the quality of each model. A Kolmogorov-Smirnov test (K-S test) with a significance level of 5 % is also used to determine if two datasets differ significantly. Let us assume two null hypotheses:

• H 1 : "Data follow a Gaussian distribution".

• H 2 : "Data follow an α-stable distribution".

The notation H 1 and H 2 are used for the rest of the paper. If the test rejects a null hypothesis H i , we stop the classification. Otherwise, if the test is valid, we then classify the samples belonging to the test base by using the belief classifier.

The results obtained with the belief functions are compared using a Bayesian approach. The prior probability p(C i ) is first calculated corresponding to the proportion of each class in the learning set. For each class, the application of Bayes theorem gives posterior probabilities:

p(C i /x) = p(x/C i )p(C i ) n u=1 p(x/C u )p(C u ) . ( 48
)
Finally, the decision is chosen by using the maximum nimber of posterior probabilities. We first present a classification of synthetic data, belonging to Gaussian distributions, by modeling features using a Gaussian distribution and α-stable pdf. We then compare the α-stable and Gaussian model with synthetic data belonging to α-stable distributions. The same approach is finally applied to real data from a mono-beam echo-sounder.

Application to synthetic data generated by Gaussian distributions

In this subsection, we show that the α-stable and Gaussian model have the same behavior during the classification when synthetic data are generated by Gaussian distributions. We simulated three Gaussian pdfs in R × R (c.f. Table 1). Each distribution holds 3000 samples. A mesh between [-4, 4] × [-4, 4] enables the level curves of pdf to be plotted (Figure 7).

Presentation of the data

µ V First Gaussian pdf (C 1 ) 2 3 1 1.5 1.5 3 Second Gaussian pdf (C 2 ) 1 1 3 1.5 1.5 1 Third Gaussian pdf (C 3 ) -1 1
1 0 0 3

Results

We generated data and split samples randomly into two sets: one third for the learning base and the rest for the test base. The learning base was used to estimate each mean and variance for the Gaussian distributions and each parameter α, β, γ and δ for the α-stable distributions [START_REF] Fama | The behavior of stock-market prices[END_REF][START_REF] Koutrouvelis | Regression-type estimation of the parameters of stable laws[END_REF][START_REF] Mcculloch | Simple consistent estimators of stable distribution parameters[END_REF]. We needed to make a mesh between [-4, 4] × [-4, 4] to calculate plausibility functions in R 2 with the α-stable model. Consequently, samples were chosen possessing features with dimensions within [-4, 4].

H

The one dimension case.

For each class C i , we consider each feature as a source of information. The estimated probability density functions can be seen as pignistic probability functions. For each pignistic probability function, we associate a mass function which will subsequently be combined with other mass functions.

The K-S test is valid for the two null hypotheses H 1 and H 2 (Table 2). We can assert that samples are drawn from the same distribution. The proposed method performs well in the particular case of synthetic data belonging to Gaussian distributions. However, the classification rates (Table 3) are very low because there is confusion between classes (Figure 8). The classification accuracy obtained with the Bayesian approach is significantly the same as the belief functions. Indeed, the Bayesian performs well provided that the probability density functions are wellestimated. The classification results depend on the prior probabilitiy estimation: if a prior probability of a class is underestimated, the confusion between classes increases.

We extended our study to a vector of two dimensions to try to improve the classification rate.

The two dimension case.

We generalized the study by considering a vector of two features. The learning base was used to estimate each mean and matrix of covariance for the Gaussian distributions (Figure 9) and each parameter α, σ, θ and δ for the α-stable distributions [START_REF] Nolan | Estimation of stable spectral measures[END_REF] (Figure 10). As in one dimension, the two null hypotheses H 1 and H 2 were verified (Table 4). The classification results (Table 3) are somewhat less than impressive because there is confusion between classes (between C 1 and C 2 cf. Figure 7). Furthermore, the work in two dimensions enables better results to be obtained than the combination. The Bayesian approach gives significantly the same results as the theory of belief functions. Indeed, the estimations of model and the estimation of prior probabilities are well-estimated.

We extended our study by considering synthetic data generated by α-stable distributions.

Bayesian approach Belief approach

Classification accuracies (%) 95 % confidence intervals Classification accuracies (%) 95 % confidence intervals

H

Application to synthetic data generated by α-stable distributions

In this subsection, we show the interest in modeling features in one and two dimensions with a single α-stable distribution. 4: Statistical values for a Kolmogorov-Smirnov test with a significance level of 5 % (p-value: the critical value to reject the null hypothesis, ksstat: the greatest discrepancy between the observed and expected cumulative frequencies).

α σ θ δ First α-stable pdf (C 1 ) 1.5 1 1 0 π 2 0 0 Second α-stable pdf (C 2 ) 1.5 1 1 0 π 2 2 1.4 Third α-stable pdf (C 3 ) 1.5 1 1 0 π 2 1 0.5
Table 5: Numerical values for the α-stable pdf.

Presentation of the data

Three classes of artificial data sets were generated from 2D α-stable distributions [START_REF] Modarres | A method for simulating stable random vectors[END_REF] (c.f. Table 5 and [START_REF] Nolan | Calculation of multidimensional stable densities[END_REF] for the significance of σ and θ) with 3000 samples. A mesh between [-4, 4] × [-4, 4] enables the level curves of pdf to be plotted (Figure 11).

Results

The one dimension case. We can assert that the estimation with the α-stable hypothesis is better than the Gaussian hypothesis (Figure 12). There is a difference between the mode of the Gaussian model and the mode of α-stable model. Indeed, the mean, calculated from samples, for the Gaussian model does not correspond to the mode when samples do not belong to a Gaussian law. Moreover, the values of pdfs are not the same compared to the real pdfs. 6: Statistical values for a Kolmogorov-Smirnov test with a significance level of 5 % (p-value: the critical value to reject the null hypothesis, ksstat: the greatest discrepancy between the observed and expected cumulative frequencies).

p-value ksstat

The K-S test is not satisfied for the null hypothesis H 1 for each class and each component. Consequently, we stopped the classification step for the Gaussian model. However, the K-S test is valid for the null hypothesis H 2 (Table 6). We can observe that the results obtained with the theory of belief functions are significantly the same as the results obtained using the Bayesian approach. This phenomenon can be explained by the fact that prior probabilities were well-estimated. Consequently, we changed voluntarily the values of prior probabilities. We fixed arbitrarily the prior probabilities: p(C 1 ) = 1/6, p(C 2 ) = 2/3 and p(C 3 ) = 1/6. The classification accuracy decreased to 53.10 % (Table 7). Indeed, the learning base and the test base can sometimes be significantly different. For example, the prior probability of class sand can be greater than the prior probabilities of classes rock and silt, although it does not represent the truth, and the classification rate of class sand can be favoured. This phenomenon introduces confusion between classes and the classification rate can decrease. Consequently, the learning step is an important step for the Bayesian approach and can have a bad effect on the classification rate whereas the theory of belief functions is less sensitive to the learning step. The two dimension case.

Bayesian approach

We considered a vector of two features. In Figure 13, we draw the iso-contour lines of data estimated with the Gaussian model and in Figure 14 the iso-contour lines of data estimated with the α-stable model As in one dimension, the null hypothesis H 1 is rejected at significance level 5 %, but the null hypothesis H 2 is valid (Table 8).

The classification results (Table 9) are somewhat less than impressive because there is confusion between classes. Furthermore, we obtain the same classification accuracies with the vector of two dimensions compared to the combination of two features. However, the classification decreases when the prior probabilities are not well-estimated (Table 7). We also compared these results with a GMM. It is well-known that the GMM can accurately model an arbitrary continuous distribution. This model estimation can be performed efficiently using the Expectation-Maximization (EM) algorithm. We compared classification accuracies by increasing the number of components. We drew the classification accuracy (Figure 15) in relation to the number of components. Let us assume the null hypothesis:

• G n : "Data follow a Gaussian mixture model with n components".

The K-S test rejects the null hypothesis G n for a number of components n < 4. The representations of iso-contour lines (Figure 16, Figure 17 and Figure 18) confirm the K-S test. Furthermore, the classification accuracies are lower than classification rates with a single α-stable pdf.

However, we can not increase indefinitely the number of components because a phenomenon of overfitting can appear (for example the representation of iso-contour by using a GMM with 5 components). Consequently, the number of 4 components seems to be a good compromise between the overfitting and good estimation. The classification accuracies are roughly the same between the single α-stable model and the GMM with 4 components. p-value ksstat C 1 0.0683 0.1003 C 1 0.3694 0.0688 C 3 0.1154 0.0932 Table 8: Statistical values for a Kolmogorov-Smirnov test with a significance level of 5 % (p-value: the critical value to reject the null hypothesis, ksstat: the greatest discrepancy between the observed and expected cumulative frequencies).

Application to real data from a mono-beam echo-sounder 5.3.1. Presentation of data

We used a data set from a mono-beam echo-sounder given by the Service Hydrographique et Océanique de la Marine (SHOM). Raw data represents an echo signal amplitude according to time. In [START_REF] Pouliquen | Sea-bed identification using echo-sounder signals[END_REF], the authors classify seven types of sea floor by comparing the time envelope of echo signal amplitude with a set of theoretical references curves. In our application, raw data (Figure 19) were processed using the Quester Tangent Corporation (QTC) software [START_REF] Caughey | Sea bottom classification from echo sounding data[END_REF] to

Conclusions

In this paper, we have shown the advantages in using the α-stable distributions to model data with heavy-tails. We considered the imprecision and uncertainty of data for modeling pdf and the theory of belief functions offers a way to model these constraints. We have examined the fundamental definitions of this theory and have shown a way to calculate plausibility functions when the knowledge of sensors is an α-stable pdf. We finally validated our model by making a comparison with the model examined by Caron et al. [START_REF] Caron | Least Committed basic belief density induced by a multivariate Gaussian pdf[END_REF] where the pdf is a Gaussian.

Synthetic and real data were finally classified using the theory of belief functions. We estimated the learning base, calculated plausibility functions, combined plausibility functions and calculated maximum pignistic probability to obtain the classification accuracy. The K-S test shows that synthetic data can be modeled by an α-stable model. The proposed method performs well for synthetic data. The consideration of two features improves the classification rate compared to the combination of two features for the Gaussian data. This study was then extended to a real application by using features extracted from a mono-beam echo-sounder. The Gaussian model is limited when modeling features compared to the α-stable model because the α-stable model has more degrees of freedom. To sum up, the interest in using an α-stable model is that data, in both one or two dimensions, can be modeled. The problem with the α-stable model is the computation time. Indeed, we proceed numerically to calculate plausibility functions (discretization for the bivariate α-stable pdf). Furthermore, it is difficult to generalize in dimension d > 2 without increasing CPU time.

We present preliminary results about classification problem by modeling features with α-stable probability density functions in dimension d ≤ 2. In future works, we plan to make experiments in the high-dimensional case with d ≥ 3. One solution in perspective to improve the computation time is to optimize the code in another language such that C/C++. Another perspective to this work is the use of other features from sensors to improve the classification accuracy, images from side-scan sonar.

  with x = ζ(y) satisfying Bet f (ζ(y)) = Bet f (y) and θ(y):

  ksstat p-value ksstat C 1 0.1509 0.0873 0.1509 0.0873 C 1 0.0490 0.0650 0.1997 0.0836 C 3 0.1399 0.0870 0.1665 0.0841 Table

  Influence of parameter α with β = 0, γ = 1 and δ = 0. Influence of parameter β with α = 1.5, γ = 1 and δ = 0. Influence of parameter γ with α = 1.5, β = 0 and δ = 0. Influence of parameter δ with α = 1.5, β = 0 and γ = 1.

Figure 1 :Figure 2 :

 12 Figure 1: Different α-stable pdf.

Figure 3 :

 3 Figure 3: Representation of probability density function of speed and the pignistic probabilities calculated from the Fiche et al. and Caron et al. [43] approaches.

Figure 4 :Figure 5 :

 45 Figure 4: Isoprobability surface of f .

Figure 6 :

 6 Figure 6: Scheme of classification step in one and d dimensions.

Figure 7 :

 7 Figure 7: Iso-contour lines of generated Gaussian pdf.

Figure 8 :

 8 Figure 8: Empirical pdf and its estimations (The first row corresponds to the first feature and the second row corresponds to the second feature).

Figure 9 :

 9 Figure 9: Iso-contour lines of data estimated with a Gaussian model.

Figure 10 :

 10 Figure 10: Iso-contour lines of data estimated with an α-stable model.

Figure 11 :

 11 Figure 11: Iso-contour lines of generated α-stable pdf.

Figure 12 :

 12 Figure 12: Empirical pdf and its estimations (The first row corresponds to the first feature and the second row corresponds to the second feature).

Figure 13 :

 13 Figure 13: Iso-contour lines of data estimated with a Gaussian model.

Figure 14 :

 14 Figure 14: Iso-contour lines of data estimated with an α-stable model.

  interval for the α-stable model Confidence intervals for the Gaussian model

Figure 15 :

 15 Figure 15: Classification rates according the number of components.

Figure 16 :

 16 Figure 16: Representation of iso-contour by modeling data with a mixture of Gaussian for the class C 1 (from left to right: 2, 3, 4 and 5 components).

Figure 17 :

 17 Figure 17: Representation of iso-contour by modeling data with a mixture of Gaussian for the class C 2 (from left to right: 2, 3, 4 and 5 components).

Figure 18 :Figure 19 :

 1819 Figure 18: Representation of iso-contour by modeling data with a mixture of Gaussian for the class C 3 (from left to right: 2, 3, 4 and 5 components).

Figure 20 :

 20 Figure 20: Running sample variances.

Figure 21 :

 21 Figure21: Empirical pdf and its estimations (The first row corresponds to the feature called "third quantile calculated on echo signal amplitude" and the second row corresponds to the feature called "25th quantile calculated on cumulative energy").

Figure 22 :

 22 Figure 22: Level curves of three empirical α-stable distributions with 10000 samples (The first column refers to the class rock, the second column refers to the class sand and the third column refers to the class silt).

Figure 23 :

 23 Figure 23: Level curves with an α-stable model estimation (The first column refers to the class rock, the second column refers to the class sand and the third column refers to the class silt).

Figure 24 :

 24 Figure 24: Level curves with a Gaussian model estimation (The first column refers to the class rock, the second column refers to the class sand and the third column refers to the class silt).

Table 1 :

 1 Numerical values for the Gaussian pdf.

Table 2 :

 2 Statistical values for a Kolmogorov-Smirnov test with a significance level of 5 % (p-value: the critical value to reject the null hypothesis, ksstat: the greatest discrepancy between the observed and expected cumulative frequencies).

			1	H 2
			p-value ksstat p-value ksstat
	C 1	first feature	0.7021 0.0542 0.9241 0.0422
		second feature 0.8799 0.0452 0.9241 0.0422
	C 2	first feature	0.8693 0.0464 0.3736 0.0712
		second feature 0.6855 0.0557 0.3736 0.0712
	C 3	first feature	0.8437 0.0464 0.4150 0.0667
		second feature 0.6602 0.0551 0.7866 0.0493

Table 3 :

 3 Classification accuracies and 95 % confidence intervals.

Table

	C 1	first feature	0.1758 0.0851
		second feature 0.7614 0.0517
	C 2	first feature	0.8487 0.0458
		second feature 0.4794 0.0630
	C 3	first feature	0.5505 0.0621
		second feature 0.2753 0.0776

Table 7 :

 7 Classification accuracies and 95 % confidence intervals with the prior probabilities p(C 1 ) = 1/6, p(C 2 ) = 2/3 and p(C 3 ) = 1/6.

i.e. the pdf is unimodal with a mode µ, continuous and strictly monotonous increasing (decreasing) at left (right) of the mode.

δ d refers to the Dirac's measure.

λ refers to the Lebesgue's measure.

Acknowledgements

The authors would like to thank the Service Hydrographique et Océanique de la Marine (SHOM) for the data and G. Le Chenadec for his advice concerning the data.

Belief approach

Classification accuracies (%) 95 % confidence intervals Classification accuracies (%) 95 % confidence intervals obtain some features, which were normalized between [0,1]. These data had already been used for the navigation of an Autonomous Underwater Vehicle (AUV) [START_REF] Martin | Autonomous Underwater Vehicle sensors fusion by the theory of belief functions for Rapid Environment Assessment[END_REF]. The frame of discernment is Θ = {silt, rock, sand}, with 4853 samples from silt, 6017 from rock and 7338 from sand. The selection of features plays an important role in classification because the classification accuracy can be different according to the features. In [START_REF] Karoui | Seabed segmentation using optimized statistics of sonar textures[END_REF], the author gave an overview of methods based on features to classify data and choose features which discriminate all classes. Some features calculated by the QTC software can be modeled by an α-stable pdf. There is a graphical test called converging variance test [START_REF] Granger | Infinite Variance" and Research Strategy in Time Series Analysis[END_REF] which enables us to determine if samples belong to an α-stable distribution. This test calculates the estimation of variances by varying the number of samples. If the estimation of variance diverges, the samples can belong to an α-stable distribution. We give an example of this test for the feature called "third quantile calculated on echo signal amplitude" of the class rock (Figure 20). We chose the features called the "third quantile calculated on echo signal amplitude" and the "25th quantile calculated on cumulative energy".

Results

We randomly select 5000 samples for the data set. Half the samples were used for the learning base and the rest for the test base. We first consider features in one dimension, i.e. two mass functions were calculated and combined to obtain a single mass function.

The one dimension case.

We can observe that the assumption of the α-stable model can easily accommodate the data compared to the Gaussian model (Figure 21).

The α-stable hypothesis is valid with the K-S test whereas the K-S test rejects the Gaussian hypothesis (Table 10). The approach using the theory of belief functions (classification accuracy of 82.68 %) give better results than the Bayesian approach (classification accuracy of 80.64 %) but not significantly. This phenomenon can be explained by the fact that the estimated models and prior probabilities are well-estimated.

The two dimensions case.

We extended the study by using a vector of two features. The level curves of empirical data (Figure 22) show there is confusion between the class silt and the class sand. The α-stable model (Figure 23) considers the confusion between classes compared to the Gaussian model (Figure 24