Vladimir Vovk

Ivan Petej
email: ivan.petej@gmail.com

Valentina Fedorova

From Conformal to Probabilistic Prediction

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

Introduction

In essence, conformal predictors output systems of p-values: to each potential label of a test object a conformal predictor assigns the corresponding p-value, and a low p-value is interpreted as the label being unlikely. It has been argued, especially by Bayesian statisticians, that p-values are more difficult to interpret than probabilities; besides, in decision problems probabilities can be easily combined with utilities to obtain decisions that are optimal from the point of view of Bayesian decision theory. In this paper we will apply the idea of transforming p-values into probabilities (used in a completely different context in, e.g., [START_REF] Vovk | A logic of probability, with application to the foundations of statistics (with discussion)[END_REF], Sect. 9, and [START_REF] Sellke | Calibration of p-values for testing precise null hypotheses[END_REF]) to conformal prediction: the p-values produced by conformal predictors will be transformed into probabilities.

The approach of this paper is as follows. It was observed in [START_REF] Vovk | Criteria of efficiency for conformal prediction[END_REF] that some criteria of efficiency for conformal prediction (called "probabilistic criteria") encourage using the conditional probability Q(y | x) as the conformity score for an observation (x, y), Q being the data-generating distribution. In this paper we extend this observation to label-conditional predictors (Sect. 2).

Next we imagine that we are given a conformal predictor Γ that is nearly optimal with respect to a probabilistic criterion (such a conformal predictor might be an outcome of a thorough empirical study of various conformal predictors using a probabilistic criterion of efficiency). Essentially, this means that in the limit of a very large training set the p-value that Γ outputs for an observation (x, y) is a monotonic transformation of the conditional probability Q(y | x) (Theorem 1 in Sect. 3).

Finally, we transform the p-values back into conditional probabilities using the distribution of p-values in the test set (Sect. 5). Following [START_REF] Vovk | A logic of probability, with application to the foundations of statistics (with discussion)[END_REF] and [START_REF] Sellke | Calibration of p-values for testing precise null hypotheses[END_REF], we will say that at this step we calibrate the p-values into probabilities, In Sect. 6 we give an example of a realistic situation where use of the techniques developed in this paper improves on a standard approach. The perfor-mance of the probabilistic predictors considered in that section is measured using standard loss functions, logarithmic and Brier (Sect. 4).

Comparisons with related work

It should be noted that in the process of transforming p-values into probabilities suggested in this paper we lose a valuable feature of conformal prediction, its automatic validity. Our hope, however, is that the advantages of conformal prediction will translate into accurate probabilistic predictions.

There is another method of probabilistic prediction that is related to conformal prediction, Venn prediction (see, e.g., [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF], Chap. 6, or [START_REF] Vovk | Venn-Abers predictors[END_REF]). This method does have a guaranteed property of validity (perhaps the simplest being Theorem 1 in [START_REF] Vovk | Venn-Abers predictors[END_REF]); however, the price to pay is that it outputs multiprobabilistic predictions rather than sharp probabilistic predictions. There are natural ways of transforming multiprobabilistic predictions into sharp probabilistic predictions (see, e.g., [START_REF] Vovk | Venn-Abers predictors[END_REF], Sect. 4), but such transformations, again, lead to the loss of the formal property of validity.

As preparation, we study label-conditional conformal prediction. For a general discussion of conditionality in conformal prediction, see [START_REF] Vovk | Conditional validity of inductive conformal predictors[END_REF]. Objectconditional conformal prediction has been studied in [START_REF] Lei | Distribution free prediction bands for nonparametric regression[END_REF] (in the case of regression).

Criteria of efficiency for label-conditional conformal predictors and transducers

Let X be a measurable space (the object space) and Y be a finite set equipped with the discrete σ-algebra (the label space); the observation space is defined to be Z := X × Y. A conformity measure is a measurable function A that assigns to every sequence (z 1 , . . . , z n ) ∈ Z * of observations a same-length sequence (α 1 , . . . , α n ) of real numbers and that is equivariant with respect to permutations: for any n and any permutation π of {1, . . . , n},

(α 1 , . . . , α n ) = A(z 1 , . . . , z n ) =⇒ α π(1) , . . . , α π(n) = A z π(1) , . . . , z π(n) .
The label-conditional conformal predictor determined by A is defined by

Γ (z 1 , . . . , z l , x) := {y | p y > } , (1)
where (z 1 , . . . , z l ) ∈ Z * is a training sequence, x is a test object, ∈ (0, 1) is a given significance level, and for each y ∈ Y the corresponding label-conditional p-value p y is defined by

p y := i = 1, . . . , l + 1 | y i = y & α y i < α y l+1 |{i = 1, . . . , l + 1 | y i = y}| + τ i = 1, . . . , l + 1 | y i = y & α y i = α y l+1 |{i = 1, . . . , l + 1 | y i = y}| , (2)
where τ is a random number distributed uniformly on the interval [0, 1] and the corresponding sequence of conformity scores is defined by (α y 1 , . . . , α y l , α y l+1 ) := A(z 1 , . . . , z l , (x, y)).

It is clear that the system of prediction sets (1) output by a conformal predictor is nested, namely decreasing in .

The label-conditional conformal transducer determined by A outputs the system of p-values (p y | y ∈ Y) defined by [START_REF] Grenander | On the theory of mortality measurement[END_REF] for each training sequence (z 1 , . . . , z l ) of observations and each test object x.

Four criteria of efficiency

Suppose that, besides the training sequence, we are also given a test sequence, and would like to measure on it the performance of a label-conditional conformal predictor or transducer. As usual, let us define the performance on the test set to be the average performance (or, equivalently, the sum of performances) on the individual test observations. Following [START_REF] Vovk | Criteria of efficiency for conformal prediction[END_REF], we will discuss the following four criteria of efficiency for individual test observations; all the criteria will work in the same direction: the smaller the better.

-The sum y∈Y p y of the p-values; referred to as the S criterion. This is applicable to conformal transducers (i.e., the criterion is -independent). -The size |Γ | of the prediction set at a significance level ; this is the N criterion. It is applicable to conformal predictors ( -dependent). -The sum of the p-values apart from that for the true label: the OF ("observed fuzziness") criterion. -The number of false labels included in the prediction set Γ at a significance level ; this is the OE ("observed excess") criterion.

The last two criteria are simple modifications of the first two (leading to smoother and more expressive pictures). Equivalently, the S criterion can be defined as the arithmetic mean 1 |Y| y∈Y p y of the p-values; the proof of Theorem 1 below will show that, in fact, we can replace arithmetic mean by any mean ([3], Sect. 3.1), including geometric, harmonic, etc.

Optimal idealized conformity measures for a known probability distribution

In this section we consider the idealized case where the probability distribution Q generating independent observations z 1 , z 2 , . . . is known (as in [START_REF] Vovk | Criteria of efficiency for conformal prediction[END_REF]). The main result of this section, Theorem 1, is the label-conditional counterpart of Theorem 1 in [START_REF] Vovk | Criteria of efficiency for conformal prediction[END_REF]; the proof of our Theorem 1 is also modelled on the proof of Theorem 1 in [START_REF] Vovk | Criteria of efficiency for conformal prediction[END_REF]. In this section we assume, for simplicity, that the set Z is finite and that Q({z}) > 0 for all z ∈ Z.

An idealized conformity measure is a function A(z, Q) of z ∈ Z and Q ∈ P(Z) (where P(Z) is the set of all probability measures on Z). We will sometimes write the corresponding conformity scores as A(z), as Q will be clear from the context. The idealized smoothed label-conditional conformal predictor corresponding to A outputs the following prediction set Γ (x) for each object x ∈ X and each significance level ∈ (0, 1). For each potential label y ∈ Y for x define the corresponding label-conditional p-value as

p y = p(x, y) := Q({(x , y) | x ∈ X & A((x , y), Q) < A((x, y), Q)}) Q Y ({y}) + τ Q({(x , y) | x ∈ X & A((x , y), Q) = A((x, y), Q)}) Q Y ({y}) (3)
(this is the idealized analogue of ( 2)), where

Q Y is the marginal distribution of Q on Y and τ is a random number distributed uniformly on [0, 1]. The prediction set is Γ (x) := {y ∈ Y | p(x, y) > } . (4)
The idealized smoothed label-conditional conformal transducer corresponding to A outputs for each object x ∈ X the system of p-values (p y | y ∈ Y) defined by (3); in the idealized case we will usually use the alternative notation p(x, y) for p y .

Four idealized criteria of efficiency

In this subsection we will apply the four criteria of efficiency that we discussed in the previous section to the idealized case of infinite training and test sequences; since the sequences are infinite, they carry all information about the data-generating distribution Q. We will write Γ A (x) for the Γ (x) in (4) and p A (x, y) for the p(x, y) in ( 3) to indicate the dependence on the choice of the conformity measure A. Let U be the uniform probability measure on the interval [0, 1]. An idealized conformity measure A is:

-S-optimal if E (x,τ )∼Q X ×U y p A (x, y) ≤ E (x,τ )∼Q X ×U y p B (x, y) for any idealized conformity measure B, where The conditional probability (CP) idealized conformity measure is

