Dissimilarity Metric Learning in the Belief Function Framework

Chunfeng Lian, Su Ruan * , and Thierry Denoeux

Abstract-The Evidential K-Nearest-Neighbor (EK-NN) method provided a global treatment of imperfect knowledge regarding the class membership of training patterns. It has outperformed traditional K-NN rules in many applications, but still shares some of their basic limitations, e.g., 1) classification accuracy depends heavily on how to quantify the dissimilarity between different patterns and 2) no guarantee for satisfactory performance when training patterns contain unreliable (imprecise and/or uncertain) input features. In this paper, we propose to address these issues by learning a suitable metric, using a low-dimensional transformation of the input space, so as to maximize both the accuracy and efficiency of the EK-NN classification. To this end, a novel loss function to learn the dissimilarity metric is constructed. It consists of two terms: the first one quantifies the imprecision regarding the class membership of each training pattern; while, by means of feature selection, the second one controls the influence of unreliable input features on the output linear transformation. The proposed method has been compared with some other metric learning methods on several synthetic and real data sets. It consistently led to comparable performance with regard to testing accuracy and class structure visualization.

Index Terms-Dempster-Shafer Theory, Dissimilarity Metric Learning, Evidential K-NN, Feature Transformation, Feature Selection, Pattern Classification, Dimensionality Reduction.

I. INTRODUCTION

T HE K-nearest neighbor (K-NN) rule [1] is one of the most well-known pattern classification algorithms. As a case-based learning method without need of any prior assumptions [START_REF] Denoeux | Classification using belief functions: relationship between case-based and model-based approaches[END_REF], the K-NN classifier has been widely used in practice thanks to its simplicity. The original voting K-NN [START_REF] Cover | Nearest neighbor pattern classification[END_REF] assigns an object into the class represented by its majority nearest neighbors in the training set, while the information concerning the dissimilarity (distance) between the object and its neighbors is neglected. Then, the weighted K-NN [START_REF] Dudani | The distance-weighted k-nearest-neighbor rule[END_REF] has been proposed, in which this dissimilarity is imported into the classification procedure. However, in the case of uncertain and imprecise data, many samples may be corrupted with noise or located in highly overlapping areas; consequently, it becomes difficult for these classical K-NN classifiers to obtain satisfactory classification results.

Asterisk indicates corresponding author. Chunfeng Lian is with the Sorbonne Universités, Université de Technologie de Compiègne, CNRS, Heudiasyc, Compiègne, France, and also with the Université de Rouen, QuantIF, LITIS, Rouen, France (e-mail: chunfeng.lian@utc.fr).

* Su Ruan is with the Université de Rouen, QuantIF, LITIS, Rouen, France (e-mail: su.ruan@univ-rouen.fr).

Thierry Denoeux is with the Sorbonne Universités, Université de Technologie de Compiègne, CNRS, Heudiasyc, Compiègne, France (e-mail: thierry.denoeux@hds.utc.fr).

The Dempster-Shafer theory (DST), also known as the theory of belief functions or Evidence theory, is an extension of both probability theory and the set-membership approach [START_REF] Shafer | A mathematical theory of evidence[END_REF]. As a powerful framework for modeling and reasoning with uncertain and/or imprecise information [START_REF] Yager | On the Dempster-Shafer framework and new combination rules[END_REF]- [START_REF] Smets | The transferable belief model[END_REF], DST has shown remarkable applications in divers fields, such as model parameter estimation [START_REF] Denoeux | Maximum likelihood estimation from uncertain data in the belief function framework[END_REF]- [START_REF] Su | Parametric regression analysis of imprecise and uncertain data in the fuzzy belief function framework[END_REF], unsupervised learning [START_REF] Denoeux | EVCLUS: evidential clustering of proximity data[END_REF]- [START_REF] Makni | Introducing spatial neighbourhood in evidential C-means for segmentation of multi-source images: application to prostate multi-parametric mri[END_REF], supervised learning [START_REF] Denoeux | A K-nearest neighbor classification rule based on Dempster-Shafer theory[END_REF]- [START_REF] Jiao | Belief rulebased classification system: Extension of FRBCS in belief functions framework[END_REF], ensemble learning [START_REF] Altınc | Ensembling evidential k-nearest neighbor classifiers through multi-modal perturbation[END_REF], [START_REF] Bi | The combination of multiple classifiers using an evidential reasoning approach[END_REF], information fusion [START_REF] Milisavljevic | Sensor fusion in anti-personnel mine detection using a two-level belief function model[END_REF]- [START_REF] Lelandais | Fusion of multi-tracer PET images for dose painting[END_REF], etc. To endow the K-NN method with the capability to handle uncertain information, Denoeux has extended it in the belief function framework. An Evidential K-NN (EK-NN) rule has been proposed in [START_REF] Denoeux | A K-nearest neighbor classification rule based on Dempster-Shafer theory[END_REF], and further optimized in [START_REF] Zouhal | An evidence-theoretic K-NN rule with parameter optimization[END_REF]. The EK-NN rule provides a global treatment of partial knowledge regarding the class membership of training patterns. Ambiguity and distance reject options are also taken into account based on the concepts of lower and upper expected losses [START_REF] Denoeux | Analysis of evidence-theoretic decision rules for pattern classification[END_REF].

The EK-NN method has outperformed other traditional K-NN methods in many situations when using the same information [START_REF] Zouhal | An evidence-theoretic K-NN rule with parameter optimization[END_REF], whereas they still have some identical features: 1) the performances of the K-NN rules are strongly influenced by the chosen dissimilarity between different patterns. Better than directly using the simple Euclidean distance measure (such as in the original EK-NN), an adaptive dissimilarity metric tailored for the application should ensure better classification performance; 2) the efficiency of the K-NN rules substantially decrease when the dimensionality of the input data increases.

We propose a solution based on dissimilarity metric learning to deal with these inherent drawbacks of the K-NN classifications. Given an input space X, the metric learning problem can be formulated as finding a transformation matrix A, such that the dissimilarity between any two patterns can be defined in the transformed space Z = AX [START_REF] Yang | Distance metric learning: A comprehensive survey[END_REF]. Various studies have demonstrated that a properly learnt dissimilarity measure can dramatically boost the performance of the distance-based learning methods [START_REF] Lowe | Similarity metric learning for a variable-kernel classifier[END_REF]- [START_REF] Kedem | Non-linear metric learning[END_REF]. Even with a linear transformation of the input space [START_REF] Goldberger | Neighbourhood components analysis[END_REF]- [START_REF] Weinberger | Distance metric learning for large margin nearest neighbor classification[END_REF], the K-NN classification can reach significant improvement. In [START_REF] Goldberger | Neighbourhood components analysis[END_REF], Goldberger et al. proposed a metric learning method called Neighborhood Component Analysis (NCA), which maximizes the expected leave-oneout classification accuracy from a stochastic version of the K-NN classification. Based on a softmax probability distribution defined in the transformed space, NCA labels each query instance by the majority vote of all training samples. As a main advantage of NCA, a continuous and differentiable cost function in respect of the linear transformation matrix A is deduced. This cost function can be minimized by gradient descent. The learnt matrix A can also be forced to be low-rank, thus accelerating K-NN test and facilitating class structure visualization. Though the cost used in NCA is differentiable, it seems to be sensitive to the initialization. Inspired by NCA, Weinberger et al. proposed a Large Margin Nearest Neighbor (LMNN) method to learn a Mahalanobis distance metric for K-NN classification [START_REF] Weinberger | Distance metric learning for large margin nearest neighbor classification[END_REF]. LMNN attempts to classify the K nearest neighbors as the same class label, under the constraint that different classes should be separated by a large margin. The learning problem is formulated as a semi-definite programming problem. The corresponding cost function consists of two terms; the first term penalizes large dissimilarities between instances with the same class label in a predefined neighborhood; while as a hinge loss, the second term penalizes small dissimilarities between instances with different class labels in the whole training pool. As a convex function in respect of the matrix A, the cost function of LMNN can be optimized efficiently.

Different from the global learning methods such as Linear Discriminant Analysis (LDA) and Principal Component Analysis (PCA), both NCA and LMNN can adapt to the local structure of the application at hand. By learning a local dissimilarity metric, they effectively improved the K-NN classification accuracy in many situations. However, since they were not designed specifically for tackling data that contains unreliable input features, their performance may severely decline with this kind of imperfect information.

