Lianmeng Jiao
email: lianmeng.jiao@utc.fr

Thierry Denoeux
email: thierry.denoeux@hds.utc.fr

Quan

Quan Pan
email: quanpan@nwpu.edu.cn

Fusion of Pairwise Nearest-Neighbor Classifiers Based on Pairwise

Keywords: pattern classification, nearest-neighbor classifier, pairwise-weighted distance metric, Dempster-Shafer theory

I. INTRODUCTION

The nearest-neighbor (NN) rule, first proposed by Fix and Hodges [START_REF] Fix | Discriminatory analysis, nonparametric discrimination: consistency properties[END_REF], is one of the most popular and successful pattern classification techniques. Given a set of 𝑁 labeled samples (or phototypes) 𝒯 = {(x (1) , 𝜔 (1) ), ⋅ ⋅ ⋅ , (x (𝑁 ) , 𝜔 (𝑁 ) )} with input vector x (𝑖) ∈ ℛ 𝐷 and class label 𝜔 (𝑖) ∈ Ω = {𝜔 1 , ⋅ ⋅ ⋅ , 𝜔 𝑀 }, the NN rule classifies a query pattern y ∈ ℛ 𝐷 to the class of its nearest neighbor in the training set 𝑇 . The basic rationale of the NN rule is both simple and intuitive: patterns close in feature space are likely to belong to the same class. The good behavior of the NN rule with unbounded number of prototypes is well known [START_REF] Cover | Nearest neighbor pattern classification[END_REF]. However, in many practical pattern classification applications only a small number of prototypes is available. Typically, in such a scarce-prototype situation, the ideal asymptotical behavior of NN classifier degrades dramatically [START_REF] Devroye | A Probabilistic Theory of Pattern Recognition[END_REF]. This motivates the growing interest in finding variants of the NN rule and adequate distance measures (or metrics) that help improve the NN classification performance in small data set situations.

As the core of the NN rule, the distance metric plays a crucial role in determining the classification performance. To overcome the limitations of the original Euclidean (L2) distance metric, a number of adaptive methods have recently been proposed to address the distance metric learning issue. According to the structure of the metric, these methods can be mainly divided into two categories: global distance metric learning [START_REF] Bar-Hillel | Learning distance functions using equivalence relations[END_REF], [START_REF] Zhang | Parametric distance metric learning with label information[END_REF], [START_REF] Eick | Using clustering to learn distance functions for supervised similarity assessment[END_REF], and local distance metric learning [START_REF] Paredes | A class-dependent weighted dissimilarity measure for nearest neighbor classication problems[END_REF], [START_REF]Learning weighted metrics to minimize nearest-neighbor classification error[END_REF], [START_REF] Wang | Improving nearest neighbor rule with a simple adaptive distance measure[END_REF], [START_REF] Jahromi | A method of learning weighted similarity function to improve the performance of nearest neighbor[END_REF], [START_REF] Jiao | An evidential k-nearest neighbor classification method with weighted attributes[END_REF]. The first approach learns the distance metric in a global sense, i.e., to share the same simple-weighted (SW) distance metric for all the prototypes:

𝑑 𝑆𝑊 (x, y) = ⎷ 𝐷 ∑ 𝑗=1 𝜆 2 𝑗 (𝑥 𝑗 -𝑦 𝑗 ) 2 , (1)
where, x is a prototype in the training set, y is a query pattern to be classified, and 𝜆 𝑗 is the weight of the 𝑗-th feature. Although the above global distance metric is intuitively appealing, it is too coarse as the feature weights of the distance metric are irrelevant with the class labels of the patterns. This issue becomes more severe when some features behave distinctly for different classes (for example, one feature may be more discriminative for some classes, but irrelevant for others). So, many methods have been developed to learn a distance metric in a local setting, i.e., the feature weights may be different for different patterns. The most representative one is the class-dependent weighted (CDW) distance metric [START_REF] Paredes | A class-dependent weighted dissimilarity measure for nearest neighbor classication problems[END_REF], [START_REF]Learning weighted metrics to minimize nearest-neighbor classification error[END_REF], which is related to the class index of the prototype:

𝑑 𝐶𝐷𝑊 (x, y) = ⎷ 𝐷 ∑ 𝑗=1 𝜆 2 𝑐,𝑗 (𝑥 𝑗 -𝑦 𝑗 ) 2 , ( 2
)
where 𝑐 is the class index of prototype x. Though the above CDW distance metric provides more freedom than the SW one, however, as illustrated in the following example, this distance metric is insufficient to reflect the local specifics in feature space for query patterns in different classes. Fig. 1 illustrates a simple three-class classification problem, where the data in each class are uniformly distributed. (x (1) , 𝐴), (x (2) , 𝐵) and (x (3) , 𝐶) are two-dimensional data points in training set 𝑇 . y 1 and y 2 are the query data to be classified. Considering the classification of data y 1 between Class A and Class B, when calculating the distance of y 1 to x (1) and x (2) , intuitively, to avoid classifying it as Class B mistakenly, feature X should be x (1) (9,6)

x (2) (11,9)

x (3) [START_REF] Suen | The combination of multiple classifiers by a neural network approach[END_REF][START_REF] Jiao | An evidential k-nearest neighbor classification method with weighted attributes[END_REF] y 1 [START_REF] Wang | Improving nearest neighbor rule with a simple adaptive distance measure[END_REF][START_REF] Wang | Improving nearest neighbor rule with a simple adaptive distance measure[END_REF] • y 2 [START_REF] Jiao | An evidential k-nearest neighbor classification method with weighted attributes[END_REF][START_REF] Jiao | An evidential k-nearest neighbor classification method with weighted attributes[END_REF] Fig. 1. A three-class classification example.

given a larger weight. However, when classifying data y 2 as Class B or Class C, feature Y should be given a larger weight to determine the distance of y 2 to x (2) and x (3) . In other words, the feature weights should be related to the labels of class pairs to be classified.

Motivated by the above consideration, in this paper we propose a pairwise-weighted (PW) distance metric related to the labels of the class pairs to be classified. Because only two classes are considered for each PW distance metric, the feature weights can be learnt in a more local way. Based on each PW distance metric, a pairwise NN (PNN) classifier can be designed to separate two classes. Then a polychotomous NN classification problem can be solved by fusing several PNN classifiers.

A variety of schemes has been proposed for deriving a combined decision from individual ones, such as majority voting [START_REF] Ruta | Classifier selection for majority voting[END_REF], Bayes combination [START_REF] Kuncheva | Decision templates for multiple classifier fusion: an experimental comparison[END_REF], multilayered perceptrons [START_REF] Suen | The combination of multiple classifiers by a neural network approach[END_REF]. Considering that the output of each PNN classifier may have uncertainty, in this paper, the PNN classifiers are combined within the framework of Dempster-Shafer theory (DST) [START_REF] Dempster | Upper and lower probabilities induced by multivalued mapping[END_REF], [START_REF] Shafer | A Mathematical Theory of Evidence[END_REF] due to its well capability of representing and combining uncertain information [START_REF] Jiao | Combining sources of evidence with reliability and importance for decision making[END_REF].

The rest of the paper is organized as follows. In Section II, the PW distance metric is defined and a parameter optimization procedure is designed to learn the involved feature weights. Based on the proposed PW distance metric, the corresponding PNN classifiers are combined within the framework of DST in Section III and then two experiments are given to evaluate the performance of the proposed method in Section IV. At last, Section V concludes the paper.

II. PAIRWISE-WEIGHTED DISTANCE METRIC

A. Definition

Definition 1 (Pairwise-weighted distance metric): Suppose x and y are two 𝐷-dimensional patterns whose class labels belong to Ω 𝑝,𝑞 = {𝜔 𝑝 , 𝜔 𝑞 }, the pairwise-weighted (PW) distance metric between x and y is defined as

𝑑 𝑃 𝑊 (x, y) = ⎷ 𝐷 ∑ 𝑗=1 𝜆 2 𝑝,𝑞,𝑗 (𝑥 𝑗 -𝑦 𝑗 ) 2 , ( 3
)
where 𝜆 𝑝,𝑞,𝑗 is a constant that weights the role of the 𝑗-th feature in the distance metric concerning the class pair Ω 𝑝,𝑞 .

This definition includes, as particular cases, the distance metrics revisited in the introduction. If 𝜆 𝑝,𝑞,𝑗 = 1 for all 𝑝 = 1, ⋅ ⋅ ⋅ , 𝑀, 𝑞 = 1, ⋅ ⋅ ⋅ , 𝑀, 𝑗 = 1, ⋅ ⋅ ⋅ , 𝐷, the above distance metric is just the L2 distance metric. Besides, the SW and CDW distance metrics correspond to the cases where the metric weights are irrelevant of the class labels or only dependent on the class label of the first pattern, respectively. Therefore, the above weighted distance metric provides a more general dissimilarity measure than the L2, SW or CDW because the weights depend on the labels of two considered classes.

Remark 1: Compared with the existing distance metrics reviewed above, the PW one provides more flexibility to design the feature weights so that the local specifics in feature space can be well characterized. Take the three-class classification problem studied in the introduction for example. In Fig. 1, using the PW distance metric, to classify Class B from Class A, 𝜆 𝐵,𝐴,𝑋 (the first two subscripts denote the class labels, the third subscript denotes the feature index) can takes much larger value than 𝜆 𝐵,𝐴,𝑌 , while one can assign smaller value for 𝜆 𝐵,𝐶,𝑋 than 𝜆 𝐵,𝐶,𝑌 to classify Class B from Class C.

B. Feature Weights Learning

In the previous subsection, the definition of the PW distance metric was given and the advantages of this proposed distance metric were also analyzed. The only open parameters in the PW distance metric are the feature weights related to the labels of the two considered classes. In this subsection we aim to learn the feature weights 𝜆 𝑝,𝑞,𝑗 (1 ≤ 𝑝 < 𝑞 ≤ 𝑀 , 𝑗 = 1, ⋅ ⋅ ⋅ , 𝐷) from the training set via optimizing some criteria. A simple way of defining the criteria for the desired metric is to keep the data pairs from the same class close to each other while separating those data pairs from different classes far with each other [START_REF] Xing | Distance metric learning with application to clustering with side-information[END_REF]. We divide the training set 𝒯 into 𝑀 subsets 𝒯 𝑘 , 𝑘 = 1, ⋅ ⋅ ⋅ , 𝑀, with each 𝒯 𝑘 containing all the 𝑛 𝑘 training data belonging to the same class 𝜔 𝑘 :

𝒯 𝑘 = {(x (𝑖) , 𝜔 𝑘 )|𝑖 ∈ 𝐼 𝑘 },
where 𝐼 𝑘 is the set of indices of the training data x (𝑖) belonging to the class 𝜔 𝑘 .

In the following, we consider learning the feature weights 𝜆 𝑝,𝑞,𝑗 (𝑗 = 1, ⋅ ⋅ ⋅ , 𝐷) from training subsets 𝒯 𝑝 and 𝒯 𝑞 . Let the set of data pairs from the same class be denoted by

𝒮 = { (x (𝑚) , x (𝑛) )|𝑚, 𝑛 ∈ 𝐼 𝑝 ; 𝑚 < 𝑛 } ∪ { (x (𝑚) , x (𝑛) )|𝑚, 𝑛 ∈ 𝐼 𝑞 ; 𝑚 < 𝑛 } ,
and the set of data pairs from different classes by