Q X is the marginal distribution of Q on X; -N-optimal if E (x,τ )∼Q X ×U |Γ A (x)| ≤ E (x,
A((x, y), Q) := Q(y | x).
An idealized conformity measure A is a (label-conditional) refinement of an idealized conformity measure B if

B((x 1 , y)) < B((x 2 , y)) =⇒ A((x 1 , y)) < A((x 2 , y)) (5)
for all x 1 , x 2 ∈ Z and all y ∈ Y. (Notice that this definition, being labelconditional, is different from the one given in [START_REF] Vovk | Criteria of efficiency for conformal prediction[END_REF].) Let R(CP) be the set of all refinements of the CP idealized conformity measure. If C is a criterion of efficiency (one of the four discussed above), we let O(C) stand for the set of all C-optimal idealized conformity measures.

Theorem 1. O(S) = O(OF) = O(N) = O(OE) = R(CP).
Proof. We start from proving R(CP) = O(N). Fix a significance level . A smoothed confidence predictor at level is defined as a random set of observations (x, y) ∈ Z; in other words, to each observation (x, y) is assigned the probability P (x, y) that the observation will be outside the prediction set. Under the restriction that the sum of the probabilities Q(x, y) of observations (x, y) outside the prediction set (defined as x Q(x, y)P (x, y) in the smoothed case) is bounded by Q Y (y) for a fixed y, the N criterion requires us to make the sum of Q X (x) for (x, y) outside the prediction set (defined as x Q X P (x, y) in the smoothed case) as large as possible. It is clear that the set should consist of the observations with the smallest Q(y | x) (by the usual Neyman-Pearson argument: cf. [START_REF] Lehmann | Testing Statistical Hypotheses[END_REF], Sect. 3.2).

Next we show that O(N) ⊆ O(S). Let an idealized conformity measure A be N-optimal. By definition,

E x,τ |Γ A (x)| ≤ E x,τ |Γ B (x)|
for any idealized conformity measure B and any significance level . Integrating over ∈ (0, 1) and swapping the order of integrals and expectations,

E x,τ 1 0 |Γ A (x)| d ≤ E x,τ 1 0 |Γ B (x)| d . (6)
Since

|Γ (x)| = y∈Y 1 {p(x,y)> } ,
we can rewrite (6), after swapping the order of summation and integration, as where we have used the fact that p(x, y) is distributed uniformly on [0, 1] when ((x, y), τ ) ∼ Q × U (see [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF] and [START_REF] Vovk | Criteria of efficiency for conformal prediction[END_REF]). Finally, we notice that O(N) = O(OE). Indeed, for any significance level ,

E x,τ y∈Y 1 0 1 {p A (x,y)> } d ≤ E x,τ y∈Y 1 0 1 {p B (x,y)> } d .
E x,τ |Γ (x)| = E (x,y),τ |Γ (x) \ {y}| + (1 -),
again using the fact that p(x, y) is distributed uniformly on [0, 1] and so P (x,y),τ (y ∈ Γ (x)) = 1 -.

Criteria of efficiency for probabilistic predictors

Given a training set (z 1 , . . . , z l ) and a test object x, a probabilistic predictor outputs a probability measure P ∈ P(Y), which is interpreted as its probabilistic prediction for the label y of x; we let P(Y) stand for the set of all probability measures on Y. The two standard way of measuring the performance of P on the actual label y are the logarithmic (or log) loss -ln P ({y}) and the Brier loss

y ∈Y 1 {y =y} -P ({y }) 2 ,
where 1 E stands for the indicator of an event E: 1 E = 0 if E happens and 1 E = 0 otherwise. The efficiency of probabilistic predictors will be measured by these two loss functions. Suppose we have a test sequence (z l+1 , . . . , z l+k ), where z i = (x i , y i ) for i = l + 1, . . . , l + k, and we want to evaluate the performance of a probabilistic predictor (trained on a training sequence z 1 , . . . , z l ) on it. In the next section we will use the average log loss

1 k |Y| l+k i=l+1 y ∈Y 1 {y =yi} -P i ({y }) 2 ,
where P i ∈ P(Y) is the probabilistic prediction for x i . Notice that in the binary case, |Y| = 2, the average log loss coincides with the mean log error (used in, e.g., [START_REF] Vovk | Venn-Abers predictors[END_REF], [START_REF] Vovk | Criteria of efficiency for conformal prediction[END_REF]) and the standardized Brier loss coincides with the root mean square error (used in, e.g., [START_REF] Vovk | Venn-Abers predictors[END_REF], ( 13)).

Calibration of p-values into conditional probabilities

The argument of this section will be somewhat heuristic, and we will not try to formalize it in this paper. Fix y ∈ Y. Suppose that q := P (y | x) has an absolutely continuous distribution with density f when x ∼ Q X . (In other words, f is the density of the image of Q X under the mapping x → P (y | x).) For the CP idealized conformity measure, we can rewrite (3) as

p(q) := q 0 q f (q )dq D , (7)
where D := Q Y ({y}); alternatively, we can set D := 1 0 q f (q )dq to the normalizing constant ensuring that p(1) = 1. To see how ( 7) is a special case of (3) for the CP idealized conformity measure, notice that the probability that Y = y and P (Y | X) ∈ (q , q + dq ), where (X, Y ) ∼ f , is q f (q )dq . In [START_REF] Sellke | Calibration of p-values for testing precise null hypotheses[END_REF] we write p(q) rather than p y since p y depends on y only via q.

We are more interested in the inverse function q(p), which is defined by the condition p = q(p) 0 q f (q )dq D .

When q ∼ f , we have

P(p(q) ≤ a) = P(q ≤ q(a)) = q(a) 0 f (q )dq .
Therefore, when q ∼ f , we have P(a ≤ p(q) ≤ a + da) = q(a+da) q(a)

f (q )dq ≈ 1 q(a) q(a+da) q(a) q f (q )dq = Dda q(a) , and so

q(c) ≈ D P(c ≤ p(q) ≤ c + dc) dc .
This gives rise to the algorithm given as Algorithm 1, which uses real p-values (2) instead of the ideal p-values (3). The algorithm is transductive in that it uses a training sequence of labelled observations and a calibration sequence of unlabelled objects (in the next section we use the test sequence as the calibration sequence); the latter is used for calibrating p-values into conditional probabilities. Given all the p-values for the calibration sequence with postulated label y, find the corresponding antitonic density g(p) (remember that the function q(p) is known to be monotonic, namely isotonic) using Grenander's estimator (see [START_REF] Grenander | On the theory of mortality measurement[END_REF] or, e.g., [START_REF] Devroye | A Course in Density Estimation[END_REF], Chap. 8). Use D/g(p) as the calibration function, where D := g(1) is chosen in such a way that a p-value of 1 is calibrated into a conditional probability of 1. (Alternatively, we could set D to the fraction of observations labelled as y in the training sequence; this approximates setting D := Q Y ({y}).) The probabilities produced by this procedure are not guaranteed to lead to a probability measure: the sum over y can be different from 1 (and this phenomenon has been observed in our experiments). Therefore, in the last line of Algorithm 1 we normalize the calibrated p-values to obtain genuine probabilities.

Experiments

In our experiments we use the standard USPS data set of hand-written digits. The size of the training set is 7291, and the size of the test set is 2007; however, instead of using the original split of the data into the two parts, we randomly split all available data (the union of the original training and test sets) into a training set of size 7291 and test set of size 2007. (Therefore, our results somewhat depend on the seed used by the random number generator, but the dependence is minor and does not affect our conclusions at all; we always report results for seed 0.) A powerful algorithm for the USPS data set is the 1-Nearest Neighbour (1-NN) algorithm using tangent distance [START_REF] Simard | Efficient pattern recognition using a new transformation distance[END_REF]. However, it is not obvious how this algorithm could be transformed into a probabilistic predictor. On the other hand, there is a very natural and standard way of extracting probabilities from support vector machines, which we will refer to it as Platt's algorithm in this paper: it is the combination of the method proposed by Platt [START_REF] Platt | Probabilities for SV machines[END_REF] with pairwise coupling [START_REF] Wu | Probability estimates for multi-class classification by pairwise coupling[END_REF] (unlike our algorithm, which is applicable to multi-class problems directly, Platt's method is directly applicable only to binary problems). In this section we will apply our method to the 1-NN algorithm with tangent distance and compare the results to Platt's algorithm as implemented in the function svm from the e1071 R package (for our multi-class problem this function calculates probabilities using the combination of Platt's binary method and pairwise coupling).

There is a standard way of turning a distance into a conformal predictor ([13], Sect. 3.1): namely, the conformity score α i of the ith observation in a sequence of observations can be defined as min j:yj =yi d(x i , x j ) min j =i:yj =yi d(x i , x j ) ,

where d is the distance; the intuition is that an object is considered conforming if it is close to an object labelled in the same way and far from any object labelled in a different way. Table 1 compares the performance of the conformal-type probabilistic predictor based on the 1-NN conformity measure [START_REF] Simard | Efficient pattern recognition using a new transformation distance[END_REF], where d is tangent distance, with the performance of Platt's algorithm with the optimal values of its parameters. The conformal predictor is parameter-free but Platt's algorithm depends on the choice of the kernel. We chose the polynomial kernel of degree 3 (since it is known to produce the best results: see [START_REF] Vapnik | Statistical Learning Theory[END_REF], Sect. 12.2) and the cost parameter C := 2.9 in the case of the average log loss and C := 3.4 in the case of the standardized Brier loss (the optimal values in our experiments). (Reporting the performance of Platt's algorithm with optimal parameter values may look like data snooping, but it is fine in this context since we are helping our competitor.) Table 2 reports the performance of Platt's algorithm as function of the degree of the polynomial kernel with the cost parameter set at C := 10 (the dependence on C is relatively mild, and C = 10 gives good performance for all degrees that we consider).

  τ )∼Q X ×U |Γ B (x)| for any idealized conformity measure B and any significance level ; -OF-optimal ifE ((x,y),τ )∼Q×U y =y p A (x, y ) ≤ E ((x,y),τ )∼Q×U y =y p A (x, y )for any idealized conformity measure B; -OE-optimal ifE ((x,y),τ )∼Q×U |Γ A (x) \ {y}| ≤ E ((x,y),τ )∼Q×U |Γ B (x) \ {y}|for any idealized conformity measure B and any significance level .

Since 1 0 1

 11 {p(x,y)> } d = p(x, y), we finally obtain E x,τ y∈Y p A (x, y) ≤ E x,τ y∈Y p B (x, y). Since this holds for any idealized conformity measure B, A is S-optimal. The argument in the previous paragraph in fact shows that O(S) = O(N) = R(CP). The equality O(S) = O(OF) follows from E x,τ y p(x, y) = E (x,y),τ

Algorithm 1

 1 Conformal-type probabilistic predictorInput: training sequence (z1, . . . , z l ) ∈ Z l Input: calibration sequence (x l+1 , . . . , x l+k ) ∈ X k Input: test object x0 Output: probabilistic prediction P ∈ P(Y) for the label of x0 for y ∈ Y do for each xi in the calibration sequence find the p-value p y i by (2) (with l + i in place of l + 1) let gy be the antitonic density on [0, 1] fitted to p y l+1 , . . . , p y l+k find the p-value p y 0 by (2) (with 0 in place of l + 1) for each y ∈ Y, set P ({y}) := gy(1)/gy(p y 0 ) end for set P ({y}) := P ({y})/ y P ({y }) for each y ∈ Y and the standardized Brier loss

Table 1 .

 1 The performance of the two algorithms, Platt's (with the optimal values of parameters) and the conformal-type probabilistic predictor based on 1-Nearest Neighbour with tangent distance

	algorithm average log loss standardized Brier loss
	optimized Platt	0.06431	0.05089
	conformal-type 1-NN	0.04958	0.04359

Table 2 .

 2 The performance of Platt's algorithm with the polynomial kernels of various degrees for the cost parameter C = 10

	degree average log loss standardized Brier loss
	1	0.12681	0.07342
	2	0.09967	0.06109
	3	0.06855	0.05237
	4	0.11041	0.06227
	5	0.09794	0.06040

Acknowledgments. We thank the reviewer for useful comments. In our experiments we used the R package e1071 (by David Meyer, Evgenia Dimitriadou, Kurt Hornik, Andreas Weingessel, Friedrich Leisch, Chih-Chung Chang, and Chih-Chen Lin) and the implementation of tangent distance by Daniel Keysers. This work was partially supported by EPSRC (grant EP/K033344/1, first author) and Royal Holloway, University of London (third author).


Lars Carlsson
email: lars.a.carlsson@astrazeneca.com

Martin Eklund
email: martin.eklund@farmbio.uu.seulfn@lundbeck.com