In this paper, our goal is to maximize the accuracy and efficiency of the EK-NN classifier on data that contains unreliable input features. To this end, we propose to learn an adaptive dissimilarity metric from this kind of imperfect data in the belief function framework. By using samples in the training pool as independent items of evidence, the belief regarding the class membership of each instance is modeled and refined using DST. A specific cost function consisting of two terms is constructed for learning a low-dimensional transformation matrix A. The first term attempts to minimize the imprecision regarding the class membership of each instance. The 2,1norm regularization of A acts as the second term, considering its good property for feature selection as already shown in muti-task learning [START_REF] Evgeniou | Multi-task feature learning[END_REF], multiclass classification [START_REF] Xiang | Discriminative least squares regression for multiclass classification and feature selection[END_REF], semisupervised learning [START_REF] Han | Semisupervised feature selection via spline regression for video semantic recognition[END_REF], etc. By means of feature selection, it aims to manage the influence of unreliable input features on the output transformation. The proposed cost function is solved efficiently by a first order method (namely the proximal forward-backward splitting algorithm [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]). The influence of the sparsity regularization is tuned according to the application at hand. Finally, a low-dimensional transformation of the input space is realized to widely separate instances from different classes, therefore increasing the classification accuracy and reducing the searching time of the EK-NN classifier simultaneously.

The rest of this paper is organized as follows. The background on DST and the EK-NN classification rule is recalled in Section II. The proposed metric learning method based on DST is then introduced in Section III. In Section IV, the proposed method is tested on both synthetic and real data sets, and some comparison with other methods is presented. Finally, we conclude paper in Section V.

II. BACKGROUND

The necessary background on DST and the EK-NN classification rule are briefly reviewed in Sections II-A and II-B, respectively.

A. Dempster-Shafer Theory

DST is also known as the theory of belief functions or evidence theory. As a generalization of both probability theory and the set-membership approaches, DST has two main components, i.e., quantification of a piece of evidence and combination of different items of evidence.

1) Evidence Quantification: DST is a formal framework for reasoning under uncertainty based on the modeling of evidence [START_REF] Shafer | A mathematical theory of evidence[END_REF]. Let ω be a variable taking values in a finite domain

Ω = {ω 1 , • • • , ω c }, called the frame of discernment.
An item of evidence regarding the actual value of ω can be represented by a mass function m on Ω, defined from the powerset

2 Ω to the interval [0, 1], such that A⊆Ω m(A) = 1. (1)
Each number m(A) denotes a degree of belief attached to the hypothesis that "ω ∈ A". Function m is said to be normalized if m(∅) = 0, which is assumed in this paper. Any subset A with m(A) > 0 is called a focal element of mass function m. If all focal elements are singletons, m is said to be Bayesian; it is then equivalent to a probability distribution. A mass function m with only one focal element is said to be categorical and is equivalent to a set. Corresponding to a normalized mass function m, we can associate belief and plausibility functions from 2 Ω to [0, 1] defined as:

Bel(A) = B⊆A m(B); (2)
P l(A) = B∩A =∅ m(B). (3)
Quantity Bel(A) (also known as credibility) can be interpreted as the degree to which the evidence supports A, while P l(A) can be interpreted as the degree to which the evidence is not contradictory to A. Functions Bel and P l are linked by the relation P l(A) = 1 -Bel(A). They are in one-to-one correspondence with mass function m.

2) Evidence Combination: In DST, beliefs are elaborated by aggregating different items of evidence. The basic mechanism for evidence combination is Dempster's rule of combination [START_REF] Shafer | A mathematical theory of evidence[END_REF]. Since Dempster's rule cannot well manage high conflicts between different pieces of evidence, various alternatives to it have been developed to tackle this problem under different situations, e.g., the TBM conjunctive and disjunctive combination rules [START_REF] Smets | The combination of evidence in the transferable belief model[END_REF], Yager's rule [START_REF] Yager | On the Dempster-Shafer framework and new combination rules[END_REF], Dubois-Prade's rule [START_REF] Dubois | Representation and combination of uncertainty with belief functions and possibility measures[END_REF], the weighted average [START_REF] Murphy | Combining belief functions when evidence conflicts[END_REF], [START_REF] Liu | Combination of sources of evidence with different discounting factors based on a new dissimilarity measure[END_REF], and the cautious and bold disjunctive rules [START_REF] Denoeux | Conjunctive and disjunctive combination of belief functions induced by nondistinct bodies of evidence[END_REF] etc. In addition, the discounting strategy has been used in some other methods to deal with the conflicts, and a new dissimilarity measure consisting of both the conflict and distance has been introduced in [START_REF] Liu | Combination of sources of evidence with different discounting factors based on a new dissimilarity measure[END_REF] to determine the discounting factor of each source of evidence to be combined. The conflicts have also been used for detecting the change occurrences in the fusion of multi-temple information [START_REF] Liu | Dynamic evidential reasoning for change detection in remote sensing images[END_REF], like in the change detection of remote sensing. Nevertheless, these alternative methods usually increase the complexity for applications, Dempster's rule still remains the most popular one for combining independent evidence.

Let m 1 and m 2 be two mass functions derived from independent items of evidence. They can be fused via Dempster's rule to induce a new mass function m 1 ⊕ m 2 defined as

(m 1 ⊕ m 2 )(A) = 1 1 -Q B∩C=A m 1 (B)m 2 (C), (4)
where Q = B∩C=∅ m 1 (B)m 2 (C) measures the degree of conflict between evidence m 1 and m 2 .

B. Evidential K-NN Classification Rule

An EK-NN classifier was proposed in [START_REF] Denoeux | A K-nearest neighbor classification rule based on Dempster-Shafer theory[END_REF] based on DST. Depending on the informativeness of the training samples with respect to the class membership of the query pattern, the EK-NN classifier computes a mass function over the whole frame of classes, and provides a global treatment of imperfect training knowledge with uncertainty.

Let

{(X i , Y i )|i = 1, • • • , N } be a collection of N training pairs, in which X i = [x 1 , • • • , x V ] T is the ith training sample with V features and Y i ∈ {ω 1 , • • • , ω c }
is the corresponding class label. Given a query instance X t , its class membership can be determined through the following steps:

• Each neighbor of X t is considered as an item of evidence that supports certain hypotheses regarding the class membership of X t . Let X j be one of its K nearest neighbors with class label Y j = ω q . The mass function induced by X j , which supports the assertion that X t also belongs to

ω q is m t,j ({ω q }) = α exp(-γ q d 2 t,j ), m t,j (Ω) = 1 -α exp(-γ q d 2 t,j ), (5)
where d t,j is the distance between X j and X t , while α and γ are two tuning parameters that can be optimized via the method proposed in [START_REF] Zouhal | An evidence-theoretic K-NN rule with parameter optimization[END_REF].

• Dempster's rule (4) is then executed to combine all neighbors' knowledge and obtain a global mass function for X t . The lower and upper bounds for the belief of any specific hypothesis are then quantified via the credibility (2) and plausibility (3) values, respectively. In the case of {0,1} losses, the final decision on the class label of X t can be made alternatively through maximizing the credibility, the plausibility, or the pignistic probability, as defined by Smets [START_REF] Smets | The transferable belief model[END_REF]. As an adaptive version of the EK-NN classifier, a neural network classifier based on DST has been proposed in [START_REF] Denoeux | A neural network classifier based on Dempster-Shafer theory[END_REF]. Some other alternatives to the EK-NN method have also been developed. For instance, the credal classification methods [START_REF] Liu | A new belief-based K-nearest neighbor classification method[END_REF]- [START_REF] Liu | A new incomplete pattern classification method based on evidential reasoning[END_REF] have been proposed by Liu et al. to deal with the overlapping classes in different cases. These methods permit the objects to be associated with not only the single classes but also meta-classes (i.e., disjunction of several classes) with different masses of belief, thus endowing the ability to specify the imprecision of classification.

III. EVIDENTIAL DISSIMILARITY METRIC LEARNING

A new approach, called evidential dissimilarity metric learning (EDML), is proposed in this section. By learning an adaptive dissimilarity measure on training samples that contain unreliable input features, EDML aims to maximize the performance of the EK-NN classifier.

A. Criterion of EDML Let {(X i , Y i )|i = 1, • • • , N } be a collection of N training pairs, in which X i = [x 1 , • • • , x V ]
T is the ith observation with V input features, and Y i is the corresponding class label taking values in a frame of discernment Ω = {ω 1 , • • • , ω c }. Assume the dissimilarity between instances X i and X j is quantified by a squared distance measure:

d 2 (X i , X j ) = (X i -X j ) T A T A(X i -X j ). (6)
Then, EDML attempts to find an optimal matrix A ∈ R v×V under the constraint v V . Such a linear transformation of the input space can boost the performance of the EK-NN classifier, since important features will have a strong impact when calculating the distance; classification is also faster in the low-dimensional transformed space.