𝒟 = { (x (𝑚) , x (𝑛) )|𝑚 ∈ 𝐼 𝑝 ; 𝑛 ∈ 𝐼 𝑞 } .
Following the idea presented in [START_REF] Yang | Distance metric learning: a comprehensive survey[END_REF], a logistic regression model can be assumed when estimating the probability for any data pair (x (𝑚) , x (𝑛) ) to share the same class

Pr(+|(x (𝑚) , x (𝑛) )) = 1 1 + exp ( 𝑑 2 𝑃 𝑊 (x (𝑚) , x (𝑛) ) -𝜇 𝑝,𝑞
) , (4) and then the probability for any data pair (x (𝑚) , x (𝑛)

) to share different classes is Pr(-|(x (𝑚) , x (𝑛) )) = 1 - 1 1+exp(𝑑 2 𝑃 𝑊 (x (𝑚) ,x (𝑛) )-𝜇𝑝,𝑞) = 1 1+exp(-𝑑 2 𝑃 𝑊 (x (𝑚) ,x (𝑛) )+𝜇𝑝,𝑞) , ( 5
) where, "+" and "-" denote the data pair (x (𝑚) , x (𝑛) ) belonging to the same class and different classes, respectively. Parameter 𝜇 𝑝,𝑞 is the threshold. The data pair (x (𝑚) , x (𝑛) ) will be assigned more probability to be in the same class when their PW distance is much less then the threshold 𝜇 𝑝,𝑞 . In contrast, if their PW distance is much larger then the threshold 𝜇 𝑝,𝑞 , they will be given more probability to have different classes.

Then, the overall log-likelihood for both the data pairs in 𝒮 and 𝒟 can be written as

ℒ 𝑔 ( {𝜆 𝑝,𝑞,𝑗 } 𝐷 𝑗=1 , 𝜇 𝑝,𝑞 ) = log Pr(+|𝒮) + log Pr(-|𝒟) = - ∑ (x (𝑚) ,x (𝑛) )∈𝒮 log ( 1 + exp ( 𝑑 2 𝑃 𝑊 (x (𝑚) , x (𝑛) ) -𝜇 𝑝,𝑞 )) - ∑ (x (𝑚) ,x (𝑛) )∈𝒟 log ( 1 + exp ( -𝑑 2 𝑃 𝑊 (x (𝑚) , x (𝑛) ) + 𝜇 𝑝,𝑞 )) (6)
With the maximum likelihood estimation, the feature weights learning problem can be formulated as the following optimization problem max

The above constrained nonlinear optimization problem can be solved using Newton's method or some existing optimization software packages, such as the Matlab optimization toolbox. In a similar way, we can learn the feature weights 𝜆 𝑝,𝑞,𝑗 , 𝑗 = 1, ⋅ ⋅ ⋅ , 𝐷, for all other class pairs Ω 𝑝,𝑞 for 1 ≤ 𝑝 < 𝑞 ≤ 𝑀 , which are prepared for the following classifying process.

III. FUSION OF PNN CLASSIFIERS IN DST FRAMEWORK

Base on the proposed PW distance metric concerning class pair Ω 𝑝,𝑞 , a pairwise NN (PNN) classifier can be designed to separate the two classes Ω 𝑝 and Ω 𝑞 . For an 𝑀 -class classification problem, 𝑀 (𝑀 -1)/2 PNN classifers 𝒞 𝑝,𝑞 (1 ≤ 𝑝 < 𝑞 ≤ 𝑀 ) can be designed and the final classification result can be obtained by combining the output of these PNN classifers. A popular method of combining the output of pairwise classifiers is the voting strategy, where each classifier gives a vote for the predicted class and the class with the largest number of votes is predicted. However, a classifier 𝒞 𝑝,𝑞 is trained to distinguish only between classes Ω 𝑝 and Ω 𝑞 , thus its vote for an query pattern from a different class should be taken with care. In this section, we aim to overcome this difficulty by modeling the uncertainty related to the output of PNN classifiers in the framework Dempster-Shafer theory (DST).

A. Basics of DST

In DST, a problem domain is represented by a finite set Ω = {𝜔 1 , 𝜔 2 , ⋅ ⋅ ⋅ , 𝜔 𝑛 } of mutually exclusive and exhaustive hypotheses called the frame of discernment. A mass function or basic belief assignment (BBA) expressing the belief committed to the elements of 2 Ω by a given source of evidence is a mapping function m(⋅):

2 Ω → [0, 1], such that m(∅) = 0 and ∑ 𝐴∈2 Ω m(𝐴) = 1. (8)
Elements 𝐴 ∈ 2 Ω having m(𝐴) > 0 are called focal elements of the BBA m(⋅). The BBA m(𝐴) measures the degree of belief exactly assigned to a proposition 𝐴 and represents how strongly the proposition is supported by evidence. The belief assigned to Ω, or m(Ω), is referred to as the degree of global ignorance.

Shafer [START_REF] Shafer | A Mathematical Theory of Evidence[END_REF] For decision-making, the maximum plausibility (as defined in Eq.( 9)) rule is usually utilized to make the final decision.

Suppose m = m 1 ⊕ m 2 , then Pl({𝜔 𝑖 }) ∝ Pl 1 ({𝜔 𝑖 })Pl 2 ({𝜔 𝑖 }), ∀𝜔 𝑖 ∈ Ω. (11)
That is, when combining several mass functions, we do not need to compute the complete mass function using Dempster's rule of combination. Instead, we can compute the combined plausibility using Eq.( 11) to make the decision equivalently.

In order to manipulate the belief functions more effectively, some probabilistic operations (conditioning, deconditioning, etc) are introduced to DST framework [START_REF] Smets | Belief functions: the disjunctive rule of combination and the generalized Bayesian theorem[END_REF]. Conditional beliefs represent knowledge which is valid provided that a hypothesis is satisfied. Let m Ω (⋅) be a BBA on Ω, 𝑆 ⊆ Ω an hypothesis and m Ω 𝑆 (⋅) the categorical BBA such that m Ω 𝑆 (𝑆) = 1. Then the conditional belief function m

Ω [𝑆](⋅) is m Ω [𝑆] = m Ω ⊕ m Ω 𝑆 . ( 12
)
The above operation is referred to as Dempster's rule of conditioning. In contrary, if we want to recover m Ω (⋅) from the conditional belief function m Ω [𝑆](⋅), the following deconditioning operation can be used

m Ω (𝐴 ∪ 𝑆) = m Ω [𝑆](𝐴), for 𝐴 ∈ 2 Ω ∖ ∅, (13)
where, 𝑆 denotes the complement of set 𝑆 with respect to set Ω.

B. Combination of PNN Classifiers

Our aim is to use the Dempster-Shafer theory to model the uncertainty inherent in the pairwise classification. Now, we have a set of PNN classifiers 𝒞 𝑝,𝑞 (1 ≤ 𝑝 < 𝑞 ≤ 𝑀 ) and we first look for the representation of their output in terms of the DST.

For the output of each PNN classifier 𝒞 𝑝,𝑞 , there are two kinds of uncertainty. The first one is related to the fact that the real class label of query pattern y may actually not belong to class pair Ω 𝑝,𝑞 (called outer-pair uncertainty). The second one is that even the real class label of query pattern y belongs to class pair Ω 𝑝,𝑞 , affected by the noises of the training patterns, the result of the classifier is not always accurate (called innerpair uncertainty). So, the frame of discernment should be set as Ω = {𝜔 1 , ⋅ ⋅ ⋅ , 𝜔 𝑀 }, which can characterize both kinds of uncertainty.

For the PNN classifier 𝒞 𝑝,𝑞 , suppose the class label of the nearest neighbor is 𝜔 𝑝 ∈ Ω 𝑝,𝑞 . It can be seen as a piece of evidence that support the query pattern y belongs to 𝜔 𝑝 . However, considering the outer-pair uncertainty, this piece of evidence is conditioned on the hypothesis 𝜔 𝑝 ∈ Ω 𝑝,𝑞 :

m Ω [Ω 𝑝,𝑞 ]({𝜔 𝑝 }) = 1. (14)
Further, due to the inner-pair uncertainty, this piece of evidence does not by itself provide 100% reliability. In DST formalism, this can be expressed by saying that only some part of the belief is committed to