Ulf Norinder

Aggregated conformal prediction

We present the aggregated conformal predictor (ACP), an extension to the traditional inductive conformal prediction (ICP) where several inductive conformal predictors are applied on the same training set and their individual predictions are aggregated to form a single prediction on an example. The results from applying ACP on two pharmaceutical data sets (CDK5 and GNRHR) indicate that the ACP has advantages over traditional ICP. ACP reduces the variance of the prediction region estimates and improves efficiency. Still, it is more conservative in terms of validity than ICP, indicating that there is room for further improvement of efficiency without compromising validity.

Introduction

Quantitative Structure-Activity Relationship (QSAR) modeling for predicting properties, e.g. solubility or toxicity, of chemical compounds using statistical learning techniques is a widespread approach within the pharmaceutical industry to prioritize compounds for experimental testing or to alert for potential toxicity. In making informed decisions based on predictions from QSAR models, the confidence in such predictions is of vital importance and conformal predictors have been successfully applied to the drug discovery setting [START_REF] Eklund | Application of Conformal Prediction in QSAR[END_REF]. In particular, we have shown that a Mondrian inductive conformal predictor is efficient (i.e. informative) and almost valid when applied to binary categorical data; exact validity was not achieved due to deviations from the exchangeability assumption [START_REF] Eklund | The application of conformal prediction to the drug discovery process[END_REF]. Furthermore, we have demonstrated that an inductive conformal predictor (ICP) applied to regression data from the pharmaceutical industry was valid. However, for regression the prediction regions were in many cases as wide or wider than the possible ranges of the true responses (i.e. the range of the experimental assay that generates the true label). An important problem is thus to improve the efficiency of the ICPs when applied in the QSAR domain.

There are many ways to improve efficiency of conformal predictions, for example by using improved nonconformity scores, choosing different machine-learning methods, or using a transductive approach [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF]. The transductive approach is the most appealing in terms of validity, but is often computationally costly and the nonconformity scores can be difficult to compute.

Another interesting approach to improve efficiency is the cross conformal predictor (CCP) [START_REF] Vovk | Cross-conformal predictors[END_REF]. Here the training data is divided into separate non-overlapping folds and each fold is used as a calibration set and the remainder of the data is used as a proper training set. This division allows for more data to be used for calibration and p-values are averaged over all folds. Similarly, the bootstrap conformal predictor (BCP) bootstraps datasets and uses the out-of-bag examples as a calibration set. p-values are then averages across all bootstrap replications.

In this paper we attempt to generalize the BCP and the CCP in what we term the aggregated inductive conformal predictor (ACP). We will empirically assess the ACP using data from the pharmaceutical domain and we show through a theoretical argument and experiments that ACP seems to have advantages over the standard ICP.

Aggregated conformal predictor

Consider the standard prerequisite for a description of a conformal predictor (CP), a bag of examples {z 1 , . . . , z i , . . . , z l } drawn from an exchangeable distribution Q. Each example z i = (x i , y i ) can be described by its object x i ∈ X and its label y i ∈ Y. The labels can be either categorical or continuous. For an inductive conformal predictor (ICP) the bag {z i } is partitioned into two different bags, one holding the proper training examples {z 1 , . . . , z m } and the other holding the calibration examples {z m+1 , . . . , z l }. The ICP p-value is then computed as

p = | {j = m + 1, ..., l : α j ≥ α l+1 } | l -m + 1
The prediction region of an ICP is determined by the "borderline" p-value, p t , i.e. the smallest value p can obtain and still satisfy p > . We can thus view p t as the th sample quantile (estimated from above)

p t = U -1 l-m ( ),
where U l-m is the empirical cumulative probability distribution of the p-values defined by

U l-m (p) = 1 l -m + 1 l j=m+1 I(p j < p),
where I(•) = 1 if p j < p and I(•) = 0 otherwise. We now introduce definitions of Exchangeable resampling and Consistent resampling, after which we define what we mean by an aggregated conformal predictor (ACP).

Definition 1 (Exchangeable resampling). Let {z * 1 , . . . , z * n } be a bag of examples resampled from the empirical distribution Q l . We call this resampling

exchangeable if P {(z * 1 , . . . , z * n )} = P {(z * π(1) , . . . , z * π(n) )},
where π is any permutation of {1, ..., n}.

Definition 2 (Consistent resampling). Let T = T (z 1 , ..., z l , Q) be a statistic and T * = T (z * 1 , ..., z * n , Q l ) be an exchangeably resampled version of T . Further, let G l and G * l be the probability distributions of T and T * , respectively. We call the sampling process consistent (with respect to T ) if

sup z | G l -G * l |→ 0 as l → ∞ and n → ∞.
Definition 3 (ACP: Aggregated Conformal Predictor). The following procedure is repeated B times, for b = 1, ..., B: Resample a bag {z * 1 , . . . , z * n b } of examples from {z 1 , . . . , z l } using a consistent resampling procedure with respect to α t . Compute the ICP p-value using the resampled bag,

p * b = | {j = m b + 1, ..., n b : α * j ≥ α * n b +1 } | n b -m b + 1 , (1)
where α * j are the nonconformity scores computed using {z * 1 , . . . , z * n b } (m b and n b are indexed with b to make explicit that they may differ for different values of b). We define the ACP p-value as

p B = 1 B B b=1 p * b ( 2
)
and the corresponding prediction region as

Γ (z 1 , . . . , z l , x l+1 ) := {y|p B > }. (3)
A smoothed ACP can be defined analogously. We note that the cross-conformal predictor and the bootstrapped conformal predictor suggested by Vovk [START_REF] Vovk | Cross-conformal predictor[END_REF] are two examples of ACPs.

Proposition 1. The aggregated conformal predictor is conservatively valid.

Proof. Since we use an exchangeable resampling procedure to construct the ACP and since an ICP is conservatively valid (Proposition 4.1, [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF]), each resampled ICP in the ACP is conservatively valid (by symmetry). From this follows that the ACP also is conservatively valid.

Remark 1. Proposition 1 only holds unconditionally. The situation is different conditional on the particular dataset we have observed.

How does the ACP improve on the ICP?

For a p in an ICP, p t is a hard threshold in the sense that the label y corresponding to p is either inside or outside the prediction region. Heuristically, the ACP averages over thresholds varying around p t (since p * t based on a resampled bag {z * 1 , ..., z * l } fluctuates around p t ), resulting in a smoothed threshold estimate with decreased variance compared to the estimate p t .

We can use the method used by Bühlmann and Yu [START_REF] Bühlmann | Analyzing bagging[END_REF] to analyze this in a bit more detail. Consider the function

δ(p) = I(p t < p),
which indicates whether a p is smaller or larger than the borderline value p t in an ICP (and thus if its corresponding label y either is inside or outside the prediction region). The sample quantile estimate p t follows a normal distribution with mean = U -1 ( ) and variance

σ 2 = (1 -) (l -m + 1) [f (q)] 2 = 1 - (l -m + 1)
, where F is the population cumulative distribution function with density function f [START_REF] Ruppert | Statistics and Data Analysis for Financial Engineering[END_REF]. For a p in the neighborhood of

p = p(c) = + cσ √ l -m + 1 (4)
we have the approximation

δ(p(c)) ≈ I(W < c), W ∼ N (0, 1), ( 5
)
where W is the limiting random quantity from the asymptotic distribution of p t (because of the construction of p(c) in Equation ( 4)). For a fixed c, this is a hard threshold function of W . It follows that

E [δ(p(c))] → P (W < c) = Φ(c) as l → ∞ and l -m → ∞ Var [δ(p(c))] → Φ(c)(1 -Φ(c)) as l → ∞ and l -m → ∞, (6)
where Φ(•) is the cumulative distribution function of the standard normal distribution. Note that the variance does not converge to zero; δ(p(c)) assumes the values 0 and 1 with a positive probability even as l tends to infinity. However, for the ACP the situation looks different,

δ B (p(c)) = E * [I(p * t ≤ p(c))] = E * I √ l -m + 1(p * t -p t )/σ ≤ √ l -m + 1(p(c) -p t )/σ = Φ( √ l -m + 1(p(c) -p t )) + o p (1) ≈ Φ(c -W ), W ∼ N (0, 1).
where the first approximation over the second equal sign follows because we by definition of the ACP have a consistent resampling process. Comparing with Equation ( 5) for an ICP, the ACP produces a smoothed decision function of Z and therefore reduces variance. Again, following Bühlmann and Yu [START_REF] Bühlmann | Analyzing bagging[END_REF], we can study the case p = p(0) = , i.e. when we are right at the population threshold and therefore has maximum variance. Then

δ B (p(0)) → Φ(-W ) ∼ U [0, 1]
and, therefore,

E [δ B (p(0))] → E [U ] = 1/2 as l → ∞ Var [δ B (p(0))] → Var(U ) = 1/12 as l → ∞. (7)
Comparing Equation [START_REF] Bühlmann | Analyzing bagging[END_REF] to Equation ( 7), we see that the variance is reduced to one third for ACP compared to ICP.

Empirical results of ACP

We used two different machine learning methods; the support vector machine (SVM) [START_REF] Vapnik | Statistical learning theory[END_REF] implemented in the Java library version of libsvm [START_REF] Chang | LIBSVM: A library for support vector machines[END_REF] with a Gaussian radial basis kernel function

K(x, x ) = exp (-γ x -x 2 ),
with γ = 0.002 and C = 50, and Random Forest (RF) [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF], for which the default settings were used and each ensemble contained 100 trees. Following Equation ( 16) in [START_REF] Papadopoulos | Reliable prediction intervals with regression neural networks[END_REF], we defined the nonconformity measure used in combination with SVM according to

α i = |y i -ŷi | exp (μ i ) , ( 8
)
where μi is the prediction of the value ln (|y i -ŷi |) produced by a support vector regression machine trained on the proper training sets. After training the underlying SVM of the ICP, we calculate the residuals |y j -ŷj | for all proper training examples j = 1, ..., m and train an SVM on the pairs (x i , ln(|y i -ŷi |)). Measure (8) normalizes the absolute prediction error with the predicted accuracy of the SVM on a given example. The nonconformity measure used with RF was

α i = |y i -ŷi | νi ,
where νi is the RF prediction of the value |y i -ŷi | with the same settings as for the other RF model. Two public dataset (CDK5 and GNRHR) from the pharmaceutical domain was used [START_REF] Chen | Beyond the scope of free-wilson analysis: Building interpretable qsar models with machine learning algorithms[END_REF]. The datasets consist of 230 and 198 examples, respectively. Each example (chemical compound) was described (characterised) by so-called signature descriptors [START_REF] Faulon | The signature molecular descriptor. 4. Canonizing molecules using extended valence sequences[END_REF] in the same way as described in [START_REF] Eklund | The application of conformal prediction to the drug discovery process[END_REF].