To learn such a matrix A, we regard each X i as a query instance. Then, the squared distance between X i and X j (i.e., d 2 (X i , X j )) is used in [START_REF] Yager | On the Dempster-Shafer framework and new combination rules[END_REF], so as to represent the partial knowledge concerning the class membership of X i that offered by training sample (X j , Y j = ω q ). Parameters α and γ used in [START_REF] Yager | On the Dempster-Shafer framework and new combination rules[END_REF] are restricted to be one for simplification.

Let Γ q (q = 1, . . . , c) be the set of training samples (except X i ) belonging to the same class ω q . Since the corresponding mass functions point to the same hypothesis (i.e., Y i = ω q ), they can be combined via Dempster's rule (namely (4)) to deduce a global mass function for all training samples in Γ q :

m Γq i ({ω q }) = 1 -j∈Γq [1 -exp {-d(X i , X j )}] m Γq i (Ω) = j∈Γq [1 -exp {-d(X i , X j )}] . (7)
For q = 1, . . . , c, the global mass function m

Γq i quantifies the evidence refined from the training pool that support the assertion Y i = ω q . The mass of belief m Γq i (Ω) measures the imprecision of this evidence. In other words, it can be regarded as the calculation of the unreliability of the hypothesis Y i = ω q . If the actual value of Y i is ω q , the corresponding imprecision should then close to zero, i.e., m Γq i (Ω) ≈ 0; in contrast, imprecision pertaining to other hypotheses should close to one, i.e., m Γr i (Ω) ≈ 1, for ∀r = q. According to this assumption, we propose to represent the prediction loss for training sample (X i , Y i ) as

loss i (A) = c q=1 t i,q •    1 -m Γq i ({ω q }) • c r =q m Γr i (Ω)    2 , (8)
where t i,q is the qth element of a binary vector t i = {t i,1 , . . . , t i,c }, with t i,q = 1 if and only if

Y i = ω q . When Y i = ω q is true, minimizing loss i (A) can force both m Γq i ({ω q }) = 1 -m Γq i (Ω) and c r =q m Γr i (Ω)
to approach one as far as possible, thus achieving the goal to maximize the reliability of the right hypothesis (Y i = ω q ) but minimize the reliability of other assertions. As the result, the learnt matrix A can lead X i only close to samples from the same class in the transformed space, thus protecting the classification performance of the EK-NN method.

Therefore, for all training samples, the loss function in respect of the transformation matrix A can be finally defined as

l(A) = 1 N N i=1 loss i (A) + λ||A|| 2,1 , (9)
where loss i (A) represents the learning cost for training sample (X i , Y i ) that quantified by [START_REF] Denoeux | Maximum likelihood estimation from uncertain data in the belief function framework[END_REF]. The 2,1 -norm sparsity regularization

||A|| 2,1 = V j=1 v i=1 A 2 i,j 1/2 (10)
is imported to select input features. By forcing columns of the transformation matrix A to be zero during the learning procedure, this sparsity term only selects the most reliable input features to calculate the linear transformation, thus controlling the influence of unreliable input features on the output lowdimensional transformed space. Scalar λ is a hyper-parameter that controls the influence of this regularization. Generally speaking, a too small λ may fail to limit the influence of unreliable input features; while, a too large λ may also delete informative features.

B. Optimization

Since loss i (8) is differentiable in respect of matrix A while ||A|| 2,1 [START_REF] Su | Parametric regression analysis of imprecise and uncertain data in the fuzzy belief function framework[END_REF] is partly smooth (it is non-smooth when and only when A = 0), the proximal Forward-Backward splitting (FBS) algorithms [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF], [START_REF] Raguet | A generalized forward-backward splitting[END_REF], which belong to the class of first order methods, are efficient alternatives to solve the proposed loss function [START_REF] Ramasso | Making use of partial knowledge about hidden states in HMMs: an approach based on belief functions[END_REF]. More specifically, as an improved version of the classical FBS methods, the Beck-Teboulle proximal gradient algorithm [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] is used in this paper considering its computational simplicity and fast convergence rate.

In general, each iteration of the FBS algorithms can be broken up into a gradient descent step using 1 N N i=1 loss i (A), followed by a proximal operation using ||A|| 2,1 . According to ( 6)-( 8), the derivative of loss i concerning A (i.e., ∂loss i /∂A) can be deduced as

∂loss i ∂A = c q=1 2t i,q 1 -m Γq i ({ω q }) c r =q m Γr i (Ω) - ∂m Γq i ({ω q }) ∂A c r =q m Γr i (Ω) -m Γq i ({ω q }) c r =q ∂m Γr i (Ω) ∂A c s =r,q m Γs i (Ω) . (11)
In which, value m Γq i is calculated via [START_REF] Smets | The transferable belief model[END_REF], and for ∀q = 1, . . . , c, ∂m

Γq i ({ω q }) ∂A = - j∈Γq ∂m ij (Ω) ∂A l∈Γq\j m il (Ω); (12)
Algorithm 1: Beck-Teboulle proximal gradient algorithm [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] Initialize A (0) ∈ R v×V and β > 0, set H (0) = A (0) and t (0) = 1 ; for n = 0, 1, 2, . . . do

G (n) = H (n) -β -1 N N i=1 ∂lossi ∂A | A=H (n) ; A (n+1) = prox β -1 ,||A||2,1 G (n) = arg min A * λ||A * || 2,1 + β 2 ||A * -G (n) || 2 ; t (n+1) = [ 4t (n) 2 + 1 + 1]/2 ; δ (n) = 1 + [t (n) -1]/t (n+1) ; H (n+1) = A (n) + δ (n) [A (n+1) -A (n) ] ; end
While, mass m ij is determined using ( 5) and ( 6), and

∂m ij (Ω) ∂A = 2m ij ({ω q })A(X i -X j )(X i -X j ) T . ( 13
)
Based on ( 11)-( 13), the Beck-Teboulle proximal gradient algorithm executes as the form shown in Algorithm 1, so as to deduce an optimal or at least sub-optimal low-dimensional transformation matrix A. To facilitate the optimization procedure, classical metric learning methods (e.g., PCA) can also be used to generate the initialization (i.e. A (0) ) for the proposed method. The learnt matrix A is then applied in [START_REF] Dubois | Representation and combination of uncertainty with belief functions and possibility measures[END_REF] to measure the dissimilarity between different instances, and finally used in the EK-NN classification.

Complexity of the gradient calculation: Given m ij and m Γq i , for i = 1, . . . , N , the calculation of ∂loss i /∂A can be performed in the following order:

1) calculation of ∂m ij (Ω)/∂A for j = 1, . . . , N -1 using (13). Since (X i -X j )(X i -X j ) T can be determined beforehand, this step requires (N -1)vV 2 arithmetic operations; 2) calculation of ∂m Γq i /∂A for q = 1, . . . , c using [START_REF] Masson | ECM: An evidential version of the fuzzy C-means algorithm[END_REF], which requires c q=1 |Γ q |vV = (N -1)vV arithmetic operations; 3) calculation of the last term in [START_REF] Denoeux | EVCLUS: evidential clustering of proximity data[END_REF], which requires cvV operations. Since t i has and only has one nonzero element, the calculation of (11) needs cvV operations. Based on above steps, the complexity for calculating

∂loss i /∂A is O vV 2 N + vV (N + c) .

IV. EXPERIMENTAL RESULTS

The presented experiments consist of four parts. In the first part, the proposed method, namely EDML, was evaluated on a synthetic data set. The proportion of unreliable (noisy and imprecise) features in this synthetic data was varied to assess the robustness of EDML under different situations. In the second part, EDML was evaluated on several real data sets. The corresponding classification accuracy was compared with some other metric learning methods. The parameters used in the proposed method were also studied. Finally, we further compared the two-dimensional visualization performance of different metric learning methods, so as to evaluate whether the proposed method can effectively separate instances from different classes in low-dimensional subspaces.

A. Performance on Synthetic Data

The studied synthetic data sets were generated using a process similar to the one described in [START_REF] Weston | Use of the zero norm with linear models and kernel methods[END_REF]. The feature space contains n r relevant features uniformly and independently distributed between -1 and +1. To obtain a high non-linear discriminant surface, the output label for a given instance is defined as

y = ω 1 if max i (x i ) > 2 1-1 nr -1, ω 2 otherwise, (14)
where x i is the ith relevant feature. Besides the relevant features, there are n u irrelevant (noisy) features also uniformly distributed between -1 and +1, without any relation with the class label; and also n i imprecise features copied as the cubic of the relevant features. The numbers of relevant, irrelevant and imprecise features were set, respectively, as n r = 2, n u ∈ {6, 16, 26, 36, 46} and n i = 2 to simulate five different situations. Under each situation, we generated 150 training instances and the same number of testing instances. PCA, NCA, LMNN and the proposed EDML methods were executed to learn a twodimensional dissimilarity metric A (i.e., ∈ R nr×(nr+nu+ni) ) on the training set. The obtained metric A was then used in the EK-NN to classify both the training and testing samples. As two baselines, results obtained by the SVM and Evidential Neural Network (ENN) [START_REF] Denoeux | A neural network classifier based on Dempster-Shafer theory[END_REF] classifiers joint with PCA were also included for comparison.