𝜔 𝑝 { m Ω [Ω 𝑝,𝑞 ]({𝜔 𝑝 }) = 𝛼 𝑝,𝑞 m Ω [Ω 𝑝,𝑞 ]({𝜔 𝑞 }) = 1 -𝛼 𝑝,𝑞 , (15)
where, 𝛼 𝑝,𝑞 ∈ [0, 1] represents the reliability that the PNN classifier 𝒞 𝑝,𝑞 provides the right classification result between classes 𝜔 𝑝 and 𝜔 𝑞 .

In determining the reliability factor 𝛼 𝑝,𝑞 , the probability model in Eq.( 4) for data pair from the same class can be used

𝛼 𝑝,𝑞 = 1 1 + exp (𝑑 2 𝑃 𝑊 -𝜇 𝑝,𝑞 ) , ( 16
)
where, 𝑑 𝑃 𝑊 denotes the PW distance between the query pattern y and its nearest neighbor and 𝜇 𝑝,𝑞 is the threshold learned in Eq. [START_REF] Paredes | A class-dependent weighted dissimilarity measure for nearest neighbor classication problems[END_REF].

In order to combine the output of different PNN classifiers in a uniform framework, the conditional BBA constructed as Eq.( 15) should be deconditioned using Eq.( 13) as

{ m Ω 𝑝,𝑞 ({𝜔 𝑞 }) = m Ω [Ω 𝑝,𝑞 ]({𝜔 𝑝 }) = 𝛼 𝑝,𝑞 m Ω 𝑝,𝑞 ({𝜔 𝑝 }) = m Ω [Ω 𝑝,𝑞 ]({𝜔 𝑞 }) = 1 -𝛼 𝑝,𝑞 . ( 17
)
Because the mass function and plausibility function are in one-to-one correspondence, so we can compute the plausibility function Pl 𝑝,𝑞 (⋅) from the above constructed and deconditioned BBA m Ω 𝑝,𝑞 (⋅) using Eq.( 9) as

Pl 𝑝,𝑞 ({𝜔 𝑖 }) = { 𝛼 𝑝,𝑞 , if 𝑖 = 𝑝 1 -𝛼 𝑝,𝑞 , if 𝑖 = 𝑞 1,
others .

(

) 18
In order to decrease the computation complexity, instead of combining the 𝑀 (𝑀 -1)/2 BBAs m Ω 𝑝,𝑞 (⋅) (1 ≤ 𝑝 < 𝑞 ≤ 𝑀 ) using Dempsters rule of combination, we can compute the combined plausibility function Pl(⋅) directly using Eq.( 11) to make the decision as follows

Pl({𝜔 𝑖 }) ∝ Pl ′ ({𝜔 𝑖 }) = ∏ 1≤𝑝<𝑞≤𝑀 Pl 𝑝,𝑞 ({𝜔 𝑖 }), ∀𝜔 𝑖 ∈ Ω. (19)
Note that the combined plausibility function Pl(⋅) is proportional to Pl ′ (⋅), so the maximum plausibility rule can be used for Pl ′ (⋅) equivalently to make a hard decision. The class label of query pattern y is assigned to the class with maximum plausibility criteria.

Remark 2:

As can be seen from above, the combination process of PNN classifiers is quite time-saving. So, when classifying a query pattern, the time is mainly consumed in the classification process of multiple PNN classifiers. Even though the number of PNN classifiers is of 𝑀 2 order (with 𝑀 (𝑀 -1)/2 classifiers), each classifier only uses the training samples from the corresponding classes (about 2𝑁/𝑀 samples averagely). Hence the total number of the computed samples is about 𝑁 (𝑀 -1), which is just 𝑀 -1 times larger than the original NN classifier. To most classification problems, such as the benchmark data sets studied in next section, the number of considered classes is not very large, so the computation cost of the proposed method is not a big problem.

IV. NUMERICAL EXPERIMENTS

The capabilities of the proposed NN classification method based on the PW distance metric and Dempster-Shafer theory (denoted as PNN-DST) was empirically assessed through two different types of experiments. In the first experiment, a synthetic data set was used to show the behavior of the proposed method in controlled setting. In the second one, several standard benchmark data sets from the well-known UCI Repository of Machine Learning Databases [START_REF] Merz | UCI Repository of Machine Learning Databases[END_REF] were considered, with the aim to show that the proposed technique is uniformly adequate for a variety of tasks involving different data conditions: large/small data sets, high/low dimensionality, etc.

A. Synthetic Data

A three-class classification problem was used to compare our method with the original NN classifier based on L2 distance metric (L2-NN) [START_REF] Fix | Discriminatory analysis, nonparametric discrimination: consistency properties[END_REF] and the NN classifier based on CDW distance metric (CDW-NN) [START_REF] Paredes | A class-dependent weighted dissimilarity measure for nearest neighbor classication problems[END_REF], [START_REF]Learning weighted metrics to minimize nearest-neighbor classification error[END_REF]. The following class-conditional normal distributions were assumed. training sets randomly drawn from the above distributions. A fixed test set of 300 query patterns, independently drawn from the same distributions, is used for error statistics. As can be seen from the result, the CDW-NN classifier is just slightly better than the original L2-NN classifier for small training sets. This is mainly because the CDW distance metric only characterizes the local specifics of the features with respect to the prototypes, while, in this simulation scenario, for each class, the statistical variances of the prototypes in different features (i.e. X axis, Y axis) are almost the same. Under these circumstances, the CDW distance metric asymptotically approximates to the L2 distance metric as the amount of training data increases. The proposed PNN-DST classifier produces the lowest error rate as the PW distance metric provides more local specifics of the features for each PNN sub-classifier and further in the combination process the output uncertainty of those PNN sub-classifiers is well addressed. Besides, the performance improvement is more significant for small training sets, in which case the ideal asymptotical behavior of NN classifier degrades dramatically.

To better illustrate the superiority of the proposed PNN-DST classifier, the classification results of different methods for one Monte Carlo run with 30 training data and 30 test data are given in Fig. 3. As can be seen in Fig. 3

B. Benchmark Data Sets

In this second experiment, ten well-known benchmark data sets from UCI repository were used to evaluate the performance of the PNN-DST classifier. The main characteristics of the data sets are summarized in Table I. In order to evaluate the effectiveness of the combination process with DST, apart from the L2-NN and CDW-NN methods, we also consider the method of fusing the PNN classifiers with voting strategy (denoted as PNN-VOTE).

The classification results of the ten benchmark data sets are shown in Table II. The significance of the differences between results is evaluated using a Mc Nemar test [START_REF] Dietterich | Approximate statistical tests for comparing supervised classification learning algorithms[END_REF] at level 5%. For each dataset, the best classification accuracy is underlined, and those that are significantly improved than the baseline L2-NN method are printed in bold. As can be seen from these results, the PNN-DST method, presented in this paper, outperforms the L2-NN or CDW-NN for most of the data sets. Additionally, for Ecoli, Glass, Vehicle, Waveform, Wine and Yeast data sets, the improvements are statistically significant than the baseline L2-NN method, because the local distance metric plays more crucial role in determining the NN based classification performance for these scarce-prototype and high-dimensionality cases. On the other hand, the PNN-DST method always performs better than PNN-VOTE, especially for those data sets with small number of classes, because the voting strategy will take great adventure when the total number of votes (𝑀 (𝑀 -1)/2, with 𝑀 be the number of classes) is small.

V. CONCLUDING REMARKS

In order to improve the performance of the NN based classifier in small data set situations, a new distance metric called PW distance metric, has been proposed in this paper. Compared with the existing distance metrics, the PW one provides more flexibility to design the feature weights so that the local specifics in feature space can be well characterized. A parameter optimization procedure is designed to learn the feature weights from the training data set. Based on the PW distance metric, a PNN-DST classifier is developed, which combines the output of PNN classifiers using Dempster-Shafer theory. From the results reported in the last section, we can conclude that the proposed method achieved a uniformly good performance when applied to a variety of classification tasks, including those with high dimensionality and sparse prototypes.

ACKNOWLEDGMENT

  also defines the belief function and plausibility function of 𝐴 ∈ 2 Ω as follows Bel(𝐴) = ∑ 𝐵⊆𝐴 m(𝐵) and Pl(𝐴) = ∑ 𝐵∩𝐴∕ =∅ m(𝐵). (9) Bel(𝐴) represents the exact support to 𝐴 and its subsets, and Pl(𝐴) represents all the possible support to 𝐴 and its subsets. The belief functions m(⋅), Bel(⋅) and Pl(⋅) are in one-to-one correspondence.Several distinct bodies of evidence characterized by different BBAs can be combined using Dempster's rule of combination ⊕ . Mathematically, the Dempster's rule of combination of two BBAs m 1 (⋅) and m 2 (⋅) defined on the same frame of discernment Ω Ω ;𝐵∩𝐶=∅ m1(𝐵)m2(𝐶) , for 𝐴 ∈ 2 Ω ∖ ∅.[START_REF] Jahromi | A method of learning weighted similarity function to improve the performance of nearest neighbor[END_REF] As described in[START_REF] Shafer | A Mathematical Theory of Evidence[END_REF], Dempster's rule of combination is both commutative and associative.

  Class A: 𝜇 𝐴 = (6, 6), Σ 𝐴 = 𝑑𝑖𝑎𝑔(2, 2); Class B: 𝜇 𝐵 = (14, 6), Σ 𝐵 = 𝑑𝑖𝑎𝑔(2, 2); Class C: 𝜇 𝐶 = (14, 14), Σ 𝐶 = 𝑑𝑖𝑎𝑔(2, 2).

Fig. 2 Fig. 2 .

 22 Fig.2 shows classification error rates with different training set sizes, ranging from 5 to 100 prototypes per class. For each size, each classification algorithm runs 100 times with different

  (a), the test data y 1 , y 2 are quite close to the boundary between Class B and Class C, and in this scarce-prototype condition, it is quite difficult to make the right classification. The L2-NN and CDW-NN just misclassify these two data points as shown in Fig.3(b) and (c), because based on L2 and CDW distance metrics their nearest neighbors are the training data x 1 labeled Class C and x 2 labeled Class B, respectively. However, as shown in Fig.3(d), the test data y 2 is correctly classified based on the PNN-DST classifier. The PNN-DST classifier classifies the test data y 2 by combining the results of three PNN sub-classifiers: Pl 𝐴,𝐵 ({𝐴}) = 0.27, Pl 𝐴,𝐵 ({𝐵}) = 0.73, Pl 𝐴,𝐵 ({𝐶}) = 1; Pl 𝐴,𝐶 ({𝐴}) = 0.29, Pl 𝐴,𝐶 ({𝐵}) = 1, Pl 𝐴,𝐶 ({𝐶}) = 0.71; Pl 𝐵,𝐶 ({𝐴}) = 1, Pl 𝐵,𝐶 ({𝐵}) = 0.31, Pl 𝐵,𝐶 ({𝐶}) = 0.69. Then after the fusion of multiple PNN sub-classifiers within DST framework, we can get the combined result: Pl ′ ({𝐴}) = 0.08, Pl ′ ({𝐵}) = 0.23, Pl ′ ({𝐶}) = 0.49. Finally, based on the maximum plausibility rule we get the final classification result Class C with hard partition.

  This work is partially supported by China Natural Science Foundation (No. 61135001, 61374159) and the Doctorate Foundation of Northwestern Polytechnical University (No. CX201319).

2 Fig. 3 .

 23 Fig. 3. Classification results for different methods with 30 training data and 30 test data. (a) Training data and test data. (b) Classification results by L2-NN. (c) Classification results by CDW-NN. (d) Classification results by PNN-DST. (The blue makers represent the training data, and the red makers represent the test data, with circle for class A and square for class B and triangle for class C, respectively. The filled makers represent the data mistakenly classified.)

TABLE I .

 I STATISTICS OF THE BENCHMARK DATA SETS USED IN THE EXPERIMENT.

	Data set	# of instances	# of features	# of classes	# of training patterns	# of test patterns
	Balance	625	4	3	500	125
	Ecoli	336	7	8	200	136
	Glass	214	9	6	139	75
	Letter	20,000	16	26	15,000	5,000
	Satimage	6,435	36	6	4,435	2,000
	Segment	2,310	19	7	1,400	910
	Vehicle	846	18	4	646	200
	Waveform	5,000	21	3	3,500	1,500
	Wine	178	13	3	75	103
	Yeast	1,484	8	10	1000	484

TABLE II .

 II CLASSIFICATION ACCURACY (IN %) OF OUR METHOD IN COMPARISON WITH OTHER NN-BASED METHODS.

	Data set	L2-NN	CDW-NN	PNN-VOTE	PNN-DST
	Balance	74.40	79.20 a	75.20	77.60
	Ecoli	80.88	83.82	86.76 b	88.24
	Glass	69.33	72.00	72.00	74.67
	Letter	94.26	94.72	95.40	95.80
	Satimage	89.30	88.55	91.55	92.55
	Segment	95.05	96.15	96.70	96.92
	Vehicle	69.00	69.50	70.50	74.50
	Waveform	80.53	83.87	84.53	85.87
	Wine	76.70	89.32	87.38	90.29
	Yeast	51.65	55.79	57.64	57.85

a The results underlined correspond to the best accuracy.

b

The results typeset in boldface are significantly better than the baseline L2-NN method at level 5%.


K Autchariyapanitkul

S Chanaim

S Sriboonchitta

T Denoeux

Predicting stock returns in the capital asset pricing model using quantile regression and belief functions

Keywords: Asymmetric Laplace distribution, Capital Asset Pricing, Dempster-Shafer Theory, Financial data, Likelihood-based belief functions, Quantile regression

We consider an inference method for prediction based on belief functions in quantile regression with an asymmetric Laplace distribution. We apply this method to the capital asset pricing model to estimate the beta coefficient and measure volatility under various market conditions at given quantiles. Likelihoodbased belief functions are constructed from historical data of the securities in the S&P500 market. The results give us evidence on the systematic risk, in the form of a consonant belief function specified from the asymmetric Laplace distribution likelihood function given recorded data. Finally, we use the method to forecast the return of an individual stock.

Introduction

The Capital Asset Pricing Model (CAPM) is one of the most useful models in investment. In this model, asset returns are usually assumed to be jointly normally distributed random variables. However, this is not always the case. The CAPM assumes that the variance of returns adequately measures risk. This may be true if returns are normally distributed. In this paper, we propose to use quantile regression with an asymmetric Laplace distribution (ALD), coupled with an inference method based on belief functions, to estimate the parameters of the model and predict stock returns.

Quantile regression can characterize the entire conditional distribution of the outcome variable and is more robust to outliers and misspecification of the error distribution. It can also handle heteroscedasticity, as shown by Koenger [START_REF] Koenker | Robust tests for heteroscedasticity based on regression quantiles[END_REF]. For the application of quantile regression to the CAPM, the reader is referred to, e.g., [START_REF] Barnes | A Quantile Regression Analysis of the Cross Section of Stock Market Returns[END_REF]. It was found that the market price of beta risk is significant in both tails of the conditional distribution of returns.

We present a likelihood-based approach to the estimation of regression quantiles based on the asymmetric Laplace distribution. In [START_REF] Sánchez | Likelihood Based Inference for Quantile Regression Using the Asymmetric Laplace Distribution[END_REF], this distribution is used to model the distribution of currency exchange rates and is shown to capture the peakedness, leptokurticity (Fat tails) and skewness inherent in such data. Similarly, it is shown in [START_REF] Linden | A Model for Stock Return Distribution[END_REF] that the Laplace distribution has a geometric stability to represent the weekly and monthly distributions of stock returns and also models the high peak, fat tails and skewness of the returns.

Here, we use the Dempster-Shafer theory of belief functions introduced by Dempster [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF] and Shafer [START_REF] Shafer | A mathematical theory of evidence[END_REF]. In this approach, a piece of evidence is modeled by a belief function, which can be viewed as the distribution of a random set. This method is applied to estimation using the likelihood-based approach introduced in [START_REF] Shafer | A mathematical theory of evidence[END_REF] and recently justified in [START_REF] Denoeux | Likelihood-based belief function: justification and some extensions to lowquality data[END_REF], and to prediction using the method introduced in [START_REF] Kanjanatarakul | Forecasting using belief functions: an application to marketing econometrics[END_REF]. The main contribution of this paper is thus to propose an alternative method for drawing inferences about conditional quantiles via a likelihood-based belief function approach.

The remainder of the paper is organized as follows. Section 2 provides the background on quantile regression with asymmetric Laplace distribution and Section 3 introduces the prediction machinery using belief functions. Section 4 discusses the empirical solutions to the forecasting problem. The last section summarizes the paper.

Quantile regression with an Asymmetric Laplace distribution

Let Y be a response variable and X a vector of explanatory variables. In linear quantile regression, the conditional α-quantile q α (Y |X) of Y given X is assumed to be linearly related to X through the equation q α (Y |X) = X β α , where β α is a vector of unknown parameters and X is the transpose of X. Denoting the error by ε α , we can write the quantile regression model as

Y = X β α + ε α . (1)
We have

q α (Y |X) = q α [(X β α + ε α )|X] = q α (X β α + ε α |X) = X β α + q α (ε α |X), (2)
since given X, X β α is a constant. Thus, q α (ε α |X) = 0, which is the counterpart of the standard condition E(ε|X) = 0 in the mean linear regression model. If ε α is independent of X, the α-quantile of the noise ε α is zero, that is, 0

-∞ dF ε α (u) = α. For q α (Y |X) = X β α , we see that β α minimize E[ρ α (Y -X β )] over β
, where ρ α (•) is the so-called check (or loss) function defined by

ρ α (u) = u(α -1 (u<0) ), (3)
with 1 (u<0) denoting the usual indicator function. Thus, given i.i.d

(X i ,Y i ), a plausible estimator of β α is β α = arg min 1 n { n ∑ i=1 ρ α (Y i -X i β )}. ( 4
)
This estimator is called the Least Absolute Deviation (LAD) estimator. Suppose that the error term ε α has an ALD with mean 0 and standard deviation σ α :

f σ α (ε α ) = α(1 -α) σ α exp -ρ α ε α σ α . (5)
Then, minimizing the absolute deviation is equivalent to maximizing the likelihood and the LAD estimator of β α is a maximum likelihood estimator (MLE). The likelihood function for β α , σ α after observing the data

D = (X 1 ,Y 1 ), . . . , (X n ,Y n ) is L D (β α , σ α ) = α n (1 -α) n σ n α exp - n ∑ i=1 ρ α Y i -X i β α σ . ( 6
)
3 Statistical inference and prediction using Belief functions

Likelihood-based belief functions

Suppose we observe a realization x of the random vector X with probability density function (pdf) p θ (x), where θ ∈ Θ is an unknown parameter. In this paper, we use the method proposed by Shafer [START_REF] Shafer | A mathematical theory of evidence[END_REF], which can be derived from the Likelihood Principle (LP) and the Least Commitment Principle (LCP) [START_REF] Denoeux | Likelihood-based belief function: justification and some extensions to lowquality data[END_REF]. According to the LP, all the information about Θ is represented by the likelihood function defined by L x (θ ) = p θ (x) for all θ ∈ Θ . In statistics, the likelihood ratio has the meaning of a "relative plausibility", which can be written as:

pl x (θ 1 ) pl x (θ 2 ) = L x (θ 1 ) L x (θ 2 ) , (7)
for all (θ 1 , θ 2 ) ∈ Θ 2 or, equivalently, pl x (θ ) = cL x (θ ), for all θ ∈ Θ and some positive constant c. The LCP then implies that the highest possible value should be given to constant c [START_REF] Denoeux | Likelihood-based belief function: justification and some extensions to lowquality data[END_REF], which leads us to equating the contour function pl x with the relative likelihood:

pl x (θ ) = L x (θ ) sup θ ∈Θ L x (θ ) . ( 8
)
The information about θ is represented by the consonant belief function Bel Θ x with contour function pl x , i.e., with corresponding plausibility function Pl Θ

x (A) = sup θ ∈A pl x (θ ), for all A ⊆ Θ . The focal sets of Bel Θ

x are the levels sets of pl x defined as follows:

Γ x (ω) = {θ ∈ Θ |pl x (θ ) ≥ ω}, (9)
for θ ∈ [0, 1]. These sets are also called plausibility regions. The consonant belief function Bel Θ x is equivalent to the random set induced by the Lebesgue measure λ on [0,1] and the multi-valued mapping Γ x from [0, 1] → 2 Θ (see, [START_REF] Kanjanatarakul | Forecasting using belief functions: an application to marketing econometrics[END_REF]). We remark that the MLE of θ is the value of θ with highest plausibility.

Prediction using belief functions

Let X be a random variable with parametric density function f θ (x) for θ ∈ Θ and assume that we have observed X = x. Given the belief function Bel Θ

x about θ , we can predict the future value of a random variable Y whose pdf g θ (y) also depends on θ . In the approach introduced in [START_REF] Kanjanatarakul | Forecasting using belief functions: an application to marketing econometrics[END_REF], Y is written as a function of the parameter θ and an unobserved auxiliary variable u ∈ U with known probability measure µ not depending on θ :

Y = ϕ(θ , u). (10)
Using Equations ( 9) and (10), we can compose the multi-valued mapping Γ x from [0, 1] → 2 Θ with ϕ to get a new multi-valued mapping Γ x from [0, 1] × U to 2 Y defined as

Γ x : [0, 1] × U → 2 Y (ω, u) → ϕ(Γ x (ω), u).
We can then define the predictive belief (Bel Y x ) and plausibility (Pl Y x ) functions on Y as

Bel Y x (A) = (λ ⊗ µ)({(ω, u) ∈ [0, 1]|ϕ(Γ x (ω, u) ⊆ A}). ( 11a
)
Pl Y x (A) = (λ ⊗ µ)({(ω, u) ∈ [0, 1]|ϕ(Γ x (ω, u) ∩ A = / 0}), (11b)
for all A ⊆ Y.

4 Application to stock market prediction

Model

The CAPM measures the sensitivity of the expected excess return on security to expected market risk premium. The equation of CAPM is a linear function of the security market line:

E(R A ) -R F = β 0 + β 1 E(R M -R F ), (12)
where E(R A ) is the expected return of the asset, R M is the expected market portfolio return, R F is the risk free rate, β 0 is the intercept and β 1 is the equity beta, representing market risk. Suppose we have observed the historical returns of stock R A = (r a1 , • • • , r an ) and returns from market R M = (r m1 , • • • , r mn ). The errors will be assumed to be iid with density function [START_REF] Kanjanatarakul | Forecasting using belief functions: an application to marketing econometrics[END_REF]. The likelihood function is given by (6).

Experimental results

The data contain the weekly log returns of the integrated oil and gas company, Chesapeake Energy (CHK), during 2010-2013. The ML estimates of the parameters are shown in Table 1 for different values of α. We used a nonlinear optimization algorithm to maximize the likelihood (6) with respect to θ = (β (0,α) , β (1,α) , σ α ). The plausibility function on θ is then defined by the relative likelihood [START_REF] Shafer | A mathematical theory of evidence[END_REF] and the marginal contour function on a specific parameter is obtained by take the supremum with respect to the others parameters, e.g.,

pl R A (β (0,α) ) = sup β (1,α) ,σ pl R A (β (0,α) ). ( 13
)
Figure 1 displays two-dimensional marginal contour functions and Figure 2 shows the marginal contour functions for parameters β 0 , β 1 and σ . The 0.15 threshold corresponds to an approximate 95% confidence interval and gives us an interval of plausible values of each of the three parameters. To predict the expected return of the asset r a,n+1 for a new market portfolio return r m,n+1 , we compute the minimum and maximum of r a,n+1 given r m,n+1 at fixed α by

r a,n+1 = β 0,α + β 1,α r m,n+1 + σ α F -1 ε α (u), (14)
under the constraint pl R A (β i,α , σ α ) > ω, where F -1 ε α is the inverse cumulative distribution function (cdf) of the asymmetric Laplace distribution ALD(α, 0, 1) and u, ω are independent random variables with the same uniform distribution U([0, 1]). Given (14), we randomize independently N pairs of the random number

(ω i , u i ), i = 1, 2, • • • , N re- sulting in N intervals [r L a (ω i , u i ), r U a (ω i , u i )].
For any A ⊆ R, the stock returns Bel r ai (A) and Pl r ai (A) can be estimated by equation (11). The estimated lower and upper expectations of r a,n+1 are then: Bel R A . Figure 4 shows the upper and lower predictive quantiles of the stock returns (see [START_REF] Kanjanatarakul | Forecasting using belief functions: an application to marketing econometrics[END_REF]), defined by the inequalities pl(R A q L α ) = α and pl(R A q U α ) = α . As shown in [START_REF] Kanjanatarakul | Forecasting using belief functions: an application to marketing econometrics[END_REF], the following inequalities hold:

R L A = N ∑ i=1 r L a (ω i , u i ) N (15a) R U A = N ∑ i=1 r U a (ω i , u i ) N . (15b)
Bel(q L α R A q U 1-α ) 1 -2α . (16)
Lower-upper expectations and quantiles are other representations of prediction uncertainty, taking into account both parameter estimation uncertainty and randomness. For these data, the very narrow gap between the lower and upper cdfs shows that estimation uncertainty is small as compared to random uncertainty. In practice, these results can be used to increase the performance of the investment portfolio.

Conclusions

In this paper, we studied the method of quantile CAPM with ALD for stocks in S&P500 in the belief function framework. We used the Dempster-Shafer theory of belief functions to model the uncertainty referring to the statistical prediction based on historical data and a financial model. This method consists of two steps. First, a consonant belief function representing the uncertainty on the parameter vector θ is defined from the normalized likelihood function given the past data. The return of individual stock R A is then expressed as ϕ(θ , u), where u is a stochastic variable with known distribution and the beliefs on θ and u are transferred through ϕ, resulting in a belief function on R A . This approach has been adapted to the prediction of stock returns. A possible extension of this work is to consider uncertainty on the independent variable r m , which can also be expressed as a belief function and combined with other uncertainties to compute a belief function on R A . This issue will be addressed in future work.

Fig. 1 .

 1 Fig. 1. Marginal contour functions pl x (θ 1 , θ 2 ) in two-dimensional parameter subspaces.

Figures 3 5 Fig. 2 .

 352 Figures 3 displays the lower and upper cdfs Bel R A ((-∞, R A ]) and Pl R A ((-∞, R A ]) at given r m = 0.05. This function give us the summary of the predictive belief function

Fig. 3 .

 3 Fig. 3. Upper and lower cumulative distribution function given r m = 0.05.

  α'=.05

Fig. 4 .

 4 Fig. 4. Lower and upper quantiles of stock returns at α = 0.5.

Table 1 .

 1 Parameter estimation results. Standard errors are given in parentheses.

	Stock Name	Parameters	α = 0.40	α = 0.50	α = 0.60
	CHK	β 0	-0.011 (0.002) -0.004 (0.003) 0.005 (0.000)
		β 1	1.379 (0.163) 1.442 (0.005) 1.304 (0.017)
		σ	0.002 (0.001) 0.016 (0.001) 0.016 (0.010)


Orakanya Kanjanatarakul

Nachatchapong Kaewsompong

Songsak Sriboonchitta

Thierry Denoeux

Estimation and prediction using belief functions: Application to stochastic frontier analysis

We outline an approach to statistical inference based on belief functions. For estimation, a consonant belief functions is constructed from the likelihood function. For prediction, the method is based on an equation linking the unobserved random quantity to be predicted, to the parameter and some underlying auxiliary variable with known distribution. The approach allows us to compute a predictive belief function that reflects both estimation and random uncertainties. The method is invariant to one-to-one transformations of the parameter and compatible with Bayesian inference, in the sense that it yields the same results when provided with the same information. It does not, however, require the user to provide prior probability distributions. The method is applied to stochastic frontier analysis with crosssectional data. We demonstrate how predictive belief functions on inefficiencies can be constructed for this problem and used to assess the plausibility of various assertions.

Introduction

Many problems in econometrics can be formalized using a parametric model

(Y, Z)|x ∼ f θ ,x (y, z), (1)
where Y and Z are, respectively, observed and unobserved random vectors, x is an observed vector of covariates and f θ ,x is the conditional probability mass or density function of (Y, Z) given X = x, assumed to be known up to a parameter vector θ ∈ Θ . For instance, in the standard linear regression model, Y = (Y 1 , . . . ,Y n ) is a vector of n independent observations of the response variable, with Y i ∼ N (x i β , σ 2 ), Z = Y n+1 is an independent random value of the response variable distributed as N (x n+1 β , σ 2 ), x = (x 1 , . . . , x n+1 ) and θ = (β , σ 2 ). Having observed a realization y of Y (and the covariates x), we often wish to determine the unknown quantities in the model, i.e., the parameter θ (assumed to be fixed) and the (yet) unobserved realization z of Z. The former problem is referred to as estimation and the latter as prediction (or forecasting). These two problems have been addressed in different ways within several theoretical frameworks. The three main theories are frequentist, Bayesian and likelihoodbased inference. In the following, we briefly review these three approaches to introduce the motivation for the new method advocated in this paper.

Frequentist methods provide pre-experimental measures of the accuracy of statistical evidence. A procedure (for computing, e.g., a confidence or prediction interval) is decided before observing the data and its long-run behavior is determined by averaging over the whole sample space, assuming it is repeatedly applied to an infinite number of samples drawn from the same population. It has long been recognized that such an approach, although widely used, does not provide a reliable measure of the strength of evidence provided by specific data. The following simple example, taken from [START_REF] Berger | The likelihood principle: a review, generalizations, and statistical implications[END_REF], illustrates this fact. Suppose X 1 and X 2 are iid with probability mass function

P θ (X i = θ -1) = P θ (X i = θ + 1) = 1 2 , i = 1, 2, (2)
where θ ∈ R is an unknown parameter. Consider the following confidence set for θ ,

C(X 1 , X 2 ) = 1 2 (X 1 + X 2 ) if X 1 = X 2 X 1 -1 otherwise. ( 3
)
It is a minimum length confidence interval at level 75%. Now, let (x 1 , x 2 ) be a given realization of the random sample (X 1 , X 2 ). If x 1 = x 2 , we know for sure that θ = (x 1 + x 2 )/2 and it would be absurd to take 75% as a measure of the strength of the statistical evidence. If x 1 = x 2 , we know for sure that θ is either x 1 -1 or x 1 + 1, but we have no reason to favor any of these two hypotheses in particular. Again, it would make no sense to claim that the evidence support the hypothesis θ = x 1 -1 with 75% confidence. Although frequentist procedures do provide usable results in many cases, the above example shows that they are based on a questionable logic if they are used to assess the reliable of given statistical evidence, as they usually are. Moreover, on a more practical side, confidence and prediction intervals are often based on asymptotic assumptions and their true coverage probability, assuming it is of interest, may be quite different from the nominal one for small sample sizes. The other main approach to statistical inference is the Bayesian approach, which, in contrast to the previous approach, implements some form of post-experimental reasoning. Here, all quantities, including parameters, are treated as random variables, and the inference aims at determining the probability distribution of unknown quantities, given observed ones. With the notations introduced above, the estimation and prediction problems are to determine the posterior distributions of, respectively, θ and Z, given x and y. Of course, this is only possible if one provides a prior probability distribution π(θ ) on θ , which is the main issue with this approach. There has been a long-standing debate among statisticians about the possibility to determine such a prior when the experimenter does not know anything about the parameter before observing the data. For lack of space, we cannot reproduce all the arguments of this debate here. Our personal view is that no probability distribution is truly non-informative, which weakens the conclusions of Bayesian inference in situations where no well-justified prior can be provided.

The last classical approach to inference is grounded in the likelihood principle (LP), which states that all the information provided by the observations about the parameter is contained in the likelihood function. A complete exposition of the likelihood-based approach to statistical inference can be found in the monographs [START_REF] Edwards | Likelihood[END_REF] and [START_REF] Berger | The likelihood principle: a review, generalizations, and statistical implications[END_REF] (see also the seminal paper of Barnard et al. [START_REF] Barnard | Likelihood inference and time series[END_REF]). Birnbaum [START_REF] Birnbaum | On the foundations of statistical inference[END_REF] showed that the LP can be derived from the two generally accepted principles of sufficiency and conditionality. Frequentist inference does not comply with the LP, as confidence intervals and significance tests depend not only on the likelihood function, but also on the sample space. Bayesian statisticians accept the LP, but claim that the likelihood function does not make sense in itself and needs to be multiplied by a prior probability distribution to form the posterior distribution of the parameter given the data. The reader is referred to Refs. [START_REF] Edwards | Likelihood[END_REF] and [START_REF] Berger | The likelihood principle: a review, generalizations, and statistical implications[END_REF] for thorough discussions on this topic. Most of the literature on likelihood-based inference deals with estimation. Several authors have attempted to address the prediction problem using the notion of "predictive likelihood" [START_REF] Mathiasen | Prediction functions[END_REF], [START_REF] Bjornstad | Predictive likelihood: A review[END_REF], [START_REF] Bayarri | Difficulties and ambiguities in the definition of a likelihood function[END_REF]. For instance, the predictive profile likelihood is defined by L x (z) = sup θ f θ ,x (y, z). However, this notion is quite different conceptually from the standard notion of likelihood and, to some extent, arbitrary. While it does have interesting theoretical properties [START_REF] Mathiasen | Prediction functions[END_REF], its use poses some practical difficulties [6, page 39].

The method described in this paper builds upon the likelihood-based approach by seeing the likelihood function as describing the plausibility of each possible value of the parameter, in the sense of the Dempster-Shafer theory of belief functions [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF][10] [START_REF] Shafer | A mathematical theory of evidence[END_REF]. This approach of statistical inference was first proposed by Shafer [START_REF] Shafer | A mathematical theory of evidence[END_REF] and was later investigated by several authors (see, e.g., [START_REF] Wasserman | Belief functions and statistical evidence[END_REF] [START_REF] Aickin | Connecting Dempster-Shafer belief functions with likelihood-based inference[END_REF]). It was recently justified by Denoeux in [START_REF] Denoeux | Likelihood-based belief function: justification and some extensions to low-quality data[END_REF] and extended to prediction in [START_REF] Kanjanatarakul | Forecasting using belief functions: an application to marketing econometrics[END_REF][START_REF] Kanjanatarakul | Econometric forecasting using linear regression and belief functions[END_REF]. In this paper, we provide a general introduction to estimation and prediction using belief functions and demonstrate the application of this inference framework to the stochastic frontier model. In this model, the determination of the production frontier and disturbance parameters is an estimation problem, whereas the determination of the inefficiency terms is a prediction problem. We will show, in particular, how this method makes it possible to quantify both estimation uncertainty and random uncertainty, and to evaluate the plausibility of various hypothesis about both the production frontier and the efficiencies.

The rest of this paper is organized as follows. The general framework for inference and prediction will first be recalled in Section 2. This framework will be particularized to the stochastic frontier model in Section 3 and Section 4 will conclude the paper.

Inference and prediction using belief functions

Basic knowledge of the theory of belief functions will be assumed throughout this paper. A complete exposition in the finite case can be found in Shafer's book [START_REF] Shafer | A mathematical theory of evidence[END_REF]. The reader is referred to [START_REF] Ben Abdallah | Combining statistical and expert evidence using belief functions: Application to centennial sea level estimation taking into account climate change[END_REF] for a quick introduction on those aspects of this theory needed for statistical inference. In this section, the definition of a belief function from the likelihood function and the general prediction method introduced in [START_REF] Kanjanatarakul | Forecasting using belief functions: an application to marketing econometrics[END_REF] will be recalled in Sections 2.1 and 2.2, respectively.

Inference

Let f θ ,x (y) be the marginal probability mass or density function of the observed data Y given x. In the following, the covariates (if any) will be assumed to be fixed, so that the notation f θ ,x (y) can be simplified to f θ (y). Statistical inference has been addressed in the belief function framework by many authors, starting from Dempster's seminal work [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF]. In [START_REF] Shafer | A mathematical theory of evidence[END_REF], Shafer proposed, on intuitive grounds, a more direct approach in which a belief function Bel Θ y on Θ is built from the likelihood function. This approach was further elaborated by Wasserman [START_REF] Wasserman | Belief functions and statistical evidence[END_REF] and discussed by Aickin [START_REF] Aickin | Connecting Dempster-Shafer belief functions with likelihood-based inference[END_REF], among others. It was recently justified by Denoeux in [START_REF] Denoeux | Likelihood-based belief function: justification and some extensions to low-quality data[END_REF], from three basic principles: the likelihood principle, compatibility with Bayesian inference and the least commitment principle [START_REF] Ph | Belief functions: the disjunctive rule of combination and the generalized Bayesian theorem[END_REF]. The least committed belief function verifying the first two principles, according to the commonality ordering [START_REF] Dubois | A set-theoretic view of belief functions: logical operations and approximations by fuzzy sets[END_REF] is the consonant belief function Bel Θ y defined by the contour function

pl y (θ ) = L y (θ ) sup θ ∈Θ L y (θ ) , (4)
where L y (θ ) = f θ (y) is the likelihood function. The quantity pl y (θ ) is interpreted as the plaussibility that the true value of the parameter is θ . The corresponding plausibility and belief functions can be computed from pl y as:

Pl Θ y (A) = sup θ ∈A pl y (θ ), (5a)
Bel Θ y (A) = 1 -sup θ ∈A pl y (θ ), (5b)
for all A ⊆ Θ . The focal sets of Bel Θ y are the levels sets of pl y (θ ) defined as follows:

Γ y (ω) = {θ ∈ Θ |pl y (θ ) ≥ ω}, (6)
for ω ∈ [0, 1]. These sets may be called plausibility regions and can be interpreted as sets of parameter values whose plausibility is greater than some threshold ω. When ω is a random variable with a continuous distribution U ([0, 1]), Γ y (ω) becomes a random set equivalent to the belief function Bel Θ y , in the sense that

Bel Θ y (A) = P ω (Γ y (ω) ⊆ A) (7a) Pl Θ y (A) = P ω (Γ y (ω) ∩ A = / 0), (7b)
for all A ⊆ Θ such that the above expressions are well-defined.

(θ ) = (2π) -n/2 exp -1 2 ∑ n i=1 (y i -θ ) 2 (2π) -n/2 exp -1 2 ∑ n i=1 (y i -y) 2 (8a) = exp - n 2 (θ -y) 2 , ( 8b
)
where y is the sample mean. The plausibility and belief that θ does not exceed some value t are given by the upper and lower cumulative distribution functions (cdfs) defined, respectively, as

Pl y (θ ≤ t) = sup θ ≤t pl x (θ ) (9a) = exp -n 2 (t -y) 2 if t ≤ y 1 otherwise (9b) and Bel y (θ ≤ t) = 1 -sup θ >t pl x (θ ) (10a) = 0 if t ≤ y 1 -exp -n 2 (t -y) 2 otherwise. (10b)
The focals sets (6) are closed intervals

Γ y (ω) = y - -2 ln ω n , y + -2 ln ω n . ( 11
)
When ω has a uniform distribution on [0, 1], Γ y (ω) is a closed random interval. The cdfs of its lower and upper bounds are equal, respectively, to the lower and upper cdfs [START_REF] Dempster | A generalization of Bayesian inference (with discussion)[END_REF] and (9).

Prediction

The prediction problem can be defined as follows: having observed the realization y of Y with distribution f θ (y), we wish to make statements about some yet unobserved data Z ∈ Z whose conditional distribution f y,θ (z) given Y = y also depends on θ . The uncertainty on Z has two sources: (1) the randomness of the generation mechanism of Z given θ and y and (2) the estimation uncertainty on θ . In the approach outlined here, the latter uncertainty is represented by the belief function Bel Θ y on θ obtained by the approach described in the previous section. The random generation mechanism for Z can be represented by a sampling model such as the one used by Dempster [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF] for inference. In this model, the new data Z is expressed as a function of the parameter θ and an unobserved auxiliary random variable ξ with known probability distribution independent of θ :

Z = ϕ(θ , ξ ), ( 12
)
where ϕ is defined in such a way that the distribution of Z for fixed θ is f y,θ (z).

When Z is a real random variable, a canonical model of the form ( 12) can be obtained as Z = F -1 y,θ (ξ ), where F y,θ is the conditional cumulative distribution function (cdf) of Z given Y = y, F -1 y,θ is its generalized inverse and ξ has a continuous uniform distribution in [0, 1]. This canonical model can be extended to the case where Z is a random vector. For instance, assume that Z is a two-dimensional random vector (Z 1 , Z 2 ). We can write

Z 1 = F -1 y,θ (ξ 1 ) (13a)
Z 2 = F -1 y,θ ,Z 1 (ξ 2 ), (13b)
where F y,θ is the conditional cdf of Z 1 given Y = y, F y,θ ,Z 1 is the conditional cdf of Z 2 given Y = y and Z 1 and ξ = (ξ 1 , ξ 2 ) has a uniform distribution in [0, 1] 2 . Equation [START_REF] Dubois | A set-theoretic view of belief functions: logical operations and approximations by fuzzy sets[END_REF] gives us the distribution of Z when θ is known. If we only know that θ ∈ Γ y (ω) and the value of ξ , we can assert that Z is in the set ϕ(Γ y (ω), ξ ).

As ω and ξ are not observed but have a joint uniform distribution on [0, 1] 2 , the set ϕ(Γ y (ω), ξ ) is a random set. It induces belief and plausibility functions defined as

Bel y (Z ∈ A) = P ω,ξ (ϕ(Γ y (ω), ξ ) ⊆ A) , (14a)
Pl y (Z ∈ A) = P ω,ξ (ϕ(Γ y (ω), ξ ) ∩ A = / 0) , (14b)
for any A ⊆ Z.

Example 2 Continuing Example 1, let us assume that Z ∼ N (θ , 1) is a yet unoberved normal random variable independent of Y . It can be written as

Z = θ + Φ -1 (ξ ), ( 15
)
where Φ is the cdf of the standard normal distribution. The random set ϕ(Γ y (ω), ξ ) is then the random closed interval

ϕ(Γ y (ω), ξ ) = y - -2 ln ω n + Φ -1 (ξ ), y + -2 ln ω n + Φ -1 (ξ ) . (16)
Expressions ( 14) for the belief and plausibility of any assertion about Z can be approximated by Monte Carlo simulation [START_REF] Kanjanatarakul | Forecasting using belief functions: an application to marketing econometrics[END_REF].

As remarked by Bjornstad [START_REF] Bjornstad | Predictive likelihood: A review[END_REF], a prediction method should have at least two fundamental properties: it should be invariant to any one-to-one reparametrization of the model and it should be asymptotically consistent, in a precise sense to be defined. An additional property that seems desirable is compatibility with Bayesian inference, in the sense that it should yield the same result as the Bayesian approach when a prior distribution on the parameter is provided. Our method possesses these three properties. Parameter invariance follows from the fact that it is based on the likelihood function; compatibility with Bayes is discussed at length in [START_REF] Kanjanatarakul | Forecasting using belief functions: an application to marketing econometrics[END_REF] and consistency will be studied in greater detail in a forthcoming paper.

Application to stochastic frontier analysis

In this section, we apply the above estimation and prediction framework to the stochastic frontier model (SFM). To keep the emphasis on fundamental principles of inference, only the simplest case of cross-sectional data will be considered. The model as well as the inference method will be in introduced in Section 3.1 and an illustration with simulated data will be presented in Section 3.2.

Model and inference

The SFM [START_REF] Aigner | Formulation and estimation of stochastic frontier production function models[END_REF] defines a production relationship between a p-dimensional input vector x i and output Y i of each production unit i of the form

lnY i = β ln x i +V i -U i , ( 17
)
where β is a vector of coefficients, V i is an error term generally assumed to have a normal distribution N (0, σ 2 v ) and U i is a positive inefficiency term. Usual models for U i are the half-normal distribution |N (0, σ 2 u )| (i.e., the distribution of the absolute value of a normal variable) and the exponential distribution. The SFM is thus a linear regression model with asymmetric disturbances ε i = V i -U i . The inefficiency terms U i are not observed but are of particular interest in this setting.

Assuming U i to have a half-normal distribution, let λ = σ u /σ v and σ 2 = σ 2 u + σ 2 v be new parameters to be used in place of σ 2 u and σ 2 v . Although the variance of U i is not σ 2 u but (1 -2/π)σ 2 u , λ has an intuitive interpretation as the relative variability of the two sources of error that distinguish firms from one another [START_REF] Aigner | Formulation and estimation of stochastic frontier production function models[END_REF]. Using the notations defined in Section 1, we have Y = (Y 1 , . . . ,Y n ), Z = (U 1 , . . . ,U n ) and θ = (β , σ , λ ). The determination of the inefficiency terms is thus a prediction problem.

Parameter estimation

Assuming the two error components U i and V i to be independent, the log-likelihood function is [14, page 540]

ln L y (θ ) = -n ln σ + n 2 log 2 π - 1 2 n ∑ i=1 ε i σ 2 + n ∑ i=1 ln Φ - ε i λ σ . ( 18
)
The maximum likelihood estimate (MLE) θ can be found using an iterative nonlinear optimization procedure. Parameter β may be initialized by the least squares estimate, which is unbiased and consistent (except for the constant term) [START_REF] Greene | Econometric analysis[END_REF]. However, it may be wise to restart the procedure from several randomly chosen initial states, as the log-likelihood function may have several maxima for this problem.

Once θ has been found, the contour function ( 4) can be computed. The marginal contour function for any subset of parameters is the relative profile likelihood function. For instance, the marginal contour function of λ is

pl y (λ ) = sup β ,σ pl y (θ ). (19)

Prediction

The main purpose of stochastic frontier analysis is the determination of the inefficiency terms u i , which are not observed. The usual approach is to approximate u i by E(U i |ε i ), which is itself estimated by plugging in the MLEs and by replacing ε i by the residuals ε i . The main result is due to Jondrow et al. [START_REF] Jondrow | On the estimation of technical efficiency in the stochastic production function model[END_REF], who showed that the conditional distribution of U given ε i , in the half-normal case, is that of a normal N (µ * , σ 2 * ) variable truncated at zero, with

µ * = - σ 2 u ε i σ 2 = - ε i λ 2 1 + λ 2 (20a) σ * = σ u σ v σ = λ σ 1 + λ 2 . ( 20b
)
The conditional expectation of U i given ε i is

E(U i |ε i ) = λ σ 1 + λ 2 φ (λ ε i /σ ) 1 -Φ(λ ε i /σ ) - λ ε i σ , (21)
where φ and Φ are, respectively, the pdf and cdf of the standard normal distribution. As noted by Jondrow et al. [START_REF] Jondrow | On the estimation of technical efficiency in the stochastic production function model[END_REF], when replacing the unknown parameter values by their MLEs, we do not take into account uncertainty due to sampling variability. While this uncertainty becomes negligible when the sample size tends to infinity, it certainly is not when the sample is of small or moderate size.

To implement the approach outlined in Section 2.2 for this problem, we may write the cdf of U i as

F(u) = Φ[(u -µ * )/σ * ] -Φ(-µ * /σ * ) 1 -Φ(-µ * /σ * ) 1 [0,+∞) (u). (22)
Let ξ i = F(U i ), which has a uniform distribution U ([0, 1]). Solving the equation ξ i = F(U i ) for U i , we get

U i = µ * + σ * Φ -1 ξ i 1 -Φ - µ * σ * + Φ - µ * σ * . (23)
Replacing µ * and σ * by their expressions as functions of the parameters, we have

U i = ϕ(θ , ξ i ) = λ 1 + λ 2 -ε i λ + σ Φ -1 ξ i + Φ ε i λ σ (1 -ξ i ) (24)
with ε i = ln y i -β ln x i , which gives us an equation of the same form as [START_REF] Dubois | A set-theoretic view of belief functions: logical operations and approximations by fuzzy sets[END_REF], relating the unobserved random variable U i to the parameters and the auxiliary variable ξ i .

To approximate the belief function on Z = (U 1 , . . . ,U n ), we may use the Monte Carlo method described in [START_REF] Kanjanatarakul | Forecasting using belief functions: an application to marketing econometrics[END_REF]. More specifically, we randomly generate N n + 1tuples (ω ( j) , ξ

i ] be the resulting interval. The belief and plausibility of any statement Z ∈ A for A ⊂ R n as defined by ( 14) can be approximated by

Bel y (Z ∈ A) ≈ 1 N # j ∈ {1, . . . , N}|[u ( j) 1 , u ( j) 1 ] × . . . × [u ( j) n , u ( j) n ] ⊆ A , (26a)
Pl y (Z ∈ A) ≈ 1 N # j ∈ {1, . . . , N}|[u ( j) 1 , u ( j) 1 ] × . . . × [u ( j) n , u ( j) n ] ∩ A = / 0 , (26b)
where # denotes cardinality. We can also approximate the belief function on any linear combination ∑ n i=1 α i u i by applying the same transformation to the intervals [u ( j) i , u ( j) i ], using interval arithmetics. For example, the belief and plausibility of statements of the form u iu k ≤ c can be approximated as follows:

Bel y (u i -u k ≤ c) ≈ 1 N # j ∈ {1, . . . , N}|u ( j) i -u ( j) k ≤ c , (27a) Pl y (u i -u k ≤ c) ≈ 1 N # j ∈ {1, . . . , N}|u ( j) i -u ( j) k ≤ c . (27b)

Simulation experiments

To illustrate the behavior of our method, we simulated data from model with p = 1, β = (1, 0.5) , σ v = 0.175 and σ u = 0.3. We thus have, for this model, λ = 1.7143 and σ = 0.3473. Figure 1 displays the marginal contour functions of β 0 , β 1 , σ and λ for a simulated sample of size n = 100. These plots show graphically the plausibility of any assertion of the form θ j = θ j0 . For instance, we can see from Figure 1(d) that the plausibility of the assertion λ = 0 is around 0.6: consequently, the hypothesis that inefficiencies are all equal to zero is quite plausible, given the data. Figures 2 and3 show the true and estimated inefficiencies for 20 individuals in the above simulated sample of size n = 100 and in a simulated sample of size n = 1000, respectively. For the belief function estimation, we give the quantile intervals for α = 5% and α = 25%. The lower bound of the quantile interval [START_REF] Kanjanatarakul | Forecasting using belief functions: an application to marketing econometrics[END_REF] is the α quantile of the lower bounds u ( j) i of the prediction intervals, while the upper bound is the 1-α quantile of the upper bounds u ( j) i . The larger intervals in the case n = 100 reflect the higher estimation uncertainty.

Conclusions

We have shown how the estimation and prediction problems may be solved in the belief function framework, and illustrated these solutions in the case of the stochastic frontier model with cross-sectional observation. In the case of this model, the estimation problem concerns the determination of the model parameters describing the production frontier and the distributions of the noise and inefficiency terms, while the prediction problem consists in the determination of the unobserved inefficiency terms, which are of primary interest in this analysis. In our approach, uncertainties about the parameters and the inefficiencies are both modeled by belief functions induced by random sets. In particular, the random set formulation allows us to approximate the belief or plausibility of any assertion about the inefficiencies, using Monte Carlo simulation.

We can remark that parameters and realizations of random variables (here, inefficiencies) are treated differently in our approach, whereas there are not in Bayesian inference. In particular, the likelihood-based belief functions in the parameter space are consonant, whereas predictive belief functions are not. This difference is not due to conceptual differences between parameters and observations, which are just considered here as unknown quantities. It is due to different natures of the evidence from which the belief functions are constructed. In the former case, the evidence consists of observations that provide information on parameters governing a random process. In the latter case, evidence about the data generating process provides information about unobserved observations generated from that process.

The evidential approach to estimation and prediction outlined in this paper is invariant to one-to-one transformations of the parameters and compatible with Bayesian inference, in the sense that it yields the same result when provided with the same initial information. It is, however, more general, as it does not require the user to supply prior probability distributions. It is also easily implemented and does not require asymptotic assumptions, which makes it readily applicable to a wide range of econometric models.

The preliminary results reported in this paper need to be completed in several ways. First, a detailed comparison, based on underlying principles, with alternative approaches such as, e.g., the empirical Bayes method [START_REF] Robbins | An empirical Bayes approach to statistics[END_REF] or imprecise probabilities [START_REF] Walley | Statistical Reasoning with Imprecise Probabilities[END_REF] remains to be performed. Secondly, it would be interesting to study experimentally how users interpret the results of the belief function analysis to make decisions in real-world situations. Finally, theoretical properties of our method, such as asymptotic consistency, are currently being studied.

1 )

 1 for j = 1, . . . , N uniformly in [0, 1] n+1 . For i = 1 to n and j = 1 to N, we compute the minimum and the maximum of ϕ(θ , ξ ( j) i ) w.r.t. θ under the constraint pl y (θ ) ≥ ω ( j) .

Fig. 1

 1 Fig. 1 Marginal contour functions for a simulated sample of size n = 100.

Fig. 2

 2 Fig. 2 True and predicted inefficiencies for 20 individuals in a simulated sample of size n = 100.

Fig. 3

 3 Fig. 3 True and predicted inefficiencies for 20 individuals in a simulated sample of size n = 1000.

  Example 1 Let us consider the case where Y = (Y 1 , . . . ,Y n ) is an i.i.d. sample from a normal distribution N (θ , 1). The contour function on θ given a realization y of Y is pl y


Juliette Pierre

Benjamin Dollet

Valentin Leroy

Resonant Acoustic Propagation and Negative Density in Liquid Foams

We measured the dispersion relation for acoustic longitudinal waves in liquid foams, over a broad frequency range (60-600 kHz). Strong dispersion was found, with two non-dispersive behaviors, separated by a negative density regime. A new model, based on the coupled displacements of films, liquid channels and gas in the foam, rationalizes all the experimental findings.

Liquid foams, dispersions of gas bubbles in a liquid matrix stabilized by surfactants [1], are present in a wide range of industrial applications, from food and personal care to ore flotation and enhanced oil recovery [START_REF]Foam Engineering: Fundamentals and Applications[END_REF]. They are opaque materials, hence difficult to characterize in depth, and new ways of probing them are highly desirable. Acoustic probes are good candidates, since they are unexpensive and nonintrusive, and widely used as such in nondestructive testing [START_REF] Blitz | Ultrasonic Methods of Nondestructive Testing[END_REF]. However, surprisingly little is known about liquid foam acoustics, contrary to other multiphasic media such as porous media [START_REF] Attenborough | [END_REF], colloidal suspensions, or emulsions [5,6]. Most of the few existing experimental studies reported speeds of sound of order 50 m/s [7][8][9], close to the so-called Wood model [START_REF] Wood | A Textbook of Sound[END_REF]. The latter treats foams as an effective medium, which density and compressibility are given by the mixture law, i.e. averages of those of the gas and liquid phases, weighted by their respective volume fractions. However, much higher speeds of sound, of order 200 m/s, were also measured [START_REF] Moxon | [END_REF]12], and some studies reported a resonant behavior [13,14], reminiscent of the low frequency resonance of a single bubble in an unbounded liquid, the so-called Minnaert resonance [15], which is the key ingredient of acoustic propagation through dilute bubbly liquids [16,17]. In this Letter, we use a novel setup [18] to measure speed of sound and attenuation through liquid foams over a large range of frequencies f (60-600 kHz) and bubble radii (15-50 µm). Our results are explained by a new model, which fully reconciles all aforementioned different viewpoints on liquid foam acoustics. In particular, we show that liquid foams are natural acoustic metamaterials, exhibiting a negative effective density over a large range of frequencies and bubble sizes.

The general principle of our technique is to measure the complex transmission of short ultrasonic pulses through a foam sample of known thickness. From this complex transmission, we can determine the effective complex wavevector as a function of frequency: k = ω/v + iα, with ω = 2πf , and v and α the phase velocity and attenuation of sound. In an ageing foam, the bubble median radius R increases over time by gas diffusion between bubbles [1]. Hence, making measurements at different times, we obtain k(ω, R), the complex wavevector as a function of frequency and median radius. Some aspects of the ultrasonic technique and the foam production and characterization are briefly described below; more details can be found in [18] for the former, and [19] for the latter.

Transmission measurements were performed with two broadband air transducers. The thickness of the foam samples was set to 0.5 mm by sandwiching them between two plastic films, which were as thin as possible to allow a good transmission of ultrasound. An inversion procedure was used to deduce the effective wavevector k of sound in the foam from the complex transmission measured through the three-layer system {film-foam-film}. This procedure relied on the assumption that either the effective density ρ eff or the effective compressibility χ eff were known. Guided by the model presented below, we assumed that the effective compressibility of the liquid foam was given by the usual mixture law:

χ eff = Φχ ℓ + (1 -Φ)χ g ,
where Φ is the liquid volume fraction and χ ℓ,g the compressibility of the liquid and gaseous phase.

Liquid foams with liquid volume fractions ranging from 3% to 22% were obtained by the two-syringe method [19].

As in [19], the foaming liquid was an aqueous solution containing 10 g/L of sodium dodecyl sulfate to ensure good foamability, and 0.5 g/L of xanthane to reduce drainage. Given the small thickness of the samples (0.5 mm), the typical radius of the bubbles (less than 50 µm) and the presence of xanthane, we were able to neglect the gradient of liquid fraction due to gravity [20]. The gaseous phase was air saturated with vapor of C 6 F 14 , an insoluble gas that slows down ageing. Thus, the typical time over which the foam was changing (∼ 10 minutes) was long enough not only to consider the acoustic measurements as instantaneous, but also to perform bubble size measurements [19]. To do that, we took a small sub-sample of the foam from the syringe, and pour it on a liquid bath to obtain a bubble raft, a two-dimensional structure easy to image. Assuming that the left-over in the syringe was ageing the same way as the foam between the films, we were able to follow the time evolution of the bubble size distribution. We found lognormal distributions exp{-[ln(r/R)] 2 /2ǫ 2 }/(rǫ √ 2π) with a median radius R going from 15 to 50 µm in 90 minutes, and a polydispersity of ǫ ≃ 0.4 without significant time evolution (Fig. 1c). Measurements shown in Fig. 1a clearly indicate that the propagation of acoustic waves in a liquid foam is dispersive. The phase velocity at low frequencies is close to the value predicted by Wood law:

v = [(Φρ ℓ + (1 - Φ)ρ g )(Φχ ℓ + (1 -Φ)χ g )] -1/2 =
32 m/s for Φ = 11%. But it significantly deviates from this value as the frequency increases and the foam ages. At t = 90 min the phase velocity is of the order of 220 m/s.

The attenuation (Fig. 1b) is also frequency-dependent and evolves as the foam ages. For the longest times (i.e. largest bubbles), a peak is clearly visible. One can report the frequency of the maximum of attenuation f max as a function of R (Fig. 1d), and compare it with the Minnaert frequency f M = 3γP 0 /ρ/(2πR) [15], where P 0 is the pressure of the gas in the bubbles, and γ the ratio of the specific heat capacities. Taking f max = 1.6f M gives a reasonable law (dashed line in Fig. 1d), suggesting an effective Minnaert resonance, as proposed in previous experiments [13,14]. However, the Minnaert frequency is calculated for a single bubble in an infinite volume of liquid of density ρ, and for highly concentrated media such as foams, it seems unrealistic. Replacing ρ by the actual density of the surrounding medium Φρ would lead to 3f M for Φ = 11%, larger than the 1.6 factor we measured. Moreover, it predicts a resonance frequency that decreases with Φ, contrary to our experiments (Fig. 2). Finally, the range of radii over which we measured a maximum of attenuation (Fig. 1d) is too small to validate a scaling law in f max ∼ 1/R. Instead of limiting our analysis to the peak of attenuation, we can take advantage of all the radii by collapsing the data on a master curve: we plot kR p as a function of f R p , with an exponent p. As shown in Fig 3 cates that R is the key parameter that governs the frequency dependence of the acoustic behavior of the foam (the polydispersity being almost constant with time). Contrary to usual bubbly media, bubbles in a foam are in contact through thin films, separated by liquid channels (also called Plateau borders) and vertices [1]. Usually films contain a very small amount of water, hence films and liquid channels have very different inertia, whose influence we shall now model. In first approximation, we neglect polydispersity, and we idealize the foam as a one-dimensional array of spatial periodicity d of the following unit cell (Fig. 4a): a flexible, circular membrane of thickness e and radius a (representing the soap film), attached to a rigid ring of external radius b and mass m c (representing the liquid channel), and surrounded by air of density ρ a . Experiments show that the wavelength λ = 2π/Re(k) remains larger than 7R, even at the highest frequency investigated (Fig. 3). Hence, the model is based on the hypothesis λ ≫ d. Note that the deformations are not as high as suggested by the picture: the maximum of displacement of the film at 800 kHz for instance is about 1.5 nm, i.e. more than four orders of magnitude smaller than the diameter of the film.

To predict the acoustic wavevector, we need to model the behavior of the unit cell under compression (effective compressibility χ eff ), as well as its motion driven by the average pressure gradient (P 3 -P 1 )/d (effective density ρ eff ): k 2 = ω 2 ρ eff χ eff . Our model geometry resembles that of Lee et al. [21], who studied sound propagation through an array of elastic membranes. In particular, in the limit λ ≫ d, the effective compressibility is the same: an almost uniform pressure change is applied to the unit cell, hence its total volume changes as the sum of the liquid volume change and the gas volume change, which yields the usual mixture law for the effective compressibility.

The effective density, on the other hand, is more complicated because the flexible film and the rigid ring do not react the same way when submitted to a given pressure gradient. A portion of the film between r and r+dr is driven by the pressure force 2πrdr(P 1 -P 2 ) (see Fig. 4a) and by the vertical projection of the tension force 4πσ[(r∂h f /∂r) r+dr -(r∂h f /∂r) r ]. In harmonic regime, h f (r, t) = Re[z f (r)e -iωt ], the profile of the film is thus given by z f (r) = AJ 0 (qr) + ∆P/(ρeω 2 ), where A is its amplitude of vibration, ∆P = P 2 -P 1 the driving pressure, q = ω ρe/2σ the wavenumber of capillary waves on the film [? ], and J 0 the Bessel function of order 0 of the first kind [START_REF]Fundamental of Acoustics[END_REF]. The unknown constant A is determined by matching the motion of the film edge with that of the ring: h f (r = a, t) = h c (t). It is interesting to note that even though the pressure field is uniform in the lateral direction, it produces oscillations of the film. This is due to the film being attached to the ring, which has a different inertia. The ring is driven by the pressure difference ∆P and by the traction of the film -4πσa(1iωτ )(∂z f /∂r) r=a , in which we include a phenomenological damping time τ to account for a friction force (see SM1 for a full discussion on dissipation). We finally get the full profile

z f (r) = ∆P ρeω 2 1 - m c -m f (b 2 -a 2 )/a 2 m c + m f (1 -iωτ )H(qa) J 0 (qr) J 0 (qa) , (1)
where m f = πa 2 eρ is the mass of the film and H(qa) = 2J 1 (qa)/[qaJ 0 (qa)].

We now consider the motion of air between two consecutive films, which is driven by the pressure difference P 2 -P 3 . As λ ≫ d, we may consider that the whole air of the unit cell moves as one block with displacement z a given by -m a ω 2 z a = (P 2 -P 3 )πb 2 , where m a is the mass of the air in the cell. Moreover, the air in contact with the film and ring is also entrained by their motion. By continuity of displacement, we simply assume that z a equals the average displacement of the film and ring:

z a = (1 -x)z c + 2x a 2 a 0 rz f (r)dr, where x = a 2 /b 2
is the fraction of the surface covered by the film. We can then calculate the average displacement z a of the unit cell, and relate it to the effective density ρ eff by ω 2 ρ eff z a = (P 3 -P 1 )/d. Finally, we obtain an effective density of the foam ρ eff = (1 -Φ)ρ a + Φ ′ ρ, where the frequency-dependent effective liquid volume fraction Φ ′ is given by

Φ ′ = Φ c + Φ f (1 -iωτ )H(qa) 1 + x 2 Φ f +Φc Φ f -2x [1 -H(qa)] -iωτ xH(qa) , (2)
with Φ c and Φ f the volume fraction of liquid contained in the channels and in the films, respectively: Φ c + Φ f = Φ. Eq. ( 2) predicts two asymptotic non-dispersive behaviors. At low frequencies, H(qa) ≃ 1, hence Re Φ ′ = Φ: the mixture law for ρ eff is recovered. At high frequencies, H(qa) tends towards 0, leading to Re (1

/Φ ′ ) = (1 -x) 2 /Φ c + x 2 /Φ f . If Φ f ≪ Φ c [x/(1 -x)] 2 (i.e.
films contain a negligible part of water, but with a nonnegligible surface) this last formula reduces to x 2 /Φ f : the effective liquid density is governed by the films, as in Kann's model [12]. The transition between the two asymptotic regimes can be inspected by noting that 1 -H(qa) ∼ 0 -(qa) 2 /8. Thus, within the approximation

Φ c ≫ Φ f |1 -2/x|, Eq. (2) becomes Φ ′ ≃ Φ c /[1 -(ω/ω 0 ) 2 -ixωτ ],
similar to the response function of an harmonic oscillator with a resonance frequency ω 2 0 = 16σΦ f /(ρa 2 ex 2 Φ) and a damping factor xωτ . Interestingly, this resonance frequency depends neither on a nor on e since Φ f = n f πa 2 e, n f being the number per unit volume of films perpendicular to the direction of propagation. If we consider that there are N such films per bubble, we obtain n f = 3N (1 -Φ)/(4πR 3 ) and the resonance frequency becomes

ω 2 0 = 12N σ(1 -Φ) x 2 ρΦR 3 , (3)
which predicts the experimentally observed R -1.5 scaling law. The mechanism of this resonance is different from Minnaert's one: inertia also comes from water, but here the restoring force is due to the tension of the film, not the compressibility of air.

To go beyond and get a quantitative prediction over the full range of frequencies and sizes, we refine the model and include polydispersity, as explained in SM1. Briefly, this amounts to changing H(qa) by I = H(qa)a 2 n(a)da/ a 2 n(a)da in the expression (2) of the effective liquid fraction, with n(a) the distribution of film radii. If we assume that this distribution is lognormal with a median radius a 0 , a polydispersity ǫ f and a total number of films per unit volume n f , the model counts a total of nine parameters. Three are known: σ = 35 mN/m, Φ c = 11%, ǫ = 0.4. Four can be estimated: x, n f , a 0 and ǫ f . The surface fraction covered by films in a liquid foam was studied by Princen [24]. He found the following empirical dependence on Φ: x = 1 -3.20 [7.70 + (1 -Φ)/Φ] -1/2 , which gives x = 0.2 in our case. Then one can estimate that n f = 1620 films per mm 3 (assuming N = 1), a 0 = 15 µm (taking an average of 6 films per bubble imposes x = 1.5(a 0 /R) 2 ) and ǫ f = 0.4 (same polydispersity as for bubble radii). The two remaining parameters were fitted to the experimental data: τ governs the width of the resonance, e the high-frequency effective density. With τ = 10 µs and e = 70 nm, the model agrees well with the experimental data (Fig. 5). This value of e is compatible with the usual measurements of film thickness, of a few tens of nanometers [1]. It also justifies the Taylor expansion leading to Eq. ( 3), since it leads to (qa) 2 ≃ 0.15 at the resonance frequency (130 kHz). It is less easy to compare τ to any standard measurement, and the local study of the microscopic dissipation mechanisms at high frequency is left as a perspective of our work.

The model also captures well the dependence on Φ of the frequency of the maximum of attenuation, as shown in Fig. 2. The polydispersity has to be included to obtain a quantitative comparison (see SM1), but Eq. ( 3) is sufficient to understand the non-monotonic behavior: for an increasing liquid fraction, the increase of density (1/ √ Φ term, as in the modified-Minnaert relation) is dominated by the shrinkage of the films (1/x term).

A salient feature of our study is the existence of a negative real part of the density over an extended range of frequencies (100 to 300 kHz). The mechanism of this negative effective density is well illustrated by looking at the displacement of a film and a liquid channel predicted by (1) (see Fig. 4c). At low frequencies, they move in phase with the incoming pressure. At high frequencies, only the film moves. At intermediate frequencies, while the liquid channel has a small in-phase displacement, the film moves out of phase with an amplitude large enough to compensate for its small inertia. The net average movement of the system is then out of phase, which leads to a negative effective density [? ]. As a conclusion, we have evidenced two regimes for the sound propagation in a liquid foam, separated by a resonance. Our study reconciles the seemingly contradictory results previously reported: a low speed of sound, compatible with Wood's model, at low frequency and small bubble size; a much higher one, slightly lower than that in air, at high frequency and large bubble size; and a resonance in between, with a maximum of attenuation. This rich variety of behaviors is fully captured by a model coupling the motion of air, films and liquid channels.

Finally, we have shown for the first time that liquid foams are naturally acoustic single-negative metamaterials, with an effective negative density over an extended range of frequencies and bubble sizes. Moreover, in contrast with most existing metamaterials [25,26], which are synthetic, highly ordered materials, this property exists even though the foam is isotropic and polydisperse. Our findings could thus bring fresh insight to design metamaterials for applications such as acoustic insulation.

We acknowledge funding support from the Agence Nationale de la Recherche (ANR-11-BS09-001). We thank

FIG. 1 :

 1 FIG. 1: (a) Real part of the wavenumber as a function of the frequency, measured at three different ageing times in a liquid foam with Φ = 11%. Lines show dispersion curves for constant velocities of 32 m/s (dash) and 220 m/s (solid). (b) Imaginary part of the wavenumber. (c) Median radius of the bubbles (•) and polydispersity (•) as functions of the ageing time. (d) Frequency of the maximum of attenuation as a function of the median radius; solid line is the Minnaert frequency fM , dashed line is 1.6fM .

FIG. 2 :

 2 FIG.2: Frequency of the maximum of attenuation as a function of the liquid fraction, for median radii of 40 µm. Dashed line is the Minnaert resonance frequency with an effective density: fM / √ Φ. Solid line is the prediction of our model.

FIG. 3 :

 3 FIG. 3: Rescaled dispersion relation for a reference median radius of R0 = 40 µm. Lines are for constant velocities of 32 (dash) and 220 m/s (solid). Inset: the R -1.5 scaling (solid symbols) is better than the R -1 one (open symbols).

FIG. 4 :

 4 FIG. 4: (a) Sketch of the idealized foam structure: a periodic array of flexible films and rigid rings (representing the liquid channels), separated by a distance d. (b) Zoom on one unit cell: the flexible film is attached to a rigid ring sketched as a square. (c) Profiles of a 20 µm radius film at three different frequencies, as predicted by the model for a ∆P = 1 Pa excitation.Note that the deformations are not as high as suggested by the picture: the maximum of displacement of the film at 800 kHz for instance is about 1.5 nm, i.e. more than four orders of magnitude smaller than the diameter of the film.

FIG. 5 :

 5 FIG. 5: Comparison between the measured effective complex density (• real part, • imaginary part) and prediction of the equivalent of Eq. (2) in the polydisperse model (see SM1).

J.-C. Bacri, C. Gay, F. Graner and A. Saint-Jalmes for fruitful discussions.