The data was randomly split into two parts: A training set (80% of the original data) and a working set (20%). This procedure was repeated 50 times as to generate 50 training and working set, respectively. Furthermore, each training set was then, subsequently, randomly split into a proper training set (70% of the training set) and a calibration set (30%), similar to the 2:1 recommendation in [START_REF] Vovk | Cross-conformal predictor[END_REF]. This random selection of proper training and calibration examples was, in turn, repeated 100 times enabling the construction of 100 inductive conformal predictors for each working set. This sampling procedure is often called the m -n sampling or non-replacement subsampling in the bootstrap literature and consistent with Definition 2 [START_REF] Politis | Subsampling[END_REF].

The results are presented in Tables 1-8 and in Figures 1 and2. Table 5. The efficiency for an ACP and an ICP at one significance level for the 50 runs on the CDK5 data using SVM. Here, Q1 and Q3 are the lower and upper quartile, respectively.

Method

Discussion

The results presented in Tables 1-8 and in Figures 1 and2 indicate that the ACP methodology has advantages over traditional ICP. The former adds stability and robustness to the predictions, which is particularly clear in Figures 12where the variance in the prediction ranges is smaller than from an individual ICP. This is of considerable importance within the pharmaceutical domain where precision as well as robustness in predictions are key elements for successful application in ongoing discovery projects. Although the traditional ICP is valid on average, the variance is very large. This means that there is a relatively large proportion of very tight predictions regions that in fact are far from valid (for example, for the CDK5 data, an ICP with a confidence of 80% produces 46% error in the quartile with tightest prediction regions; the corresponding figure for the ACP is 21%. These tight prediction regions are offset by some prediction regions that are very wide (as wide or wider than the possible range of the response value) that are always valid, which produces a conformal predictor that is valid on average. Neither the too tight prediction regions that cannot be trusted nor non-informative regions are helpful for the researcher using the model. The results in Tables 12345678show that the ACP (as used in this paper) is too conservative on datasets with a relatively small number of examples, which is clear from Equations (1) and (2). This indicates that improved efficiency can be achieved without compromising the validity of the ACP, e.g. by more clever resampling (a large body of literature exists that address this problem, see e.g. [START_REF] Egloff | Quantile estimation with adaptive importance sampling[END_REF]) or smoothing of p-values on either side of . To conclude: We have introduced the ACP, a generalization of the BCP and CCP introduced by Vovk [5] and we have shown that it improves on classical ICP by reducing the variance in the estimated prediction regions (analogously to how a bagged predictor improves the prediction by reducing variance). ACP seems to represent a pragmatic and useful way forward for obtaining models and predictions with good precision and small prediction ranges.

Future ways to develop the ACP include to (i) study ACP for categorical labels (e.g. aggregation through voting); (ii) other resampling schemas.

Fig. 1 .

 1 Fig. 1. The distribution of prediction interval size for an ACP and a randomly sampled ICP for all 50 runs on the CDK5 data. The results are shown at = 0.2.

Fig. 2 .

 2 Fig. 2. The distribution of prediction interval size for an ACP and a randomly sampled ICP for all 50 runs on the GNRHR data. The results are shown at = 0.2.

Table 1 .

 1 The fraction of true labels within the prediction ranges for an ACP and a, for each run, randomly selected ICP at different significance levels for the 50 runs on the CDK data using SVM. Here, Q1 and Q3 are the lower and upper quartile, respectively.

		min	Q1 median mean	Q3 max
	ACP	0.4 0.6522 0.7391 0.7717 0.7717 0.8043 0.8913
	ICP	0.4 0.5100 0.5622 0.5997 0.5990 0.6303 0.6826
	ACP	0.3 0.7174 0.8261 0.8478 0.8474 0.8696 0.9348
	ICP	0.3 0.6191 0.6535 0.6879 0.6896 0.7202 0.7876
	ACP	0.2 0.8043 0.8913 0.9130 0.9170 0.9565 0.9783
	ICP	0.2 0.7296 0.7712 0.7923 0.7989 0.8282 0.8848
	ACP	0.1 0.8913 0.9565 0.9565 0.9643 0.9783 1.0000
	ICP	0.1 0.8378 0.8626 0.8933 0.8890 0.9101 0.9537
	Method	min	Q1 median mean	Q3 max
	ACP	0.4 0.3590 0.5962 0.6667 0.6505 0.7179 0.8205
	ICP	0.4 0.4297 0.5560 0.5960 0.5947 0.6569 0.7321
	ACP	0.3 0.5128 0.6731 0.7179 0.7227 0.7628 0.8718
	ICP	0.3 0.5215 0.6428 0.6736 0.6750 0.7208 0.7926
	ACP	0.2 0.7179 0.7949 0.8205 0.8286 0.8462 0.9744
	ICP	0.2 0.6782 0.7466 0.7749 0.7791 0.8126 0.8797
	ACP	0.1 0.8205 0.8974 0.9231 0.9258 0.9487 1.0000
	ICP	0.1 0.8072 0.8633 0.8745 0.8835 0.9104 0.9605

Table 2 .

 2 The fraction of true labels within the prediction ranges for an ACP and a, for each run, randomly selected ICP at different significance levels for the 50 runs on the GNRHR data using SVM. Here, Q1 and Q3 are the lower and upper quartile, respectively.

	Method	min	Q1 median mean	Q3 max
	ACP	0.2 0.7174 0.8478 0.8696 0.8696 0.9130 0.9565
	ICP	0.2 0.6926 0.7786 0.8071 0.7982 0.8249 0.8930

Table 3 .

 3 The fraction of true labels within the prediction ranges for an ACP and a, for each run, randomly selected ICP at one significance level for the 50 runs on the CDK data using RF. Here, Q1 and Q3 are the lower and upper quartile, respectively.

	Method	min	Q1 median mean	Q3 max
	ACP	0.2 0.7179 0.8205 0.8718 0.8564 0.8974 1.0000
	ICP	0.2 0.6854 0.7706 0.7919 0.8013 0.8271 0.9162

Table 4 .

 4 The fraction of true labels within the prediction ranges for an ACP and a, for each run, randomly selected ICP at one significance level for the 50 runs on the GNRHR data using RF. Here, Q1 and Q3 are the lower and upper quartile, respectively.

	Method	min	Q1 median mean	Q3	max
	ACP	0.2 0.5874 1.2190 1.6880 3.2000 3.3220 67.6700
	ICP	0.2 0.1550 0.9128 1.2620 3.5760 1.8970 398.3000

Table 6 .

 6 The efficiency for an ACP and an ICP at one significance level for the 50 runs on the CDK5 data using RF. Here, Q1 and Q3 are the lower and upper quartile, respectively.

	Method		min	Q1 median mean	Q3 max
	ACP	0.2 0.4752 1.0060 1.2640 1.3750 1.5870 3.9780
	ICP	0.2 0.0653 0.8659 1.1830 1.3680 1.6380 8.5950
	Method	min Q1 median mean Q3 max
	ACP		0.2 1.2770 1.625 1.733 1.907 1.903 25.830
	ICP		0.2 0.5486 1.457 1.704 1.957 2.000 62.400

Table 7 .

 7 The efficiency for an ACP and an ICP at one significance level for the 50 runs on the GNRHR data using SVM. Here, Q1 and Q3 are the lower and upper quartile, respectively.

	Method	min	Q1 median mean	Q3	max
	ACP	0.2 0.7481 1.6050 1.9620 2.1290 2.5490 6.3040
	ICP	0.2 0.07138 1.3940 1.8230 2.0980 2.4840 10.3100

Table 8 .

 8 The efficency for an ACP and an ICP at one significance level for the 50 runs on the GNRHR data using RF. Here, Q1 and Q3 are the lower and upper quartile, respectively.


Harris Papadopoulos
email: h.papadopoulos@frederick.ac.cy

A Cross-Conformal Predictor for Multi-label Classification

Unlike the typical classification setting where each instance is associated with a single class, in multi-label learning each instance is associated with multiple classes simultaneously. Therefore the learning task in this setting is to predict the subset of classes to which each instance belongs. This work examines the application of a recently developed framework called Conformal Prediction (CP) to the multi-label learning setting. CP complements the predictions of machine learning algorithms with reliable measures of confidence. As a result the proposed approach instead of just predicting the most likely subset of classes for a new unseen instance, also indicates the likelihood of each predicted subset being correct. This additional information is especially valuable in the multi-label setting where the overall uncertainty is extremely high.

Introduction

Most machine learning research on classification deals with problems in which each instance is associated with a single class y from a set of classes {Y 1 , . . . , Y c }. As opposed to this standard setting, in multi-label classification each instance can belong to multiple classes, so that each instance is associated with a set of classes ψ ⊆ {Y 1 , . . . , Y c }, called a labelset. There are many real-world problems in which such a setting is natural. For instance in the categorization of news articles an article discussing the positions of political parties on the educational system of a country can be classified as both politics and education. Although until recently the main multi-label classification application was the categorization of textual data, in the last few years an increasing number of new applications started to attract the attention of more researchers to this setting. Such applications include the semantic annotation of images and videos, the categorization of music into emotions, functional genomics, proteomics and directed marketing.

As a result of the increasing attention to the multi-label classification setting, many new machine learning techniques have been recently developed to deal with problems of this type, see e.g. [START_REF] Elisseeff | A kernel method for multi-labelled classification[END_REF][START_REF] Zhang | ML-RBF: RBF neural networks for multi-label learning[END_REF][START_REF] Zhang | Feature selection for multi-label naive bayes classification[END_REF][START_REF] Zhang | Multi-label neural networks with applications to functional genomics and text categorization[END_REF][START_REF] Zhang | ML-kNN: A lazy learning approach to multi-label learning[END_REF]. However, like most machine learning methods, these techniques do not produce any indication about the likelihood of each of their predicted labelsets being correct. Such an indication though can be very helpful in deciding how much to rely on each prediction, especially since the certainty of predictions may vary to a big degree between instances. To address this problem this paper examines the application of a recently developed framework for providing reliable confidence measures to predictions, called Conformal Prediction (CP) [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF], to the multi-label classification setting. Specifically it follows the newly proposed Cross-Conformal Prediction (CCP) [START_REF] Vovk | Cross-conformal predictors[END_REF] version of the framework, which allows it to overcome the prohibitively large computational overhead of the original CP. The proposed approach computes a p-value for each of the possible labelsets, which can be used either for accompanying each prediction with confidence measures that indicate its likelihood of being correct, or for producing sets of labelsets that are guaranteed to contain the true labelset with a frequency equal to or higher than a required level of confidence.

In the remaining paper, a description of the general idea behind CP and of the CCP version of the framework, is first provided in Section 2. Section 3 defines the proposed approach while Section 4 presents the experiments performed and the obtained results. Finally Section 5 gives the conclusions of this work.

Conformal and Cross-Conformal Prediction

Typically in classification we are given a set of training examples {z 1 , . . . , z l }, where each z i ∈ Z is a pair (x i , y i ) consisting of a vector of attributes x i ∈ IR d and the classification y i ∈ {Y 1 , . . . , Y c }. We are also given a new unclassified example x l+1 and our task is to state something about our confidence in each possible classification of x l+1 without assuming anything more than that all (x i , y i ), i = 1, 2, . . . , are independent and identically distributed.

The idea behind CP is to assume every possible classification Y j of the example x l+1 and check how likely it is that the extended set of examples {(x 1 , y 1 ), . . . , (x l , y l ), (x l+1 , Y j )}

(1) is i.i.d. This in effect will correspond to the likelihood of Y j being the true label of the example x l+1 since this is the only unknown value in [START_REF] Boutell | Learning multi-label scene classification[END_REF]. First a function A called nonconformity measure is used to map each pair (x i , y i ) in (1) to a numerical score

α i = A({(x 1 , y 1 ), . . . , (x l , y l ), (x l+1 , Y j )}, (x i , y i )), i = 1, . . . , l, (2a)
α Yj l+1 = A({(x 1 , y 1 ), . . . , (x l , y l ), (x l+1 , Y j )}, (x l+1 , Y j )), (2b)
called the nonconformity score of instance i. This score indicates how nonconforming, or strange, it is for i to belong in [START_REF] Boutell | Learning multi-label scene classification[END_REF]. In effect the nonconformity measure is based on a conventional machine learning algorithm, called the underlying algorithm of the corresponding CP and measures the degree of disagreement between the actual label y i and the prediction ŷi of the underlying algorithm, after being trained on (1). The nonconformity measure for multi-label learning used in this work is defined in Section 3.

The nonconformity score α Yj l+1 is then compared to the nonconformity scores of all other examples to find out how unusual (x l+1 , Y j ) is according to the nonconformity measure used. This comparison is performed with the function p((x 1 , y 1 ), . . . , (x l , y l ), (x l+1 , Y j )) = |{i = 1, . . . , l, l + 1 :

α Yj i ≥ α Yj l+1 }| l + 1 , (3)
the output of which is called the p-value of Y j , also denoted as p(Y j ). An important property of (3) is that ∀δ ∈ [0, 1] and for all probability distributions P on Z, P l+1 {((x 1 , y 1 ), . . . , (x l , y l ), (x l+1 , y l+1 )) :

p(y l+1 ) ≤ δ} ≤ δ; (4)
a proof can be found in [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF]. According to this property, if p(Y j ) is under some very low threshold, say 0.05, this means that Y j is highly unlikely as the probability of such an event is at most 5% if (1) is i.i.d.

There are two standard ways to use the p-values of all possible classifications for producing the output of a CP:

-Predict the classification with the highest p-value and output one minus the second highest p-value as confidence to this prediction and the p-value of the predicted classification (i.e. the highest p-value) as credibility.

-Given a confidence level 1 -δ, output the prediction set {Y j : p(Y j ) > δ}.

In the first case, confidence is an indication of how likely the prediction is of being correct compared to all other possible classifications, whereas credibility indicates how suitable the training set is for the particular instance; specifically a very low credibility value indicates that the particular instance does not seem to belong to any of the possible classifications. In the second case, the prediction set will not contain the true label of the instance with at most δ probability.

The important drawback of the above process is that since the last example in (1) changes for every possible classification, the underlying algorithm needs to be trained c times. Moreover the whole process needs to be repeated for every test example. This makes it extremely computationally inefficient in many cases and especially in a multi-label setting where there are 2 n possible labelsets for n classes, or 2 n -1 if we exclude the empty labelset.

To overcome this computational inefficiency problem an inductive version of the framework was proposed in [START_REF] Papadopoulos | Inductive confidence machines for regression[END_REF] and [START_REF] Papadopoulos | Qualified predictions for large data sets in the case of pattern recognition[END_REF] called Inductive Conformal Prediction (ICP). ICP is based on the same theoretical foundations described above, but follows a modified version of the approach, which allows it to train the underlying algorithm only once. This is achieved by dividing the training set into two smaller sets, the proper training set and the calibration set. The proper training set is then used for training the underlying algorithm of the ICP and only the examples in the calibration set are used for calculating the p-value of each possible classification for every test example.

Although ICP is much more computationally efficient than the original CP approach, the fact that it does not use the whole training set for training the underlying algorithm and for calculating its p-values results in lower informational efficiency. That is the resulting prediction sets might be larger than the ones produced by the original CP approach. Cross-Conformal Prediction, which was recently proposed in [START_REF] Vovk | Cross-conformal predictors[END_REF], tries to overcome this problem by combining ICP with cross-validation. Specifically, CCP partitions the training set in K subsets (folds) S 1 , . . . , S K and calculates the nonconformity scores of the examples in each subset S k and of (x l+1 , Y j ) for each possible classification Y j as

α i = A(∪ m =k S m , (x i , y i )), i ∈ S k , m = 1, . . . , K, (5a)
α

Yj ,k l+1 = A(∪ m =k S m , (x l+1 , Y j )), m = 1, . . . , K, (5b)
where A is the given nonconformity measure. Note that for (x l+1 , Y j ) K nonconformity scores α Yj ,k l+1 , k = 1, . . . , K are calculated, one with each of the K folds. Now the p-value for each possible classification Y j is computed as

p(Y j ) = K k=1 |{(x i , y i ) ∈ S k : α i ≥ α Yj ,k l+1 }| + 1 l + 1 . ( 6
)
The CCP version of the framework was chosen to be followed here due to its big advantage in computational efficiency over CP, since it needs to train the underlying algorithm only K times, and its advantage over ICP since it utilizes the whole training set for producing its p-values. As opposed to CP and ICP, the validity of which has been proven theoretically, at the moment there are no theoretical results about the validity of CCP. However in [START_REF] Vovk | Cross-conformal predictors[END_REF] its outputs have been shown to be empirically valid. The same is shown in the experimental results of this work, presented in Section 4.

ML-RBF Cross-Conformal Predictor

This section describes the proposed approach, which in effect comes down to the definition of a suitable nonconformity measure for multi-label classification and its use with the CCP framework. In order to use the CCP framework in the multi-label setting, the set of possible classifications {Y 1 , . . . , Y c } is replaced by the powerset P ({Y 1 , . . . , Y n }), where n is the number of the original classes of the problem. In the experiments that follow RBF Neural Networks for multilabel learning (ML-RBF) [START_REF] Zhang | ML-RBF: RBF neural networks for multi-label learning[END_REF] is used as underlying algorithm as it seems to be one of the best performing algorithms designed specifically for multi-label problems. However the proposed approach is general and can be used with any other method which gives scores for each class.

After being trained on a given training set, for every test example ML-RBF produces a score for each possible class of the task at hand. It then outputs as its prediction the labelset containing all classes with score higher than zero; higher score indicates higher chance of the class to be in the labelset. In order to use the scores produced by ML-RBF for computing the nonconformity scores of the proposed Cross-Conformal Predictor, the former were transformed to the range [0, 1] with the logistic sigmoid function

f (x) = 1 1 + e -x .
The nonconformity measure [START_REF] Tsoumakas | Mining multi-label data[END_REF] for the multi-label CCP can now be defined based on the transformed outputs of ML-RBF for (x i , ψ i ) after being trained on ∪ m =k S m , m = 1, . . . , K as

α i = n j=1 |t j i -o j i | d , ( 7
)
where n is the number of possible classes, t j i is 1 if Y j ∈ ψ i and 0 otherwise, and o j i is the transformed output of the ML-RBF corresponding to class Y j . Finally d is a parameter of the algorithm which controls the sensitivity of the nonconformity measure to small differences in o j i when |t j i -o j i | is small as opposed to when it is large. This nonconformity measure takes into account the distance of all outputs from the true values. Note that in the case of test examples t j i is the corresponding value for each assumed labelset.

The nonconformity measure [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF] takes into account only the outputs of the ML-RBF for each instance. However, it is very strange to have a pair of labels in a labelset that have never appeared together in the training set. So [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF] was extended to take into account the occurrence of each pair of labels in the training set ∪ m =k S m , m = 1, . . . , K. This extended nonconformity measure is defined as:

α i = n j=1 |t j i -o j i | d + λ 1≤j<r≤n t j i t r i µ j,r , (8)
where µ j,r is 0 if the labels Y j and Y r have been observed together in the labelset of at least one instance of the training set and 1 otherwise. In effect the additional part of this nonconformity measure adds λ to the nonconformity score of an example for each pair of labels in the labelset of the example (t j i t r i = 1) which have never been observed together in any instance of the training set (µ j,r = 1). The parameter λ adjusts the sensitivity to this part of the measure. A high value of λ makes this part of the measure dominate when a pair of labels that has never been observed before exists in the labelset of the example.

The complete proposed approach is derived by plugging in (8) as A in (5) and computing the p-value of each possible labelset with [START_REF] Vovk | Cross-conformal predictors[END_REF]. The resulting p-values can be used in either of the two ways described in Section 2.

Experiments and Results

Data Sets

To evaluate the performance of the proposed approach two data sets from different application domains were used, one from the semantic scene analysis domain and one from the bioinformatics domain. The first data set, scene [START_REF] Boutell | Learning multi-label scene classification[END_REF], is concerned with the semantic classification of pictures into one or more of the classes: beach, sunset, foliage, field, mountain and urban. It consists of 1211 training and 1196 test examples, each described by 294 features. The second data set, yeast [START_REF] Elisseeff | A kernel method for multi-labelled classification[END_REF], is concerned with predicting the functional classes of genes in the Yeast Saccharomyces cerevisiae. Each gene is described by the concatenation of microarray expression data and a phylogenetic profile, and is associated with a set of 14 functional classes. The data set contains 1500 genes as training set and 917 genes as test set, each described by 103 features. Both data sets were obtained from the website of the Mulan library [START_REF] Tsoumakas | Mining multi-label data[END_REF].

Single Prediction Evaluation

The first set of experiments evaluates the quality of the single predictions produced by the ML-RBF CCP and compares it to that of its underlying method and those of three other popular multi-label techniques, namely BP-MLL [START_REF] Zhang | Multi-label neural networks with applications to functional genomics and text categorization[END_REF], ML-kNN [START_REF] Zhang | ML-kNN: A lazy learning approach to multi-label learning[END_REF] and ML-NB [START_REF] Zhang | Feature selection for multi-label naive bayes classification[END_REF]. Four evaluation measures for multi-label classification were used. The first is Hamming Loss (HL), which is the most popular measure for multi-label problems, defined as:

HL = 1 g l+g i=l+1 |ψ i ψi | n , (9)
where {(x l+1 , ψ l+1 ), . . . , (x l+g , ψ l+g )} are the test examples, ψi is the predicted labelset for example i and is the symmetric difference between two sets. The second measure is Classification Accuracy (CA) defined as:

CA = 1 g l+g i=l+1 I(ψ i = ψi ), (10)
where I(true) = 1 and I(f alse) = 0. This measure is rather strict as it requires the predicted and true labelsets to be identical. The third and fourth measures are the macro averaged and micro averaged F-measure, which is the harmonic mean of precision and recall. The F-measure for a single label is defined as:

F (tp, tn, f p, f n) = 2tp 2tp + f p + f n , ( 11
)
where tp is the number of true positives, tn the number of true negatives, f p the number of false positives and f n the number of false negatives. In the multilabel case, if tp j , tn j , f p j and f n j are the same values for each label Y j , then its macro averaged version is defined as:

F macro = 1 n n j=1 F (tp j , tn j , f p j , f n j ). ( 12
)
The micro averaged version of the F-measure is defined as: Tables 1 and2 present the performance of the proposed approach together with that of its underlying algorithm and that of BP-MLL, ML-kNN and ML-NB on the scene and yeast data sets respectively. The best values for each measure are highlighted in bold. For ML-RBF the fraction parameter was set to 0.01 and the scaling factor to 1 as in [START_REF] Zhang | ML-RBF: RBF neural networks for multi-label learning[END_REF]. In the case of the CCP the number of folds for each data set was chosen so that each fold contained approximately 100 training instances, therefore 12 folds were used for the scene data set and 15 folds for the yeast dataset. The parameter d of the nonconformity measure (8) was set to 4, which seems to be a good choice based on the performed experiments, while for λ the two extreme values of 0 and 1 were used. Setting λ to 0 in effect corresponds to nonconformity measure [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF] and setting it to 1 makes the nonconformity score of any labelset containing a pair of classes that has never been observed in the training set always higher than others. For BP-MLL the number of hidden neurons was set to 20% of the input dimensionality, the learning rate to 0.05 and the training epochs to 100 as in [START_REF] Zhang | Multi-label neural networks with applications to functional genomics and text categorization[END_REF]. For ML-kNN the number of nearest neighbours was set to 10 and the smoothing parameter to 1 as in [START_REF] Zhang | ML-kNN: A lazy learning approach to multi-label learning[END_REF]. For ML-NB the percentage of remaining features after PCA was set to 0.3 as in [START_REF] Zhang | Feature selection for multi-label naive bayes classification[END_REF].

F micro = F   n j=1 tp j , n j=1 tn j , n j=1 f p j , n j=1 f n j   . ( 13
)
The results in these tables show that not only the proposed approach provides important additional information about the likelihood of each of its predictions being correct, but it also outperforms its underlying algorithm and the three other popular multi-label techniques. The only exception is the classification accuracy of the ML-RBF in the case of the yeast data set, but even in this case the accuracy of the proposed approach with both values of λ is very close.

Comparing the performance of the ML-RBF CCP with the two different λ values one can see that in the case of the scene data there is only a negligible difference while in the case of the yeast data the performance remains the same for all evaluation measures. Therefore the value of λ does not have any important effect on the performance of the single predictions produced by the ML-RBF CCP. It does however affect the quality of the resulting p-values as will be shown in the next subsection.

Prediction Region Evaluation

The main advantage of the proposed approach over other multi-label techniques is the production of a p-value for each possible labelset of a new unseen instance, which can be translated either to confidence and credibility measures for its prediction or to prediction sets that are guaranteed to contain the true labelset at a required confidence level. This subsection examines the informativeness and reliability of the resulting prediction sets and consequently of the computed p-values and confidence measures. More specifically given a required level of confidence 1 -δ, the ML-RBF CCP produces a set of labelsets that has at most δ chance of not containing the true labelset of the unseen instance. The informativeness of this set of labelsets can be assessed in terms of its size, while its reliability can be assessed by the percentage of cases for which it does not contain the true labelset, this percentage should be less than or very near δ. Tables 3 and4 present the results of the proposed approach in this setting for the scene and yeast data sets respectively with the nonconformity measure parameter λ set to 0 and 1. The same parameters reported in Subsection 4.2 were used. The two tables report the sizes of the prediction sets produced for the 95%, 90% and 80% confidence levels together with the observed error percentages, i.e. the percentages of prediction sets that did not contain the true labelset.

Table 3 reports the results for the scene data set in terms of the percentage of prediction sets containing only 1, 2, 3 to 4, 5 to 8 and 9 to 16 labelsets for each One can be 95% confident in about 60% of the test instances by considering less than 4 out of the possible 63 labelsets. By reducing the required confidence to 90% one can be certain in a single labelset for about 18% of the test instances and between one or two labelsets for about half the test instances. Finally with a confidence level of 80% a single labelset is given for more than 60% of the test instances. In terms of empirical reliability, only the percentages of errors for the 80% confidence level are slightly higher than the required significance level, which can be attributed to statistical fluctuations. Table 4 presents the same results for the yeast data set. In this case the prediction sets contained a much higher number of labelsets so the table reports the percentage of prediction sets containing between 2 i +1 and 2 i+1 labelsets with i = 6, . . . , 13 for each confidence level; there were no prediction sets containing less than 2 6 labelsets. The rather big size of the resulting prediction sets is not strange baring in mind the very low classification accuracy of this data set, which is only 18.65%. Comparing the results obtained with the nonconformity measure parameter λ set to 0 with those obtained with λ = 1, again shows the superiority of nonconformity measure [START_REF] Zhang | ML-RBF: RBF neural networks for multi-label learning[END_REF], as it produces smaller prediction sets. Considering the high difficulty of the particular task one can say that the resulting prediction sets are quite informative. The number of labelsets needed to satisfy the 80% confidence level is 1/16th or less (≤ 2 10 ) of all the possible labelsets for 85% of the test instances. Finally in terms of empirical reliability, the percentage of errors observed is in all cases below the required significance level.

Conclusions

This work examined the application of the conformal prediction framework to the multi-label setting. Unlike the other techniques developed for multi-label problems, the proposed approach accompanies each of its predictions with reliable measures of confidence. Experimental results on two popular multi-label data sets have shown that not only the proposed approach provides important additional information for each prediction, but it also outperforms other popular multi-label techniques. Furthermore its confidence measures have been shown to be informative and reliable. The provision of confidence measures can be very helpful in practical applications, considering the high uncertainty that exists in this setting.

Future work includes the development of additional nonconformity measures and the experimentation with more multi-label data. In addition generating separate p-values for each class and combining them for obtaining the p-value of each labelset could also be examined as an alternative. Finally the possibility of generating a ranking of the possible classes for each instance would also be a good addition to multi-label CP.

Table 1 .

 1 Performance of the proposed approach and comparison to that of other multi-label algorithms on the scene data set.

			Evaluation Measure
	Algorithm	HL	CA	Fmacro	Fmicro
	ML-RBF CCP with λ = 0	0.0928	0.6798	0.7417	0.7363
	ML-RBF CCP with λ = 1	0.0927	0.6831	0.7410	0.7358
	Orgiginal ML-RBF	0.0959	0.5468	0.6922	0.6890
	BP-MLL	0.2903	0.1630	0.0509	0.1665
	ML-kNN	0.0953	0.6012	0.7189	0.7183
	ML-NB	0.1309	0.4105	0.6230	0.6221

Table 2 .

 2 Performance of the proposed approach and comparison to that of other multi-label algorithms on the yeast data set.

			Evaluation Measure
	Algorithm	HL	CA	Fmacro	Fmicro
	ML-RBF CCP with λ = 0	0.1954	0.1821	0.3896	0.6432
	ML-RBF CCP with λ = 1	0.1954	0.1821	0.3896	0.6432
	Orgiginal ML-RBF	0.1970	0.1865	0.3891	0.6407
	BP-MLL	0.2272	0.0960	0.3047	0.6212
	ML-kNN	0.1980	0.1658	0.3567	0.6360
	ML-NB	0.2115	0.1254	0.3428	0.6152

Table 3 .

 3 Prediction set sizes and error rates at the 95%, 90% and 80% confidence levels for the scene data set.

			With λ = 0			With λ = 1
	# of	Confidence Level	Confidence Level
	labelsets	95%	90%	80%	95%	90%	80%
	1	0.00%	2.34%	54.93%	6.77%	18.23%	62.63%
	2	0.08%	25.50%	34.20%	9.62%	32.69%	29.43%
	3 to 2 2	25.50%	65.64%	10.87%	43.31%	45.07%	7.94%
	(2 2 + 1) to 2 3	72.74%	6.52%	0.00%	40.05%	4.01%	0.00%
	(2 3 + 1) to 2 4	1.67%	0.00%	0.00%	0.25%	0.00%	0.00%
	Errors	3.76%	9.28%	21.07%	3.60%	9.28%	20.99%

Table 4 .

 4 Prediction set sizes and error rates at the 95%, 90% and 80% confidence levels for the yeast data set. there was no prediction set containing more than 16 labelsets out of the possible 63. The last row of the table reports the percentage of errors observed for each confidence level. Comparing the results obtained with the nonconformity measure parameter λ set to 0 with those obtained with λ = 1, one can see that the latter produces much more informative prediction sets. Taking into account the classification accuracy of this data set (68.31%) and the large number of possible labelsets, the resulting prediction sets are quite tight.

			With λ = 0			With λ = 1
	# of	Confidence Level	Confidence Level
	labelsets	95%	90%	80%	95%	90%	80%
	(2 6 + 1) to 2 7	0.00%	0.00%	1.42%	0.00%	0.00%	1.42%
	2 7 < l ≤ 2 8	0.00%	0.11%	3.71%	0.00%	0.22%	4.58%
	2 8 < l ≤ 2 9	0.55%	3.82%	14.07%	0.76%	4.36%	18.65%
	2 9 < l ≤ 2 10	3.05%	7.09%	60.96%	3.93%	10.14%	60.32%
	2 10 < l ≤ 2 11	8.94%	34.46%	19.85%	13.85%	45.58%	15.05%
	2 11 < l ≤ 2 12	38.71%	52.89%	0.00%	57.58%	39.48%	0.00%
	2 12 < l ≤ 2 13	48.64%	1.64%	0.00%	23.88%	0.22%	0.00%
	2 13 < l ≤ 2 14	0.11%	0.00%	0.00%	0.00%	0.00%	0.00%
	Errors	4.69%	9.38%	19.85%	4.80%	9.60%	19.96%
	confidence level;


Chenzhe Zhou

Ilia Nouretdinov

Zhiyuan Luo

Alex Gammerman

SVM Venn Machine with k-Means Clustering

Keywords: Venn Machine, Support Vector Machine, k-means clustering

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

Introduction

Classification is one of the major tasks in machine learning. It gives predictions for the new objects based on known properties learned from the training data set. However, most algorithms could only give single prediction (i.e. label). Demand of probabilistic prediction has arisen in view of the fact that sometimes we appreciate probabilities more than single predictions. A simple example is the probabilistic weather forecasting.

But in some area, single probabilistic prediction has not yet been enough. The term multi-probabilities is then brought to mind, namely, we announce several probability distributions for the new label rather than a solitary one. Venn predictor (or Venn Machine) is one of the multi-probabilistic classification systems [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF]. There are many Venn predictors, each taxonomy used in the algorithm defines a Venn predictor even if the underlying algorithms are the same.

In our previous paper [START_REF] Zhou | A comparison of venn machine with platt's method in probabilistic outputs[END_REF], we introduced a Venn predictor with Support Vector Machines (SVM) as its underlying algorithm, which converts numerical predictions of SVM into a taxonomy. That approach was applicable to any method that initially supplied predictions with prediction scores such as the distance to the hyperplane in SVM. Nonetheless, the process is very simple: all available scores are firstly sorted and then divided into several groups by equallength intervals according to which interval the score lies. Each of these groups is a category. However, that approach could only be applied in binary cases. In this paper, we propose a method to generalize binary Venn Machine with SVM to a method capable for multi-class cases. Then we consider two alternative methods that may be more accurate: SVM Venn Machine with k-means clustering and Venn-ABERS predictor. These two algorithms are also applicable to any machine learning algorithms with prediction scores.

Methodology

In this section, two kinds of Venn predictors that use SVM as their underlying algorithm will be introduced together with our alternative methods. They are SVM Venn Machine with homogeneous intervals (VM-SVM-HI) generalized from the binary-only version of [START_REF] Zhou | A comparison of venn machine with platt's method in probabilistic outputs[END_REF] together with our alternative method SVM Venn Machine with k-means clustering (VM-SVM-KM) and the Venn-ABERS predictor based on SVM (VA-SVM) proposed by Vladimir Vovk [START_REF] Vovk | Venn predictors and isotonic regression[END_REF]. The former two algorithms could be implemented in both multi-class cases and binary cases, while VA-SVM could only deal with binary data sets.

Venn Machine

Venn Machine is a multi-probabilistic predictor described in [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF]. The basic idea of Venn Machine is to divide every example into its corresponding category based on certain rules and then the frequencies of labels in the chosen category are used as probabilities for the new object's label. Taxonomy is the way how the examples are divided into categories. The underlying algorithm is the algorithm used in the taxonomy.

Assuming a standard machine learning classification problem: given a training set of examples z 1 , z 2 , . . . , z n-1 . Each z i consists of a pair of object x i and label y i . The possible labels y i (y i ∈ Y) are finite. And we are also given a test object x n . Our task is to predict the label y n for the new object x n and give the estimation of the likelihood that our prediction is correct.

Supposing we have a taxonomy A n , consider a label y ∈ Y for the new object x n . A n assigns a category τ i to an example z i

τ i = A n ( z 1 , . . . , z i-1 , z i+1 , . . . , z n , z i ) ( 1
)
where n is the number of objects in the bag, τ i ∈ T is one of the finite categories and z i is the pair (x i , y i ), z n is the pair (x n , y). Moreover, we assign z i and z j to the same category if and only if

A n ( z 1 , . . . , z i-1 , z i+1 , . . . , z n , z i ) = A n ( z 1 , . . . , z j-1 , z j+1 , . . . , z n , z j ) (2)
The category τ n contains z n = (x n , y). Let p y be the empirical probability distribution of the labels in category τ n .

p y (y ) := |{(x * , y * ) ∈ τ n : y * = y }| |τ n | (3)
p y is a probability distribution on Y.

Having tried every possible label for x n , we get a Venn predictor. The predictor P n := {p y : y ∈ Y} is a multi-probabilistic predictor consists of K distributions, where K = |Y|. Then we could calculate a K × K frequency matrix P . The quality of a column is the minimum entry of the column. Let the best column which has the highest quality be j best . Then our predicted label is j best and the interval of possibility that our prediction is correct is

[ min i=1,...,K P i,j best , max i=1,...,K P i,j best ] (4)
Underlying Algorithm for Taxonomy Any algorithm that generates or predicts a numeric score for the example could be implemented in our taxonomy. However, we mainly focus on Venn predictors with SVM as the underlying algorithm in this paper. The decision function in SVM is a kind of scoring functions. Therefore, we use the values derived from the decision function of SVM (i.e. the values prior to applying a sign function) as part of our design.

Homogenous Intervals One of the simplest ways to design taxonomies is stated as follows. Firstly we use the training set to train an SVM and calculate the decision values (d(x) = w, x i +d) of all examples in the training set and the new object. Secondly the whole range of decision values obtained will be divided into several intervals of equal length. Each interval is a category and objects of which the decision values fall into the same interval are of the same category. This design was introduced in [START_REF] Zhou | A comparison of venn machine with platt's method in probabilistic outputs[END_REF] and could only used in binary case. Now we will discuss the generalization and alternative to it.

Combined Decision Function

In multi-class cases, we will have several binary SVM classifiers regardless of whether One-vs-One or One-vs-All approach is used. A scheme for multi-class SVM using One-vs-All approach was developed by Lambrou et al. in [START_REF] Lambrou | Reliable probability estimates based on support vector machines for large multiclass datasets[END_REF], which uses the largest decision value as the score. Generally, One-vs-One SVM is more efficient in accuracy than One-vs-All SVM. Therefore, we need to develop a new function to combine the outputs of all Onevs-One SVM classifiers and transform them into a single prediction score which could be used by Venn Machine. We call such function a Combined Decision Function.

For a data set with k possible labels: {0, 1, . . . , k -1}, there are k(k -1)/2 binary SVM if we use One-vs-One approach. For each possible label, there are k -1 related SVM decision functions. Then we use [START_REF] Lambrou | Reliable probability estimates based on support vector machines for large multiclass datasets[END_REF] to calculate the combined decision function D(x) for the new example x,

D(x) = ŷ + 1 k -1 k-1 i=0,i =ŷ N (f ŷi (x)) ( 5
)
where ŷ is the overall predicted label done by max-wins voting strategy in Onevs-One SVM, f ŷi (x) is the decision function of SVM classifier on ŷ-vs-i, N is a function that does the normalized transformation to [0, 1]. Another point we need to declare here is that in f ŷi (x) we always put ŷ before i which means we need to apply an opposite operation when ŷ is greater than i. Since the examples of label ŷ are treated as negative examples in i-vs-ŷ classifier of a binary SVM. This function firstly selects all k -1 related SVM and applies an opposite operation if ŷ is not treated as the positive class in the binary SVM classifier. Then it does the normalisation to transform the values into [0, 1]. Finally, we output the arithmetic mean of them added with ŷ as the combined decision value of new example x. The reason that adding ŷ to the arithmetic mean is that it could prevent the decision values of different classes stack at the same area.

Dividing Intervals by k -Means Clustering Instead of dividing the intervals homogeneously, we came up with a new dividing scheme, which uses k-means clustering [START_REF] Forgy | Cluster analysis of multivariate data: Efficiency vs interpretability of classifications[END_REF][START_REF] Lloyd | Least squares quantization in pcm[END_REF] to divide all decision values.

k-means clustering is a cluster analysis method which aims to divide n objects into k clusters in which each object belongs to the cluster with the nearest mean. Given a set of objects (x 1 , x 2 , . . . , x n ), where each object x i ∈ R d is a d-dimensional real vector, k-means clustering aims to partition the n objects into k sets (k ≤ n) S = {S 1 , S 2 , . . . , S k } so as to minimise the within-cluster sum of squares (WCSS):

arg min S k i=1 xj ∈Si x j -µ i 2 (6)
where µ i is the mean value of points in S i .

In our design, dimension d is fixed to "1", while the number of clusters is equal to the number of possible labels. So the heuristic algorithm we used could be described as below.

1. k initial means values are randomly generated within the data domain. 2. k clusters are created (or reassigned) by associating every object with the nearest mean value. 3. The centroid of each of the k clusters becomes the new mean value. 4. Steps 2 and 3 are repeated until the change of WCSS [START_REF] Lloyd | Least squares quantization in pcm[END_REF] between two states declines to be less than = 10 -4 .

Having applied k-means clustering, we divided the decision values into categories which could be used to calculate the matrix for new examples and make the probabilistic predictions as the standard Venn Machine does.

Venn-ABERS Predictor

Venn-ABERS predictor is a recently developed algorithm for multi-probabilistic prediction. It is modified from Zadrozny and Elkan's procedure of probability forecasting [START_REF] Zadrozny | Transforming classifier scores into accurate multiclass probability estimates[END_REF], which cannot be well calibrated. The modification introduced Venn predictors into the procedure to overcome the problem of potentially weak calibration as a result of the fact that Venn predictors are always well calibrated and guaranteed to be well calibrated under the exchangeability assumption. The basic idea of pre-trained Venn-ABERS predictor is that the training set is split into two parts: the proper training set and the calibration set. The proper training set is used to train the learning machine and predict the label for new examples, while the calibration set is used to calculate the probabilistic outputs for the predicted labels. The calibration set will be turned into a monotonically increasing set in this algorithm according to [START_REF] Ayer | An empirical distribution function for sampling with incomplete information[END_REF].

Before we discuss Venn-ABERS predictor, there are some notions to be introduced yet. First notion is the term "scoring algorithm". Scoring algorithm is an algorithm that trains a classifier on the training set and uses the classifier to output a prediction score s(x) for the new example x and predicts the label of x to be "1" if and only if s(x) ≥ c (c is a fixed threshold). So s is hereby called the scoring function. Many machine learning algorithms for classification are scoring algorithms. In our case, as what SVM defines, the decision function of SVM is a scoring function, since we assign a new example the positive label "+1" if and only if its decision value is greater than zero and vice versa for the negative label. The second notion is "isotonic calibrator ", which is a monotonically increasing function on the set {s(x 1 ), . . . , s(x l )} that maximizes the likelihood

l i=1 p i , where p i := g(s(x i )) if y i = 1 1 -g(s(x i )) if y i = 0 (7)
this function g is unique and can be found by using the "pool-adjacent violators algorithm" (PAVA) introduced in [START_REF] Ayer | An empirical distribution function for sampling with incomplete information[END_REF].

The workflow of Venn-ABERS predictor is as follows. Assuming a standard binary machine learning problem: a training set of examples z 1 , z 2 , . . . , z l . Each z i consists of a pair of object x i , and label y i . The possible labels are binary, that is, y ∈ {0, 1}. And we are also given a test object x. Our task is to predict the label y for the new object x and give the estimation of the likelihood that our prediction is correct.

Let us split the training set z 1 , z 2 , . . . , z l into two parts: the proper training set z 1 , z 2 , . . . , z m of size m (m < l) and the calibration set z m+1 , z m+2 , . . . , z l . And s : X → R is the scoring function of training set z 1 , z 2 , . . . , z m . Given a new example x, we have two calibrators. Let g 0 be the isotonic calibrator for (s(x m+1 ), y m+1 ), (s(x m+2 ), y m+2 ), . . . , (s(x l ), y l ), (s(x), 0) , g 1 be the calibrator for (s(x m+1 ), y m+1 ), (s(x m+2 ), y m+2 ), . . . , (s(x l ), y l ), (s(x), 1) .

To achieve the isotonic calibrator, we do the followings according to the definition of PAVA. First we arrange the pairs (s(x i ), y i ) in the increasing order according to the values of score function s(x i ). Having obtained a binary sequence consisting of labels y i , we applied PAVA to find the increasing sequence of them. The final isotonic calibrator g is a function mapping the increasing scores to the increasing sequence (i.e. probabilities). As the score increases, the object is more likely to be "1" in correlation with the increasing sequence.

Then the multi-probability prediction outputs for that the predicted label should be "1" is {p 0 , p 1 }, where p 0 := g 0 (s(x)) and p 1 := g 1 (s(x)). And for the reason that we need to predict the probability for the prediction label is correct, we should transform the bounds {p 0 , p 1 } to {1 -p 1 , 1 -p 0 } when the predicted label is "0".

Experimental Results

To compare our algorithm to SVM Venn Machine with homogeneous intervals and Venn-ABERS predictor, we used eight data sets from the real world which could be easily obtained from UCI Repository (http://archive.ics.uci.edu/ ml/) except that SVMguide1 is obtained from the website of LibSVM [START_REF] Chang | Libsvm: a library for support vector machines[END_REF]. The data sets we used in this paper could be divided into two parts based on their number of classes. The details of these data sets are summarised in Table 1.

Experimental Settings

For VM-SVM-KM, the number of clusters and the initial means, which are the two key features of k-means clustering, are often regarded as its biggest drawbacks. The number of clusters is an input parameter: an inappropriate choice of k may yield poor results. That is why, when performing k-means clustering, it is important to run diagnostic checks for determining the number of clusters in the data set. The choice of initial means might lead the convergence to a local minimum which may produce counterintuitive results. A good design of a combined decision function could make it easier to avoid these two drawbacks.

To have a more intuitive view of our combined decision function described in [START_REF] Lambrou | Reliable probability estimates based on support vector machines for large multiclass datasets[END_REF], we applied the algorithm to Satimage data set and plotted the histogram in Fig. 1, roughly representing the distribution of the decision values.

It can be seen obviously from the figure that there were 6 clusters in the data set, the exact number of the possible labels. The reason for this is that the decision function spreads out the values into (0, k) by adding the most possible labels. Furthermore, each cluster i (i = 1, 2, . . . , k) is approximately within the range of (i -1, i), which means we could choose the initial means from each range to avoid the local minimum trap as much as possible and speed up the convergence process. We conducted k-means clustering to these decision values and calculated the 6 centroids: 0.63, 1.91, 2.64, 3.33, 4.57, 5.59. The result seems to be a reasonable reflection of the histogram.

Then we could come to our decision that we set the number of clusters the same as the number of possible labels and we choose the initial means as 0.5, 1.5, . . . , k -0.5 if the possible labels are 0, 1, . . . , k -1.

Additionally, we need to notice that k-means clustering uses Euclidean distance as a metric and variance as a measure of cluster scatter, which makes it tend to produce equal-sized clusters. Since data is split halfway between cluster means, this can lead to suboptimal splits as some objects will be attributed to the incorrect cluster, especially for unbalanced data set as Satimage data set.

Except all the settings for the underlying algorithm, we need another setting for VA-SVM. It is the size of the calibration set. Having given careful consideration to both accuracy and narrowness of the bounds, we decided to take 30% of the whole data set as the calibration set, And the calibration set was stratified selected from the whole training set, which means the distribution of classes in the calibration set was the same as in the training set.

Although the size of proper training set in VA-SVM is smaller comparing to the size of training set in our algorithms, this is still a fair comparison because we use the same original training set for all algorithms, otherwise VA-SVM will need extra examples for probabilistic predictions. We also noticed that Venn-ABERS predictor is an inductive Venn predictor while Venn Machine is a transductive Venn predictor. The gap between inductive and transductive learning algorithms are not distinguishable in our offline setting. Because in offline setting, we use the fixed predictors to make predictions for testing set. Furthermore, in VA-SVM we repeat the computations of isotonic calibrators for each testing object which still involve all examples in calibration set.

Comparisons and Results

For binary cases, we applied VM-SVM-KM, VM-SVM-HI and VA-SVM to the data sets in the offline setting. While for multi-class cases, we only applied VM-SVM-KM and VM-SVM-HI to the data sets in both offline setting and online setting. Hence, there were three comparisons described as below. All the SVMs in these algorithms were using RBF kernel. Additionally, the parameters of SVM for each data set, including cost C and σ in RBF kernel, were determined by grid search on the training set and retained the same over corresponding algorithms respectively. The algorithms were compared in terms of their accuracies and probabilistic outputs in these data sets. In addition, we calculated the Brier scores (introduced in [START_REF] Brier | Verification of forecasts expressed in terms of probability[END_REF]) of the mean of the probabilistic bounds as evaluation for binary data sets.

The experimental results of VM-SVM-KM compared with VM-SVM-HI and VA-SVM are shown in Table 2. The comparison results of our algorithm against VM-SVM-HI for all multiclass data sets in the offline setting are shown in Table 3. And the results for the online setting are shown in Table 4.

In order to giving a more intuitive comparison, we also give figures on online performance for Wine data set in Fig. 2.

Discussion and Conclusion

From the results shown in Table 2 which comparing our method to VM-SVM-HI and VA-SVM, we could draw the following conclusions. First, VM-SVM-KM performed better in accuracy among these three data sets nevertheless the increases were small. Furthermore, VA-SVM used 30% of the training set as the calibration set which did not participate in the training of classifiers; hence it may lead to worse results. Second, the accuracies of VA-SVM slightly outnumbered the upper bound in WBC and Splice data sets, which could be due to the offline setting. Third, the probability bounds of VA-SVM were the narrowest while VM-SVM-HI had the widest bounds and VM-SVM-KM was in-between. This is the advantage of VA-SVM in view of our preference for narrow bounds. It is also backed by the Brier scores results: VA-SVM and VM-SVM-KM had close Brier scores while VM-SVM-HI had the worst results.

Another point is that VA-SVM does not calibrate their predicted label according to the probability, more specifically it is an algorithm that generates the probabilities from the scores only, while our algorithm gives predictions based on the highest likelihood. Except the improvement in accuracy, VM-SVM-KM is easy to configure because the number of clusters is the only input parameter of this algorithm which is equal to the number of classes.

From the results presented in Table 3 and Table 4 where the performance of VM-SVM-KM is compared with VM-SVM-HI in both offline and online setting, we can discover the following points.

First, it can be observed that all accuracies were within the probabilistic outputs in the offline setting, while in the online setting the accuracies exceeded the bounds in Segment and Vehicle data sets. Second, after implementing the k-means clustering, the accuracies are improved in both settings. However, in the offline setting, the improvements ranged from 0.8% to 3.8% depending on the data sets. In the online setting, the difference between these two algorithms became smaller, only 0.6% to 2.5%. Third, probability bounds become narrower after applying the k-means clustering, mostly benefiting from the rise of lower bounds. An intuitive comparison is shown in Fig. 2. The cumulative errors and cumulative error bounds in the figures all decreased after implementing k-means clustering, and the bounds became narrower in the meantime.

In summary, the improvement in each of the eight data sets was not significant which is due to the consistency of these data sets. Nevertheless, we still believe that SVM Venn Machine with k-means clustering is better when compared with homogeneous intervals since it could yield better accuracy and narrower bounds. However, in comparison with Venn-ABERS predictor, our algorithm is good on accuracy and weak on narrowness of the bounds. Despite that, our algorithm is easier to set up, and it predicts the most likely label while Venn-ABERS predictor only generates the probabilities.

Fig. 1 .

 1 Fig. 1. Histogram of combined decision values for the Satimage data set.

Fig. 2 .

 2 Fig. 2. Comparison of online performances for the Wine data set.

Table 1 .

 1 Main characteristics for each data set.

	Data Set	# of	# of	# of	Training	Testing
		Objects	Features	Classes	Set Size	Set Size
	WBC	683	10	2	400	283
	SVMguide1	7089	4	2	3089	4000
	Splice	3175	60	2	1000	2175
	Satimage	6435	36	6	4435	2000
	Segment	2310	19	7	1500	810
	DNA	3186	180	3	2000	1186
	Wine	178	13	3	100	78
	Vehicle	846	18	4	500	346

Table 2 .

 2 The offline accuracy and probability results on the binary data set.

	Data Set	Taxonomy	Accuracy	Prob. Outputs	Brier Score
	WBC	VM-SVM-KM	97.53%	[86.34%,98.94%]	0.0325
		VM-SVM-HI	97.22%	[83.63%,98.70%]	0.0369
		VA-SVM	97.17%	[85.67%,95.97%]	0.0315
	SVMguide1	VM-SVM-KM	96.93%	[91.27%,98.42%]	0.0362
		VM-SVM-HI	95.79%	[89.59%,98.97%]	0.0406
		VA-SVM	95.95%	[93.67%,96.29%]	0.0370
	Splice	VM-SVM-KM	90.21%	[82.44%,96.07%]	0.0884
		VM-SVM-HI	89.52%	[80.15%,97.35%]	0.0939
		VA-SVM	89.15%	[83.40%,88.32%]	0.0878

Table 3 .

 3 The offline accuracy and probability results on the multi-class data set.

	Data Set	Algorithm	Accuracy	Probabilistic Outputs
	Satimage	VM-SVM-HI	84.18%	[75.48%,96.92%]
		VM-SVM-KM	86.56%	[81.18%,93.33%]
	Segment	VM-SVM-HI	90.88%	[74.61%,96.68%]
		VM-SVM-KM	91.65%	[75.64%,95.60%]
	DNA	VM-SVM-HI	94.34%	[81.39%,98.07%]
		VM-SVM-KM	96.65%	[87.25%,99.48%]
	Wine	VM-SVM-HI	92.30%	[81.19%,97.11%]
		VM-SVM-KM	96.11%	[86.53%,98.47%]
	Vehicle	VM-SVM-HI	67.63%	[56.42%,77.48%]
		VM-SVM-KM	69.15%	[60.65%,79.20%]

Table 4 .

 4 The online accuracy and probability results on the multi-class data set.

	Data Set	Algorithm	Accuracy	Probabilistic Outputs
	Satimage	VM-SVM-HI	80.94%	[80.20%,81.69%]
		VM-SVM-KM	83.40%	[83.24%,83.86%]
	Segment	VM-SVM-HI	88.40%	[88.92%,93.20%]
		VM-SVM-KM	89.96%	[90.11%,91.50%]
	DNA	VM-SVM-HI	89.12%	[88.46%,89.86%]
		VM-SVM-KM	89.70%	[89.25%,90.48%]
	Wine	VM-SVM-HI	91.53%	[87.57%,94.35%]
		VM-SVM-KM	93.22%	[91.67%,96.87%]
	Vehicle	VM-SVM-HI	66.04%	[70.24%,72.17%]
		VM-SVM-KM	67.83%	[69.48%,71.02%]