Parameters of each method used in this experiment (four metric learning methods, i.e., PCA, NCA, LMNN and EDML, and three classifiers, i.e., EK-NN, ENN and SVM) can be summarized as follows:

• For LMNN, as suggested by [START_REF] Weinberger | Distance metric learning for large margin nearest neighbor classification[END_REF], parameters K and µ were set as K = 3 and µ = 0.5.

• For the proposed EDML, a rough grid search strategy was used to select an appropriate λ from {0.005, 0.007, 0.009} according to the training performance.

• For the EK-NN classifier, parameters α and γ were optimized via the operation proposed in [START_REF] Zouhal | An evidence-theoretic K-NN rule with parameter optimization[END_REF]. The number of nearest neighbors was set as K = 3.

• For the SVM, the gaussian kernel was used with the radial basis σ = 1.

• For the ENN classifier, the number of prototypes per class was set as 5. It is worth illustrating that the parameters of the compared methods were always kept the same in the sequel experiments.

Finally, the training and testing (more important) accuracy (in %) obtained by different metric learning methods are summarized in Table I, in which EDML (manually set λ = 0) and EDML-FS (namely EDML joint with Feature Selection) represent, respectively, the proposed method without/with the sparsity regularization. As can be seen, the proposed EDML-FS led to higher testing accuracy than other methods under all the five different situations. It is also worth noting that the difference increased following the augment of unreliable input features, which reveals that the proposed method is stable and immune to severely deteriorated input information.

B. Performance on Real Data

The proposed method was further evaluated using five real data sets of varying input features and classes. Four of these data sets (Wine, Seeds, Soybean-small and LSVT voice rehabilitation [START_REF] Tsanas | Objective automatic assessment of rehabilitative speech treatment in parkinson's disease[END_REF] data) were downloaded from the UCI Machine Learning Repository1 . The other one is the Olivetti face recognition data set2 . As a preprocessing operation for the Face data, we down-sampled the images to 38 × 31 pixels and used PCA to further reduce the dimensionality to 100. Properties of all the five data sets are briefly summarized in Table II.

The training and testing instances were randomly generated with 70/30 splitting, and repeated 50 times. Under each random split, we used PCA, NCA, LMNN and the proposed EDML methods, respectively, to learn a low-dimensional dissimilarity metric A (i.e., ∈ R v×V with v ≤ V ) on the training data; then used it in the EK-NN to classify both the training and testing instances. Parameters of compared methods were the same as that used in the last experiment (namely Section IV-A). For the proposed method, the hyperparameter λ was still determined by a rough grid search strategy according to the training performance. On average, good results were obtained with λ between [0.0005, 0.01] for the five data studied in this experiment.

The value of the output dimension v was orderly set as {2, 3, . . . , 15}. Then, the average testing accuracy with different v was calculated and is shown in Fig. 1. As can be seen, the proposed method consistently performed well on these data sets as compared with other methods. More specifically, LMNN (blue line) and EDML (magenta line) had comparable testing accuracy on Wine and Faces data sets; NCA (green line), LMNN and EDML resulted in almost the same performance (EDML was slightly better) on the soybeansmall data set; and EDML yielded the best performance on the other two data (Seeds and LSVT).

To further analyze the experimental results obtained on these real data sets, we computed the average training performance as a criterion to select the best output dimension v (from {2, 3, . . . , 15}) for the learnt dissimilarity metric A. The best training accuracy and the corresponding testing accuracy (more important) for each method are summarized in Table III, in which results obtained by the SVM and ENN joint with PCA are also presented as two baselines for comparison. As in the former subsection, EDML (manually set λ = 0) and EDML-FS represented the proposed method without/with the sparsity regularization. From Table III, it can be found that EDML-FS consistently yielded better performance than other methods on the first four data sets, especially on the LSVT data. This is mainly because the proposed method only selected the most informative features (from all the three hundred input features) to calculate the linear transformation. LMNN outperformed our method on the Face data set, but only with a slight difference. In addition, we can also see that, thanks to the 2,1 -norm sparsity regularization, EDML-FS performed better than EDML.

C. Parameter Analysis 1) Output Dimension: As discussed above, the best output dimension v (from {2, 3, . . . , 15}) for the five real data sets studied in the last subsection was determined according to the training performance. Therefore, besides the classification accuracy presented in Table III, the corresponding output dimension obtained by different methods on these real data sets was also summarized and is shown in Fig 2 .  2) Regularization Parameter: the hyper-parameter λ in the loss function [START_REF] Su | Parametric regression analysis of imprecise and uncertain data in the fuzzy belief function framework[END_REF] controls the effect of the sparsity regularization on the output low-dimensional transformation. It should be tuned specifically for each data set at hand. Generally speaking, a too small λ may fails to limit the influence of unreliable input features, while a too large λ may also removes many significant input features. On average, good results were obtained with λ between [0.0005, 0.01] for all the five data studied in the last subsection. In this experiment, the LSVT data set was used as an example to further analysis the influence of the parameter λ. The training and testing data were generated with 70/30 splitting, and repeated 20 times. We orderly selected a λ from {0, 0.0005, 0.0006, . . . , 0.001, 0.002, . . . , 0.015} to learn a low-dimensional transformation (with v = 5) of the input space. Then, the EK-NN classifier (with K = 3) was used to classify the testing instances on the transformed space. For all the 20 random splits, the average testing accuracy (in %) with regard to λ is finally shown in Fig. 3, in which the horizontal line represents the average accuracy of the EK-NN classification in the input space. As can be seen, relatively high performance on this data set is obtained with λ between [0.0005, 0.01]. The classification is less sensitive in the region [0.0005, 0.005] than in other regions of λ.

3) Number of Nearest Neighbors: we also studied the parameter K of the EK-NN classification with the dissimilarity metric learnt by the proposed method. Still on the LSVT data set, the training and testing data were generated with 70/30 splitting, and repeated 20 times. Under each random split, we used the proposed method to learn a low-dimensional transformation of the input space. The output dimension and the regularization parameter were set as v = 5 and λ = 0.002. Then, the EK-NN classifier with K = {1, 2, . . . , 30} was orderly executed to classify the testing instances in the transformed space. As for comparison, the EK-NN classifier with the same K was also directly executed in the input space to classify the testing instances. The average testing accuracy with regard to K is finally summarized in Fig. 4. It can be found that, with a metric learnt by the proposed method, the EK-NN classification always has higher accuracy on this data set than directly using the Euclidian distance in the input space. In addition, we can also see that the proposed method is robust to the parameter K.

(n) (o) (p) (q) (r) (s) (t) (u) (v) (w) (i) (j) (k) (l) (m)

D. Two-Dimensional Visualization

To further evaluate whether the proposed method can effectively separate instances from different classes in lowdimensional transformation space, we visualized the dimension reduction in 2D, as shown in Fig. 5. PCA, NCA, LMNN and the proposed method were still compared on the same data sets (one synthetic and five real data sets) used in the former subsections. The input feature space for the synthetic data was set as fifty (n r = 2,n i = 2 and n u = 46). For simplicity, only the first seven classes were studied in the Faces data. From the obtained results we can see that instances from different classes were always well separated by our method on all the six data sets. It led to the largest margin on the synthetic data, and the most satisfying separation on the Seeds data. In contrast, NCA did not separate the LSVT data perfectly; while LMNN resulted in large overlaps on the synthetic data.

V. CONCLUSION

To optimize the performance of the EK-NN classification on imperfect data sets, an approach based on Dempster-Shafer theory has been proposed to learn a dissimilarity metric specifying for the application at hand. By treating other training patterns as different sources of information, the belief concerning the class membership of each query pattern has been quantified and refined in the belief function framework. A specific loss function consisting of two terms has been developed for metric learning under uncertainty, in which the first term is used to minimize the imprecision regarding each instance's class membership, while the second term is the 2,1 -norm sparsity regularization of the low-dimensional transformation matrix. Through a feature selection procedure, it serves to limit the influence of uncertainty and/or imprecise input features. The proposed method has been evaluated on several synthetic and real data sets, consistently showing good performance with regard to classification accuracy and class structure visualization. Moreover, it has also proved that the proposed method is not sensitive to the parameter K.

Future work will focus on two main aspects. First, to improve the efficiency of the proposed method on large data sets, we will further study other more advanced optimization algorithms. Furthermore, it is worth extending the proposed method to learn non-linear transformation of the input space, so as to improve the performance of the EK-NN classifier on high complex data sets.

Fig. 1 .

 1 Fig. 1. Average testing accuracy obtained by different metric learning methods: (a) Wine data, (b) Seeds data, (c) Soybean-small data, (d) LSVT data and (e) Faces data. In each subfigure, the horizontal axis represents the output dimension (i.e., v) of the learnt transformation A, while the vertical axis represents the corresponding classification accuracy (in %).

Fig. 3 .

 3 Fig. 3. Average testing accuracy on the LSVT data set with regard to the hyper-parameter λ. The output dimension was set as v = 5. The horizontal line represents the accuracy obtained in the input space.

Fig. 4 .

 4 Fig. 4. Average accuracy of the EK-NN classification on the LSVT data set with regard to the number of nearest neighbors K. The output dimension was set as v = 5.

Fig. 5 .

 5 Fig. 5. Two-dimensional transformation results obtained by PCA, NCA, LMNN and the proposed method (orderly from the first to the forth column). (a)-(d) on synthetic data; (e)-(h) on Wine data; (i)-(l) on Seeds data; (m)-(p) on Soybean data; (q)-(t) on LSVT data; (u)-(x) on Faces data;

TABLE I CLASSIFICATION

 I ACCURACY (BOTH TRAINING AND TESTING, IN %) OF THE EK-NN BASED ON DIFFERENT METRIC LEARNING METHODS. IN THE STUDIED SYNTHETIC DATA SETS, nr = 2 AND n i = 2. EDML-FS AND EDML DENOTE, RESPECTIVELY, THE PROPOSED METHOD WITH/WITHOUT THE 2,1 -NORM SPARSITY REGULARIZATION. PERFORMANCE OF THE SVM AND ENN CLASSIFIERS JOINT WITH PCA WERE ALSO PRESENTED AS TWO BASELINES FOR COMPARISON.

	training	nu 6 16 26 36 46	SVM 92.00 83.33 81.33 77.33 76.67	ENN 91.33 85.33 83.33 76.00 74.67	PCA 90.67 84.67 74.00 76.00 68.67	NCA 99.33 100.00 100.00 100.00 100.00	LMNN 98.00 96.00 96.67 100.00 99.33	EDML 98.67 99.33 100.00 99.33 99.33	EDML-FS 99.33 100.00 100.00 100.00 100.00
	testing	6 16 26 36 46	85.33 74.00 66.67 69.33 64.67	84.00 72.00 69.33 69.67 66.00	84.00 73.33 64.67 62.00 57.33	91.33 84.00 78.67 70.00 82.67	90.00 86.67 78.67 78.00 76.67	86.00 86.00 84.67 76.67 76.67	94.67 92.00 90.00 95.33 94.00

TABLE II PROPERTIES

 II OF THE FIVE REAL DATA SETS STUDIED IN SECTION IV-B.

	data sets classes input features instances Wine 3 13 178 Seeds 3 7 210 Soybean 4 35 47 LSVT 2 309 126 Faces 40 100 400

TABLE III THE

 III BEST TRAINING AND THE CORRESPONDING TESTING ACCURACY (AVE±STD, IN %) OBTAINED BY DIFFERENT METHODS WITH v ∈ {2, 3, . . . , 15}. EDML-FS AND EDML DENOTE, RESPECTIVELY, THE PROPOSED METHOD WITH/WITHOUT THE SPARSITY REGULARIZATION. RESULTS OF THE SVM ANDENN JOINT WITH PCA WERE ALSO PRESENTED AS TWO BASELINES FOR COMPARISON.

	training	SVM 99.95±0.19 100.00±0.00 97.65±1.24 100.00±0.00 100.00±0.00 100.00±0.00 100.00±0.00 ENN PCA NCA LMNN EDML EDML-FS 99.58±0.47 98.74±0.644 92.90±1.87 97.04±1.33 95.86±1.50 99.13±0.70 98.52±1.27 Soybean 100.00±0.00 100.00±0.00 98.85±1.72 100.00±0.00 100.00±0.00 100.00±0.00 100.00±0.00 Wine Seeds LSVT 99.52±0.83 93.82±1.28 80.71±3.94 98.63±1.69 100.00±0.00 96.45±1.46 97.14±1.33 Faces 83.43±2.69 99.99±0.05 90.28±1.34 99.87±0.24 100.00±0.00 99.37±0.37 99.71±0.26
	testing	Wine Seeds Soybean 100.00±0.00 98.00±3.83 98.28±3.97 99.86±1.00 99.44±2.80 100.00±0.00 100.00±0.00 96.15±2.50 97.51±2.30 95.90±2.63 96.97±2.23 97.88±1.68 96.46±2.95 97.98±1.86 91.78±2.93 92.92±2.96 92.25±3.14 93.85±2.88 94.54±2.78 94.48±2.45 95.49±1.84 LSVT 80.42±4.82 85.21±3.94 80.92±5.76 82.57±6.35 82.13±5.40 85.03±4.98 86.09±4.63 Faces 65.88±4.05 85.23±3.08 89.48±2.20 89.18±3.50 97.63±1.66 93.40±2.13 97.08±1.55

Please see at https://archive.ics.uci.edu/ml/index.html

Please see at http://www.uk.research.att.com/facedatabase.html.

ACKNOWLEDGMENT

The authors thank all the anonymous referees for their helpful and constructive comments. This work was partly supported by the China Scholarship Council (CSC).


Aree Wiboonpongse

Jianxu Liu

Songsak Sriboonchitta
email: songsakecon@gmail.com

Thierry Denoeux

Modeling dependence between error components of the stochastic frontier model using copula. Application to Intercrop Coffee Production in Northern Thailand

Keywords: Stochastic Frontier, Copula, Technical Efficiency, Thailand

Modeling dependence between error components of the stochastic frontier model using copula: Application to intercrop coffee production in Northern Thailand

Aree Wiboonpongse, Jianxu Liu, Songsak Sriboonchitta, Thierry Denoeux

Introduction

Stochastic frontier models (SFMs) have been proved very useful to measure technical efficiency of production units. The stochastic frontier production model for a cross-section of observations was independently proposed by [START_REF] Aigner | Formulation and estimation of stochastic frontier production function models[END_REF] and [START_REF] Meeusen | Efficiency estimation from Cobb-Douglas production functions with composed error[END_REF]. This is essentially a linear regression model with two independent error components: a two-sided one that captures random variation of the production frontier across firms and a one-sided one that measures inefficiency relative to the frontier. In recent decades, most studies about production, cost or profit efficiency have used the conventional stochastic frontier model (see, e.g., [START_REF] Rahman | Total factor productivity growth and convergence in Northern Thai agriculture[END_REF][START_REF] Rahman | Production Efficiency of Jasmine Rice Producers Northern and North-eastern Thailand[END_REF][START_REF] Rahman | Impact of environmental production conditions on productivity and efficiency: A case study of wheat farmers in Bangladesh[END_REF], [START_REF] Sriboonchitta | Evaluation of Cost Efficiency of Thai Public Universities[END_REF], Wiboonpongse et al. (2012a[START_REF] Wiboonpongse | Inefficiency in Agricultural Markets: Asymmetric Rice Price Transmission in Thailand[END_REF]Wiboonpongse et al. ( , 2008a[START_REF] Wiboonpongse | Agricultural Production Growth Assessment for Agroeconomic Zones in Northern Thailand Using Statistical Data[END_REF][START_REF] Wiboonpongse | The Effects of Production Input, Technical Inefficiency and Biological Risk on Jasmine and Non-Jasmine Rice Yields in Thailand[END_REF], [START_REF] Wiboonpongse | The Effects of Production Inputs, Technical Efficiency and Other Factors on Jasmine and Non-jasmine Rice Yield in Production Year[END_REF], [START_REF] Sriboonchitta | Technical Efficiency of Rural Micro and Community Enterprises in the Upper North of Thailand[END_REF], 2005, 2004), [START_REF] Aigner | Formulation and estimation of stochastic frontier production function models[END_REF], [START_REF] Jondrow | On the estimation of technical inefficiency in the stochastic frontier production function model[END_REF], [START_REF] Coelli | An Introduction to Efficiency and Productivity Analysis[END_REF][START_REF] Coelli | Technical, allocative, cost and scale efficiencies in Bangladesh rice cultivation: A non-parametric approach[END_REF] Kombhakar and Lovel (2000)). In all these studies, it is assumed that the one-sided and two-sided error terms are independent, so that the joint distribution of these two random variables can easily be computed. Based on this assumption, the parameters of the SFM can be estimated using the Corrected Ordinary Least Squares (COLS) or Maximum Likelihood (ML) methods. However, the impact of the independence assumption on technical efficiency estimation has long remained an open issue. Fortunately, this strong assumption can be relaxed by using copula to fit the joint distribution of the two random error components more appropriately. [START_REF] Smith | Stochastic frontier models with dependent error components[END_REF] first proposed a SFM allowing for dependence between the two error components using copula functions. Copula functions can be used to capture rank correlation and tail dependence between the two error components, thus making the stochastic frontier analysis much more flexible. However, the log likelihood function in the copula-based stochastic frontier model generally does not have a closed form, which makes its maximization numerically complex. In this paper, we proposed to use the maximum simulated likelihood method, which has numerical and computational advantages over numerical integration method used by [START_REF] Smith | Stochastic frontier models with dependent error components[END_REF]. Furthermore, to explore the dependence structure of the error components in the SFM, we consider several copula families including the Student-t, Clayton, Gumbel and Joe families as well as their relevant rotated versions. The model with the best fit/complexity trade-off is selected using the AIC and BIC criteria. This approach was applied to cross-sectional data of coffee production in Thailand. A comparison between technical efficiencies computed with and without the independence assumption (considering the best copula model) reveals that the standard approach grossly overestimates efficiency, which has important implication for production analysis using the SFM.

The remainder of this paper is organized as follows. Section 2 introduces the necessary background on the SFM and copula. Section 3 presents the copula-based stochastic frontier approach. Empirical results with this model applied to coffee production data are reported in Section 4. Finally, Section 5 concludes the paper.

Background and theory

The SFM is a regression-like model with a disturbance function that is asymmetric and distinctly non-normal. This model will first be briefly in Section 2.1. Some background on copula will then be recalled in Section 2.2. These are the two building blocks of the copula-based model introduced in Section 3.

Stochastic frontier model

The SFM aims at representing the "ideal" relation between inputs x j and output y of a production process by a production function, seen as a theoretical ideal. The observed deviations from the production function could arise from two sources:

1. Productive inefficiency, resulting in a non-negative error term; 2. Firm-specific effects, which can enter the model with either signs. Assuming the relation between the logarithm of production log(y) and the input vector x to be linear, we obtain the following equation:

log ' i i i i i i yx vw     (1)
for i=1, 2,…,N, where N is the number of firms, y i is the production output for firm i, x i is a k×1 vector of inputs for firm i (including a constant component x i1 =1),  is a k×1 vector of unknown parameters, v i is a random error, and w i is a non-negative unobservable random variable representing production inefficiency. The two error components w i and v i are assumed to be independent. The technical efficiency TE i for firm i can then be defined as the ratio between actual to theoretically ideal production:

exp( ) exp( ). exp( )

ii ii i x v w TE w xv          (2)
It is comprised between 0 and 1. The inefficiency error term w is usually assumed to have a gamma, exponential or half-normal distribution (defined as the distribution of the absolute value of a normal variable) [START_REF] Cullinane | Estimating the Relative Efficiency of European Container Ports: A Stochastic Frontier Analysis[END_REF]. In contrast, the symmetric error term v is usually assumed to have a normal or logistic distribution.

Copula

A copula connects a given number of one-dimensional marginal distributions to form a joint multivariate distribution. In the following, we will restrict the presentation to bivariate copula, which will be used later. Sklar's theorem [START_REF] Sklar | Fonctions de repartition "an dimensions et leurs marges[END_REF] states that any bivariate cumulative distribution function (cdf) F(x 1 ,x 2 ) can be expressed as

1 2 1 1 2 2 ( , ) ( ( ), ( )), F x x C F x F x   (3)
where F 1 (.) and F 2 (.) and the marginal cdfs of X 1 and X 2 , , C  is a bivariate copula function and  is a parameter vector of the copula, commonly referred to as the dependence parameter vector. Bivariate copulas satisfy the following properties 1. C(u 1 ,0)=C(0,u 2 ) for all u 1 and u 2 in [0,1]; 2. C(u 1 ,1)=u 1 and C(1,u 2 )=u 2 for all u 1 and u 2 in [0,1]; 3. For all 0≤u 1 ≤u 2 ≤1 and 0≤v 1 ≤v 2 ≤1,

2 2 2 1 1 2 1 1 ( , ) ( , ) ( , ) ( , ) 0. C u v C u v C u v C u v    
If the random vector (X 1 ,X 2 ) has a joint density f(x 1 ,x 2 ), it can be expressed as a function of the copula density by [START_REF] Lee | Copula-based multivariate GARCH model with uncorrelated dependent errors[END_REF]:

22 1 2 1 2 1 2 1 2 1 2 1 1 2 2 1 2 1 2 1 2 ( , ) ( , ) ( ) ( ) ( , ) ( , ) ( ) ( ), F x x C u u F x F x f x x c u u f x f x x x u u x x              (4)
where f 1 (x 1 ) and f 2 (x 2 ) are the marginal densities, the small letter c represents the density function of the copula. The most common measure of dependence between random variables is Pearson"s correlation coefficient, which measures only linear dependence and is not very informative for asymmetric distributions [START_REF] Boyer | Pitfalls in tests for changes in correlations[END_REF]. To measure nonlinear dependence, rank correlation coefficients such as Kendall"s tau and Spearman"s rho are more suitable. They can be expressed in terms of the copula as follows (see Nelsen, 2006):

(6) Moreover, there is the possibility that the tail dependence may divide into upper tail and lower tail dependences, the definitions of which are written as [START_REF] Boyer | Pitfalls in tests for changes in correlations[END_REF] and

11 2 2 1 1 11 1 2 ( , ) : lim [ ( ) | ( )] lim , 1 U uu u C u u P X F u X F u u           (8)
where F 1 and F 2 are the marginal cumulative distribution functions of X 1 and X 2 , respectively. Some copulas, such as the Gaussian and Frank copulas, possess the characteristic 0 UL   , while most copulas can capture upper tail and/or lower tail dependence. For instance, Clayton copula can measure lower tail dependence, and Gumbel copula can measure upper tail dependence, and Joe copula families, which can represent upper tail dependence as well, while student-t copula reflects symmetric tail dependence [START_REF] Brechmann | Modeling Dependence with C-and D-Vine Copulas: The R Package CDVine[END_REF]. However, these copulas can only capture positive dependence except Gaussian, T and Frank copulas. Fortunately, these copulas may then be "rotated" and applied again. A thorough review of rotated copulas may be found in Christian [START_REF] Cech | Copula-based top-down approaches in financial risk aggregation[END_REF], Brechmann and Schepsmeier, 2013 and Luo (2010). There are three rotated forms, rotated 90 degrees, rotated 180 degrees and rotated 270 degrees. When rotating them by 180 degrees, one obtains the corresponding survival copulas. Non-rotated Clayton copula can capture lower tail dependence, while survival Clayton copula can measure upper tail dependence. Also, rotation by 90 and 270 degrees allows for the modeling of negative dependence. The distribution functions of rotated copulas by 90, 180 and 270 degrees are given as follows:

90 1 2 2 1 2 ( , ) (1 , ), C u u u C u u    (9) 180 1 2 1 2 1 2 ( , ) 1 (1 ,1 ), C u u u u C u u       ( 10
) 270 1 2 1 1 2 ( , ) ( ,1 ). C u u u C u u    . ( 11
) 22 12 ( , ) 1 2 1 2 1 2 1 2 1 2 [0,1] [0,1] 12 ( , ) 3 12 ( , ) 3 12 ( ) 3
.

XX S u u dC u u C u u du du E U U          2 12 , 1 2 1 2 1 2 [0,1] 4 ( , ) ( , ) 1 4 [ ( , )] 1, XX C u u dC u u E C U U       11 2 2 1 1 00 ( , ) : lim [ ( ) | ( )] lim L uu C u u P X F u X F u u        
Fré chet-Hoeffding bounds for copula functions are crucial to select an appropriate copula that should cover the sample space between the lower and the upper bounds and that as a parameter of copula approaches the lower (upper) bound of its permissible range. Fré chet-Hoeffding bounds for copula functions can be expressed as

1 2 1 2 1 2 1 2 1 2 ( , ) max( 1,0) ( , ) min( , ) ( , ) W u u u u C u u u u M u u       . ( 12
)
The copula approaches the Fré chet-Hoeffding lower (upper) bound. Therefore, we can make use of some copula families that cover a large dependence as much as possible, such as the (rotated) Clayton, (rotated) Gumbel and (rotated) Joe copula families, which can represent both negative and positive dependence. The copula families used in this paper are summarized below. We refer to the formulas and notations from Trivedi and Trivedi (2005) and Sriboonchitta et al. ( 2013).

(1) Independence copula

The independence copula is defined by

12 IN C u u 
. It connects the cdfs of two independent random variables.

(2) Gaussian copula This copula takes the form

      11 12 22 1 1 2 2 Ga 1 2 1 2 2 2 2 1 ( , ; ) exp d d , 2 1 21 uu x x x x C u u x x                  (13)
where 11     is Pearson"s correlation coefficient and  is the cdf the standard normal distribution function. This is the copula pertaining to a bivariate normal distribution with standard normal marginals and Pearson"s linear correlation coefficient  . Parameter  is related to the Kendall"s tau  and Spearman"s rho S  coefficients by the following equations:

1 (2 / )sin ( )      and 1 (6 / )sin ( / 2). S     . (3) T copula
As Gaussian copulas, T copulas belong to the class of elliptical copula (i.e., they are the copulas of elliptically contoured distributions). However, T copulas can capture tail dependence, and it is symmetric extreme dependence. They are defined as

11 12 22 ( ) ( ) ( 2)/2 12 2 2 12 ( , | , ) {1 } , (1 ) 21 vv T u T u v x xy y C u u v dx dy v              ( 14
)
where  is the Pearson"s correlation coefficient and T v is the cdf of a Student distribution with v degrees of freedom. When v tends to infinity, the T copula converges to the Gaussian copula. The symmetric tail dependence can be calculated as a function of parameters  and v as follows:

2 ( 1 (1 ) / (1 )). U L v Tv            (15)
(4) Clayton copula

The copula

1 1 2 1 2 ( , ) ( ( ) ( )) C u u u u     
is called Archimedean and the function  its generator. There are different Archimedean copulas according to different generators, such as Clayton, Gumbel, Frank, Joe and AMH copulas etc. Clayton copula family can reflect the lower tail dependence for 0   . It is characterized by the following formula:

1/ Cl 1 2 1 2 ( , | )=( 1) . C u u u u       . ( 16
)
The above copula can only capture a strong lower tail and positive dependence, but it can be rotated and used in capture negative dependence or reflect strong upper tail dependence (see [START_REF] Cech | Copula-based top-down approaches in financial risk aggregation[END_REF] and Luo (2010)). The corresponding Kendall"s tau measure is simply given by = / ( 2) CL     . However, the relation between the copula parameter and the Spearman"s rho is very complicated, and the lower tail dependence can be also simply calculated by

1/ 2 L     .
(5) Frank copula Copulas in the Frank family are the only Archimedean copulas that attain both lower and upper bounds, thus allowing for positive and negative dependence. The corresponding copula function is given by:

   12 Fr 1 2 exp( ) 1 exp( ) 1 1 ( , ; ) ln , 1+ exp( ) 1 uu C u u              (17)
where θ ∈(-∞, +∞)\{0}; positive (resp., negative) values of θ correspond to positive (resp., negative) dependence. The independence copula is recovered in the limit when θ →0. The rank correlation coefficients are given by:

1 () 4 1-+4 D     and 12 12 1 [ ( ) ( )], S D D       where () k D  is Debye function 0 ( ) , 1, 2. 1 k k kt kt D dt k e       (6) Gumbel copula
The bivariate Gumbel copula is given by

    1/ 1/ Gum 1 2 12 ( , | ) exp , ( ln ) ( ln ) C u u uu         (18)
where (1, )    . It is an asymmetric copula of the Archimedean family, which allows for strong upper tail dependence. The rotated Gumbel copula can be applied to capture negative dependence as well. The Kendall"s tau of Gumbel copula is given by The Joe copula is defined as follows:

1/ 1 2 1 2 1 2 ( , | ) 1 [(1 ) (1 ) (1 ) (1 ) ] , Joe C u u u u u u               (
The relationship between Spearman"s rho and parameter  does not have a closed form expression. The rotated Joe copula can describe negative dependence as well.

To summarize, different copulas have different characteristics, such as upper tail dependence, lower dependence, positive and negative dependences, etc. Therefore, the above copula families and relevant rotated copula can potentially capture the appropriate dependence between two random variables. Other popular copula families, such as the Farlie-Gumbel-Morgenstern (FGM) copula and Ali-Mikhail-Haq (AMH) copula have been discarded because they cannot achieve Fré chet bounds and can only accommodate relatively weak dependence between the margins. The ranges of dependence of Kendall"s tau and Spearman"s rho for the FGM copula are [-2/9, 2/9] and [-1/3, 1/3], respectively; the Kendall"s tau of the AMH copula is bounded to the interval [-0.1817, 0.3333], and the range of Spearman"s rho is [-0.2711, 0.4784] (see [START_REF] Kumar | Probability Distributions and Estimation of Ali-Mikhail-Haq Copula[END_REF][START_REF] Smith | Stochastic frontier models with dependent error components[END_REF].

Copula-based stochastic frontier model

The classical SFM makes the strong assumption that the error components v and w are independent. [START_REF] Smith | Stochastic frontier models with dependent error components[END_REF] proposed to relax this assumption and to model the dependence between v and w using copula. The classical model is recovered as a special case corresponding to the product copula. First, the density function of ( , ) f w v can be transformed into ( , ) fw  whose formula follows from [START_REF] Smith | Stochastic frontier models with dependent error components[END_REF].

( , ) ( , ) ( ) ( ) ( ( ), ( )) U V W V f w v f u u f u f w c F w F w            (21)
Then, the probability density function of  can be expressed as

( ) ( , ) R f f w dw      (22) or ( ) [ ( ) ( ( ), ( ))], W V W V f E f W c F W F W         (23)
where [] W E  denotes the expected value pertaining to the distribution of technical inefficiency W,  represents all parameters that are from marginals and copula function. We use the maximum likelihood method to estimate the copula-based stochastic frontier model. Assume that we have cross-sectional observations of n individuals or firms. The likelihood function is then given by 11 ( , , , ) ( ) ( )

nn w v i i i ii L f f y x              (24)
where y i is the output of individual i or firm i, x i is explanatory variable vector, w  and v  are from marginal distributions of W and V, respectively. Smith (2008) showed that there are very few density function of  for which the maximum likelihood estimate MLE has a closed form expression. Then, an alternate approach, maximum simulated likelihood method by [START_REF] Kao | Simulated Maximum Likelihood Estimation of The Linear Expenditure System with Binding Non-negativity Constraints[END_REF], is used to estimate the unknown parameters in copula-based stochastic frontier model (see Green (2010), [START_REF] Burns | The Simulated Maximum Likelihood Estimation of Stochastic Frontier Models with Correlated Error Components[END_REF]). The principle of maximum simulated likelihood is that computing the simulated log likelihood function will require R draws from the corresponding distribution of W for each observation. If W and V are assumed to be half-normal and normal distribution, then the density function of  can be transformed by using simulation to approximate the integration as follows:

0 22 0 22 0 0 0 0 0 0 2 0 00 0 ( )= ( ) ( ) ( ( ), ( )) 2 exp( / 2 ) ( ) ( ( ), ( )) 2 2 exp( ( ) / 2 ) ( ) ( ( ), ( )) 2 2 exp( / 2) ( ) ( ( ), 2 W V W V w V W V w ww V w W w V w w w V w W w f f w f w c F w F w dw w f w c F w F w dw w f w c F w F w d w w f w c F w                                                    00 0, 0, 0, 1 ( )) 1 ( ) ( ( ), ( )) Vw R V w ir i W w ir V w ir i r F w dw f w c F w F w R               ( 25
)
where w 0,ir is a sequence of R random draws from the standard half normal distribution. The simulated log-likelihood of the copula-based stochastic frontier approach can be expressed as follows:

0, 0, 0, 11 1 ( , , , ) log{ ( ) ( ( ), ( ))}, NR s w v V w ir i W w ir V w ir i ir L f w c F w F w R                 ( 26
)
where

i i i yx    .
In addition, following Battese and Coelli (1995), the variance terms are parameterized by replacing the relevant parameters,  ,  and 2  , can be defined by

/ wv     , 2 2 2
/ ( )

w w v      and 2 2 2 wv     .
The larger is  , the greater is the inefficiency component in the model [START_REF] Greene | Distinguishing between heterogeneity and inefficiency: stochastic frontier analysis of the World Health Organization"s panel data on national health care systems[END_REF]. The values of  and  may reveal whether inefficiency plays an important role in the composite error term and postulate the choice of the stochastic frontier approach in the present study [START_REF] Chen | Applying the stochastic frontier approach to measure hotel managerial efficiency in Taiwan[END_REF]. Note that the standard errors of these parameters can be calculated by applying the invariance property.

The most important analysis of copula-based stochastic frontier model is the technical efficiency that is specified with

0, 0, 0, 0, 1 0, 0 [exp( ) ] 1 exp( ) ( , ) () [exp( ) ( ) ( ( ), ( ))] [ ( ) ( ( ), ( ))] 1 exp( ) ( ) ( ( ), ( )) 1 ( ) ( ( R W V W V W V W V R w ir V w ir i W w ir V w ir i r V w ir i W w TE E W w f w du f E W f W c F W F W E f W c F W F W w f w c F w F w R f w c F w R                                         , 0, 1 ), ( )) R ir V w ir i r Fw     (27)
Since the maximum simulated likelihood estimation method in this study is applied, the technical efficiency can be calculated by simulated process as well.

Empirical results

The data

The data used in this study were collected by interviewing farmers in Chiang Mai province of Thailand. The data used in this study were from Kunnika Isarangkun na ayuttaya. The upland areas of Chiang Mai province in particular proved fertile areas for high quality coffee. The area has lower humidity, a shorter monsoon season, and a lower annual temperature, creating a micro-climate among the Chiang Mai mountains. Rain water from the mountains seeps into the soil, making it particularly rich in mineral content. The questionnaire was constructed to ask for details about the irrigated rice production at the farms. In particular, there was interest in the area grown, the yields obtained, the use of inputs, such as fertilizer and labor. Information was also obtained on social characteristics of the sample farmers. Data on a total of 111 sample farmers were obtained in the survey. Figure 1 describes the relationships between output and labor, output and fertilizer, labor and fertilizer. They show that there exist significant linear correlations by logarithm forms.

Estimating parameters of copula-based stochastic frontier model

The copula-based stochastic frontier model is presented in the form of linear translog production model,

0 1 2 3 45 log( ) log( ) log( ) (log( )) ^2 / 2 (log( )) ^2 / 2 log( ) log( ) i ii Y labor fertilizer labor fertilizer labor fertilizer v w             (28)
where Y i represents intercrop coffee output of the i th farmer. Consider the Log-term of Y is normal distribution, then, we assume that the v is normal distribution and the w is half-normal distribution. We analyze the copula-based stochastic frontier model by setting R equals 111 (the length of sample), which corresponds to 10 different copula families including independent copula, Gaussian copula, T copula, Frank copula, Clayton copula, Gumbel copula, Joe copula, rotated Clayton copula (180°) , rotated Gumbel copula (180°) and rotated Joe copula (180°). If there exist negative correlation between w and v, the rotated copulas (90° and 270°) should be used in this model. The log-likelihood function in [START_REF] Wiboonpongse | Agricultural Production Growth Assessment for Agroeconomic Zones in Northern Thailand Using Statistical Data[END_REF] was estimated by Nelder-Mead algorithm using starting values obtained from stochastic frontier model in R program.

Figure 2 describes the AIC and BIC are calculated from each copula-based stochastic frontier model. In terms of the values of AIC and BIC, the rotated Gumbel copula (180°) exhibits better explanatory ability than other dependence structures, while the independent copula performs the worst. Therefore, relaxing the assumption of independent between error components is effective, and is beneficial to improve the explanatory ability of the models. Moreover, the likelihood ratio (LR) test is used to compare independent copula with other copula models. The minimum LR statistics is from the independent copula and Gaussian copula, which equals to 3.294, and is significant at the 10% level. The maximum LR statistics exists between the independent copula and rotated Gumbel copula (180°), the value of which equals to 10.175 that is significant at the 1% level. Therefore, the results of LR statistics verify that the standard Normal-Half-Normal stochastic frontier (stochastic frontier based independent copula) model is inferior to the copula-based stochastic frontier models.

We now turn to the frontier estimates of the preferred model (rotated Gumbel copula-based stochastic frontier model). The estimated results are shown in Table 1. It shows that all parameters are significant at the 5% level, and the values of parameters β equal to 13.31, 1.03, -1.79, 0.09, 0.47 and -0.18 respectively. The variable labor has a positive and significant effect on yield, which implies a 1% increase in labor number allocated to coffee will increase production by 1.03%, but the impact of quantity of fertilizer per ha on coffee production shows a negative sign. This may be explained by the fact that in comparing the use of fertilizer farmers applied the high standard variance (433.5 kg/rai) with a range of 150 kg per acre to 2050 kg per rai. Using wide range of urea would show the inefficient use of fertilizer. It could be related to poor agronomic management which leads to nitrogen losses through NH3 volatilization and denitrification. Both the estimates of w  and v  are significantly different from zero at the 1% level. The estimated parameter of rotated Gumbel copula equals to 9.37 that is significantly different from zero at the 1% level as well, which confirms that serious dependence between error components exists, thereby justifying the use of copula-based stochastic frontier model. The parameter γ is 0.9135 that indicates the technical inefficient term has significant impact on intercrop coffee production. We transform the estimated parameter of rotated Gumbel copula to Kendall"s tau and Spearman"s rho that equal to 0.89 and 0.98, respectively. Plus, there exist strong lower tail dependence that equals to 0.92 between error components, which implies the higher technical efficiency it is, the smaller random disturbance it is. Given the specification of the preferred model, all farmers have a range of 0.19 to 0.83 efficiency scores with an average efficiency 0.53, while the technical efficiency of all farmers are from 0.45 to 0.99 according to independent copula-based stochastic frontier model, and the average technical efficiency is 0.74. Obviously, the traditional independence assumption of random error term and technical inefficiency overestimate the technical efficiencies in this study. Thus it can be seen that employing copula can appropriately relax the independence assumption, thus avoiding the interdependence between W and V makes influence to technical efficiency. Figure 4 illustrates distribution of technical efficiency scores for all sample farmers based independent and rotated Gumbel copulas, respectively. More than half of the farmers have high technical efficiencies that are great than 0.7 and less than 0.99 in production of intercrop coffee based on traditional SFM. However, the technical efficiencies of the preferable model show that there are a half of farmers who have high technical efficiencies but they are less than 0.9. Under the present technology and the preferable model, there is a large proportion (38%) of the sample farmers who have lower levels (<0.30) of technical efficiency score. This suggests that considerable amount of productivity is lost due to inefficiency.

Technical efficiencies

Conclusions

In this paper, we have summarized and extended the copula-based stochastic frontier model using maximum simulated likelihood estimation method that is used as substitute for numerical integral method thereby possibly generalizing the use of the copula-based SFA approach. The copula-based SFM allows the models of interest to capture the dependency between W and V. The information criteria, such as AIC and BIC, and LR test are applied to determine whether or not there exists any dependence between random error term and technical inefficiency. The most of copula families, basically, are used to combine with stochastic frontier model, such as Gaussian, T, Frank, Clayton, Gumbel, Joe, rotated Clayton (90°, 180° and 270°), rotated Gumbel (90°, 180° and 270°) and rotated Joe (90°, 180° and 270°) in this study. We used the copula-based stochastic frontier approach to analyze the intercrop coffee production associated with labor and fertilizer. The results showed that copula-based approaches had better performance than traditional stochastic frontier model which overestimate the technical efficiency in this study. This study found some significant results which showed that a 38% the coffee farmers had the production efficiencies are less than 0.3 from the maximum scale of 0.83 efficiency. This suggested that relatively considerable amount makers should seriously improve the efficiency of the intercrop coffee production.

Highlights

 A copula-based stochastic frontier method is used to examine technical  The maximum simulated likelihood method is employed to estimate copula-based SFM.  The copula-based approach allows us to capture dependency between U and V.  We added eight candidates of the Archimedean copulas.  The traditional stochastic frontier model overestimated the technical efficiency.

1

 1

Fig. 1 .

 1 Fig. 1. The relationship between explanatory variables and explained variable.

Fig. 2 . 1  1  2  3  4  5 

 2112345 Fig. 2. AIC and BIC of each copula-based stochastic frontier model.

Figure 3

 3 Figure 3 compares the preferred model (rotated Gumbel copula-based stochastic

Fig. 3 .

 3 Fig. 3. Technical efficiencies based independent and Rotated Gumbel copulas.

Fig. 4 .

 4 Fig. 4. Range of technical efficiencies and the percentage of sample farmers.

Faculty of Economics, Prince of Songkla University, Songkla, Thailand

Faculty of Economics, Chiang Mai University, Chiang Mai 50200, Thailand

UMR 7253 Heudiasyc, Université de Technologie de Compiègne, CNRS, France

Acknowledgement

The support of the Thailand Research Fund (TRF) for the Faculty of Economics and the Faculty of Agriculture in Chiang Mai University is here gratefully acknowledged. We also gratefully thank Professor Murray D. Smith and Professor Hung T. Nguyen for their valuable comments and suggestions.