Keywords: Inférence bayésienne, classification, lois conditionnelles compatibles, modèle de Boltzmann, constante de normalisation, pseudo-vraisemblance, méthode d'Ogata, simulation parfaite par couplage, algorithmes MCMC

 dans un schéma MCMC et proposons une méthode MCMC exacte basée sur la simulation parfaite par couplage. Lorsque l'on ne peut pas utiliser la technique de la simulation parfaite ou si celle-ci s'avère trop coûteuse, nous proposons de la remplacer par la méthode de l'échantillonnage de Gibbs. Nous illustrons les performances de cet algorithme sur divers jeu de données.

Introduction 1.Deterministic versus statistical classification

Supervised classification has long been used in both Machine Learning and Statistics to infer about the functional connection between a group of covariates (or explanatory variables) and a vector of indicators (or classes) (see, e.g., [START_REF] Mclachlan | Discriminant Analysis and Statistical Pattern Recognition[END_REF][START_REF] Ripley | Neural networks and related methods for classification (with discussion)[END_REF][START_REF] Ripley | Pattern Recognition and Neural Networks[END_REF][START_REF] Devroye | A Probabilistic Theory of Pattern Recognition[END_REF][START_REF] Hastie | The Elements of Statistical Learning[END_REF]. For instance, the method of boosting [START_REF] Freund | A decision-theoretic generalization of on-line learning and an application to boosting[END_REF] has been developed for this very purpose by the Machine Learning community and has also been assessed and extended by statisticians [START_REF] Hastie | The Elements of Statistical Learning[END_REF][START_REF] Bühlmann | Bagging, boosting and ensemble methods[END_REF]Yu, 2002, 2003;[START_REF] Bühlmann | Bagging, boosting and ensemble methods[END_REF][START_REF] Zhang | Boosting with early stopping: convergence and consistency[END_REF].

The k-nearest-neighbour method is a well-established and straightforward technique in this area with both a long past and a fairly resilient resistance to change [START_REF] Ripley | Neural networks and related methods for classification (with discussion)[END_REF][START_REF] Ripley | Pattern Recognition and Neural Networks[END_REF]. Nonetheless, while providing an instrument for classifying points into two or more classes, it lacks a corresponding assessment of its classification error. While alternative techniques like boosting offer this assessment, it is obviously of interest to provide the original k-nearest-neighbour method with this additional feature. In contrast, statistical classification methods that are based on a model such a mixture of distributions do provide an assessment of error along with the most likely classification. This more global perspective thus requires the technique to be embedded within a probabilistic framework in order to give a proper meaning to the notion of classification error. [START_REF] Holmes | A probabilistic nearest neighbour method for statistical pattern recognition[END_REF] propose a Bayesian analysis of the k-nearest-neighbour-method based on these premises, and we refer the reader to this paper for background and references. In a separate paper, [START_REF] Holmes | Likelihood inference in nearest-neighbour classification models[END_REF] defined another model based on autologistic representations and conducted a likelihood analysis of this model, in particular for selecting the value of k. While we also adopt a Bayesian approach, our paper differs from [START_REF] Holmes | A probabilistic nearest neighbour method for statistical pattern recognition[END_REF] in two important respects: first, we define a global probabilistic model that encapsulates the k-nearest-neighbour method, rather than working with incompatible conditional distributions, and, second, we derive a fully operational simulation technique adapted to our model and based either on perfect sampling or on a Gibbs sampling approximation, that allows for a reassessment of the pseudo-likelihood approximation often used in those settings.

The original k-nearest-neighbour method

Given a training set of individuals allocated each to one of G classes, the classical k-nearest-neighbour procedure is a method that allocates new individuals to the most common class in their neighbourhood among the training set, the neighbourhood being defined in terms of the covariates. More formally, based on a training dataset ((y i , x i )) n i=1 , where y i ∈ {1, . . . , G} denotes the class label of the ith point and x i ∈ R p is a vector of covariates, an unobserved class y n+1 associated with a new set of covariates x n+1 is estimated by the most common class among the k nearest neighbours of x n+1 in the training set (x i ) n i=1 . The neighbourhood is defined in the space of the covariates x i , namely

N k n+1 = 1 ≤ i ≤ n; d(x i , x n+1 ) ≤ d(•, x n+1 ) (k) ,
where d(•, x n+1 ) denotes the vector of distances to x n+1 and d(•, x n+1 ) (k) denotes the kth order statistic. The original k-nearest-neighbour method usually uses the Euclidean norm, even though the Mahalanobis distance would be more appropriate in that it rescales the covariates. Whenever ties occur, they are resolved by decreasing the number k of neighbours until the problem disappears. When some covariates are categorical, other types of distance can be used instead, as in the R package knncat of [START_REF] Buttrey | Nearest-neighbor classification with categorical variables[END_REF].

As such, and as also noted in [START_REF] Holmes | A probabilistic nearest neighbour method for statistical pattern recognition[END_REF], the method is both deterministic, given the training dataset, and not parameterised, even though the choice of k is both non-trivial and relevant to the performance of the method. Usually, k is selected via cross-validation, as the number of neighbours that minimises the cross-validation error rate. In contrast to cluster-analysis set-ups, the number G of classes in the k-nearest-neighbour procedure is fixed and given by the training set: the introduction of additional classes that are not observed in the training set has no effect on the future allocations.

To illustrate the original method and to compare it later with our approach, we use throughout a toy benchmark dataset taken from [START_REF] Ripley | Neural networks and related methods for classification (with discussion)[END_REF]. This dataset corresponds to a two-class classification problem in which each (sub)population of covariates is simulated from a bivariate normal distribution, both populations being of equal sizes. A sample of n = 250 individuals is used as the training set and the model is tested on a second group of m = 1, 000 points acting as a test dataset. Figure 1 presents the dataset 1 and Table 1 displays q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -1.0 -0.5 0.0 0.5 1.0 -0.2 0.2 0.6 1.0 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -1.0 -0.5 0.0 0.5 1.0 -0.2 0.2 0.6 1.0 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q

Goal and plan

As presented above, the k-nearest-neighbour method is merely an allocation technique that does not account for uncertainty. In order to add this feature, we need to introduce a probabilistic framework that relates the class label y i to both the covariates x i and the class labels of the neighbours of x i . Not only does this perspective provide more information about the variability of the classification, when compared with the point estimate given by the original method, but it also takes advantage of the full (Bayesian) inferential machinery to introduce parameters that measure the strength of the influence of the neighbours, and to analyse the role of the variables, of the metric used, of the number k of neighbours, and of the number of classes towards achieving higher efficiency. Once again, this statistical viewpoint was previously adopted by Holmes andAdams (2002, 2003) and we follow suit in this paper, with a modification of their original model geared towards a coherent probabilistic model, while providing new developments in computational model estimation.

In order to illustrate the appeal of adopting a probabilistic perspective, we provide in Figure 3 two graphs that are by-products of our Bayesian analysis. For Ripley's dataset, the first graph (on the left) gives the level sets of the predictive probabilities to be in the black class, while the second graph (on the right) partitions the square into three zones, namely sure allocation to the red class, sure allocation to the black class and an uncertainty zone. Those three sets are obtained by first computing 95% credible intervals for the predictive probabilities and then checking those intervals against the borderline value 0.5. If the interval contains 0.5, the point is ranked as uncertain. The paper is organised as follows. We establish the validity of the new probabilistic k-nearest-neighbour model in Section 2, pointing out the deficiencies of the models advanced by Holmes andAdams (2002, 2003), and then cover the different aspects of running Bayesian inference in this k-nearest-neighbour model in Section 3, addressing in particular the specific issue of evaluating the normalising constant of the probabilistic k-nearestneighbour model that is necessary for inferring about k and an additional parameter. We take advantage of an exact MCMC approach proposed in Section 3.4 to evaluate the limitations of the pseudo-likelihood alternative in Section 3.5 and illustrate the method on several benchmark datasets in Section 4.

2 The probabilistic k-nearest-neighbour model

Markov random field modelling

In order to build a probabilistic structure that reproduces the features of the original k-nearest-neighbour procedure and then to estimate its unknown parameters, we first need to define a joint distribution of the labels y i conditional on the covariates x i , for the training dataset. A natural approach is to take advantage of the spatial structure of the problem and to use a Markov random field model. Although we will show below that this is not possible within a coherent probabilistic setting, we could thus assume that the full conditional distribution of y i given y -i = (y 1 , . . . , y i-1 , y i+1 , . . . , y n ) and the x i 's only depends on the k nearest neighbours of x i in the training set. The parameterised structure of this conditional distribution is obviously open but we opt for the most standard choice, namely, like the Potts model, a Boltzmann distribution [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF] with potential function

∼ k i δ yi (y ) ,
where ∼ k i means that the summation is taken over the observations x belonging to the k nearest neighbours of x i , and δ a (b) denotes the Dirac function. This function actually gives the number of points from the same class y i as the point x i that are among the k nearest neighbours of x i . As in [START_REF] Holmes | Likelihood inference in nearest-neighbour classification models[END_REF], the expression for the full conditional is thus

f (y i |y -i , X, β, k) = exp β ∼ k i δ yi (y ) k G g=1 exp β ∼ k i δ g (y ) k (1)
where β > 0 and X is the (p, n) matrix {x 1 , . . . , x n } of coordinates for the training set.

In this parameterised model, β is a quantity that is obviously missing from the original k-nearest-neighbour procedure. It is only relevant from a statistical point of view as a degree of uncertainty: β = 0 corresponds to a uniform distribution over all classes, meaning independence from the neighbours, while β = +∞ leads to a point mass distribution at the prevalent class, corresponding to extreme dependence. The introduction of the scale parameter k in the denominator is useful in making β dimensionless.

There is, however, a difficulty with this expression in that, for almost all datasets X, there does not exist a joint probability distribution on y = (y 1 , . . . , y n ) with full conditionals equal to (1). This happens because the k-nearest-neighbour system is usually asymmetric: when x i is one of the k nearest neighbours of x j , x j is not necessarily one of the k nearest neighbours of x i . Therefore, the pseudo-conditional distribution (1) will not depend on x j while the equivalent for x j does depend on x i : this is obviously impossible in a coherent probabilistic framework [START_REF] Besag | Spatial interaction and the statistical analysis of lattice systems (with discussion)[END_REF][START_REF] Cressie | Statistics for Spatial Data[END_REF] One way of overcoming this fundamental difficulty is to follow [START_REF] Holmes | A probabilistic nearest neighbour method for statistical pattern recognition[END_REF] and to define directly the joint distribution

f (y|X, β, k) = n i=1 exp β ∼ k i δ yi (y ) k G g=1 exp β ∼ k i δ g (y ) k .
(2)

Unfortunately, there are drawbacks to this approach, in that, first, the function (2) is not properly normalised (a fact overlooked by [START_REF] Holmes | A probabilistic nearest neighbour method for statistical pattern recognition[END_REF], and the necessary normalising constant is intractable. Second, the full conditional distributions corresponding to this joint distribution are not given by (1). The first drawback is a common occurrence with Boltzmann models and we will deal with this difficulty in detail in Section 3. At this stage, let us point out that the most standard approach to this problem is to use pseudo-likelihood following [START_REF] Besag | Bayesian image restoration, with two applications in spatial statistics[END_REF], as in [START_REF] Heikkinen | Fully Bayesian approach to image restoration with an application in biogeography[END_REF] and [START_REF] Hoeting | Bayesian model averaging: A tutorial (with discussion)[END_REF], but we will show in Section 3.5 that this approximation can give poor results. (See, e.g., [START_REF] Friel | Bayesian inference in hidden markov random fields for binary data defined on large lattices[END_REF] for a discussion of this point.) The second and more specific drawback implies that (2) cannot be treated as a pseudo-likelihood [START_REF] Besag | Spatial interaction and the statistical analysis of lattice systems (with discussion)[END_REF][START_REF] Besag | Bayesian image restoration, with two applications in spatial statistics[END_REF]since, as stated above, the conditional distribution (1) cannot be associated with any joint distribution. That (2) misses a normalising constant can be seen from the special case in which n = 2, y = (y 1 , y 2 ) and G = 2, since

2 y1=1 2 y2=1 2 i=1 exp β ∼ k i δ yi (y ) k 2 g=1 exp β ∼ k i δ g (y ) k = 2 y1=1 2 y2=1 exp (β [δ y1 (y 2 ) + δ y2 (y 1 )] /k) 1 + e β/k 2 = 2 1 + e 2β/k 1 + e β/k 2 ,
which is clearly different from 1 and, more importantly, depends on both β and k. We note that the debate about whether or not one should use a proper joint distribution is reminiscent of the opposition between Gaussian conditional autoregressions (CAR) and Gaussian intrinsic autoregressions in [START_REF] Besag | On conditional and intrinsic autoregressions[END_REF], the latter not being associated with any joint distribution.

A symmetrised Boltzmann modelling

Given these difficulties, we therefore adopt a different strategy and define a joint model on the training set as

f (y|X, β, k) = exp β n i=1 ∼ k i δ yi (y ) k Z(β, k) , (3)
where Z(β, k) is the normalising constant of the distribution. The motivation for this modelling is that the full conditional distributions corresponding to (3) can be obtained as

f (y i |y -i , X, β, k) ∝ exp β/k ∼ k i δ yi (y ) + i∼ k δ y (y i ) , (4)
where i ∼ k means that the summation is taken over the observations x for which x i is a k-nearest neighbour.

Obviously, these conditional distributions differ from (1) if only because of the impossibility result mentioned above. The additional term in the potential function corresponds to the observations that are not among the nearest neighbours of x i but for which x i is a nearest neighbour. In this model, compared with single neighbours, mutual neighbours are given a double weight. This feature is of importance in that this coherent model defines a new classification criterion that can be treated as a competitor of the standard k-nearestneighbour objective function. Note also that the original full conditional (1) is recovered as (4) when the system of neighbours is perfectly symmetric (up to a factor 2). Once again, the normalising constant Z(β, k) is intractable, except for the most trivial cases.

In the case of unbalanced sampling, that is, if the marginal probabilities p 1 = P(y = 1), . . . , p G = P(y = G) are known and are different from the sampling probabilities p1 = n 1 /n, . . . , pG = n G /n, where n g is the number of training observations arising from class g, a natural modification of this k-nearest-neighbour model is to reweight the neighbourhood sizes by a g = p g n/n g . This leads to the modified model

f (y|X, β, k) = exp β i a yi ∼ k i δ yi (y ) k Z(β, k) .
This modification is useful in practice when we are dealing with stratified surveys. In the following, however, we assume that a g = 1 for all g = 1, . . . , G.

Predictive perspective

When based on the conditional expression (4), the predictive distribution of a new unclassified observation y n+1 given its covariate x n+1 and the training sample (y, X) is, for g = 1, . . . , G,

P(y n+1 = g|x n+1 , y, X, β, k) ∝ exp    β/k   ∼ k (n+1) δ g (y ) + (n+1)∼ k δ y (g)      , (5)
where

∼ k (n+1)
δ g (y ) and

(n+1)∼ k δ y (g)
are the numbers of observations in the training dataset from class g among the k nearest neighbours of x n+1 and among the observations for which x n+1 is a k-nearest neighbour, respectively. This predictive distribution can then be incorporated in the Bayesian inference process by considering the joint posterior of (β, k, y n+1 ) and by deriving the corresponding marginal posterior distribution of y n+1 .

While this model provides a sound statistical basis for the k-nearest-neighour methodology as well as a means of assessing the uncertainty of the allocations to classes of unclassified observations, and while it corresponds to a true, albeit unavailable, joint distribution, it can be criticised from a Bayesian point of view in that it suffers from a lack of statistical coherence (in the sense that the information contained in the sample is not used in the most efficient way) when multiple classifications are considered. Indeed, the k-nearest-neighbour methodology is invariably used in a repeated manner, either jointly on a sample (x n+1 , . . . , x n+m ) or sequentially. Rather than assuming simultaneously dependence in the training sample and independence in the unclassified sample, it would be more sensible to consider the whole collection of points as issuing from a single joint model of the form given by ( 3), but with some having their class missing at random. Always reasoning from a Bayesian point of view, addressing jointly the inference on the parameters (β, k) and on the missing classes (y n+1 , . . . , y n+m )i.e. assuming exchangeability between the training and the unclassified datapoints-certainly makes sense from a foundational perspective as a correct probabilistic evaluation and it does provide a better assessment of the uncertainty about the classifications as well as about the parameters.

Unfortunately, this more global and arguably more coherent perspective is mostly unachievable if only for computational reasons, since the set of the missing class vector (y n+1 , . . . , y n+m ) is of size G m . It is practically impossible to derive an efficient simulation algorithm that would correctly approximate the joint probability distribution of both parameters and classes, especially when the number m of unclassified points is large. We will thus adopt the more ad hoc approach of dealing separately with each unclassified point in the analysis, because this simply is the only realistic way. This perspective can also be justified by the fact that, in realistic machine learning set-ups, the unclassified data (y n+1 , . . . , y n+m ) mostly occur in a sequential environment with, furthermore, the true value of y n+1 being revealed before y n+2 is observed.

In the following sections, we mainly consider the case G = 2 as in [START_REF] Holmes | Likelihood inference in nearest-neighbour classification models[END_REF], because this is the only case where we can conduct a full comparison between different approximation schemes, but we indicate at the end of Section 3.4 how a Gibbs sampling approximation allows for a realistic extension to larger values of G, as illustrated in Section 4.

Bayesian inference and the normalisation problem

Given the joint model ( 3) for (y 1 , . . . , y n+1 ), Bayesian inference can be conducted in a standard manner [START_REF] Robert | The Bayesian Choice[END_REF], provided computational difficulties related to the unavailability of the normalising constant can be solved. Indeed, as stressed in the previous section, from a Bayesian perspective, the classification of unclassified points can be based on the marginal predictive (or posterior) distribution of y n+1 obtained by integration over the conditional posterior distribution of the parameters, namely, for g = 1, 2,

P(y n+1 = g|x n+1 , y, X) = k P(y n+1 = g|x n+1 , y, X, β, k)π(β, k|y, X) dβ , (6)
where π(β, k|y, X) ∝ f (y|X, β, k)π(β, k) is the posterior distribution of (β, k) given the training dataset (y, X).

While other choices of prior distributions are available, we choose for (k, β) a uniform prior on the compact support {1, . . . , K} × [0, β max ]. The limitation on k is imposed by the structure of the training dataset in that K is at most equal to the minimal class size, min(n 1 , n 2 ), while the limitation on β, β < β max , is customary in Boltzmann models, because of phase-transition phenomena [START_REF] Møller | Spatial Statistics and Computational Methods[END_REF]: when β is above a certain value, the model becomes "all black or all white", i.e. all y i 's are either equal to 1 or to 2. (This is illustrated in Figure 5 below by the convergence of the expectation of the number of identical neighbours to k.) The determination of β max is obviously problem-specific and needs to be solved afresh for each new dataset since it depends on the topology of the neighbourhood. It is however straighforward to implement in that a Gibbs simulation of (3) for different values of β quickly exhibits the "black-or-white" features.

MCMC steps

Were the posterior distribution π(β, k|y, X) available (up to a normalising constant), we could design an MCMC algorithm that would produce a Markov chain approximating a sample from this posterior [START_REF] Robert | Monte Carlo Statistical Methods[END_REF], for example through a Gibbs sampling scheme based on the full conditional distributions of both k and β. However, because of the associated representation (4), the conditional distribution of β is nonstandard and we need to resort to a hybrid sampling scheme in which the exact simulation from π(β|k, y, X) is replaced with a single Metropolis-Hastings step. Furthermore, use of the full conditional distribution for k can impose fairly severe computational constraints. Indeed, for a given value β (t) , computing the posterior f (y|X, β (t) , i)π(β (t) , i), for i = 1, . . . , K, requires computations of order O(KnG), once again because of the likelihood representation. A faster alternative is to use a hybrid step for both β and k: in this way, we only need to compute the full conditional distribution of k for one new value of k, modulo the normalising constant. An alternative to Gibbs sampling is to use a random walk Metropolis-Hastings algorithm: both β and k are then updated using random walk proposals. Since β ∈ (0, β max ) is constrained, we first introduce a logistic reparameterisation of β, β = β max exp(θ) (exp(θ) + 1) , and then propose a normal random walk on the θ's, θ ∼ N (θ (t) , τ 2 ). For k, we use instead a uniform proposal on the 2r neighbours of k (t) , namely {k (t) -r, . . . , k (t) -1, k (t) + 1, . . . k (t) + r} {1, . . . , K}. This proposal distribution with probabiltity density Q r (k, •), with k ∼ Q r (k (t-1) , •), thus depends on a parameter r ∈ {1, . . . , K} that needs to be calibrated so as to aim at optimal acceptance rates, as does τ 2 . The acceptance probability in the Metropolis-Hastings algorithm is thus

ρ = f (y|X, β , k )π(β , k ) Q r (k (t-1) , k ) f (y|X, β (t-1) , k (t-1) )π(β (t-1) , k (t-1) ) Q r (k , k (t-1) ) × exp(θ ) (1 + exp(θ )) 2 exp(θ (t-1) ) (1 + exp(θ (t-1) )) 2 ,
where the second ratio is the ratio of the Jacobians due to the reparameterisation.

Once the Metropolis-Hastings algorithm has produced a satisfactory sequence of (β, k)'s, the Bayesian prediction for an unobserved class y n+1 associated with x n+1 is derived from (6). In fact, if we use a 0 -1 loss function [START_REF] Robert | The Bayesian Choice[END_REF] for predicting y n+1 , namely L(ŷ n+1 , y n+1 ) = I ŷn+1 =yn+1 , the Bayes estimator ŷπ n+1 is the most probable class g according to the posterior predictive (6). The associated measure of uncertainty is then the posterior expected loss, P(y n+1 = g|x n+1 , y, X).

Explicit calculation of ( 6) is obviously impossible and this distribution must be approximated from the MCMC chain {(β, k) (1) , . . . , (β, k) (M ) } simulated above, namely by

M -1 M i=1 P y n+1 = g|x n+1 , y, X, (β, k) (i) . (7)
Unfortunately, since (3) involves the intractable constant Z(β, k), the above schemes cannot be implemented as such and we need to replace f with a more manageable target. We proceed below through three different approaches that try to overcome this difficulty, postponing the comparison till Section 3.5.

Pseudo-likelihood approximation

A first solution, dating back to [START_REF] Besag | Spatial interaction and the statistical analysis of lattice systems (with discussion)[END_REF], is to replace the true joint distribution with the pseudo-likelihood, defined as

f (y|X, β, k) = n i=1 exp β/k ∼ k i δ yi (y ) + i∼ k δ y (y i ) 2 g=1 exp β/k ∼ k i δ g (y ) + i∼ k δ y (g) (8)
and made up of the product of the (true) conditional distributions associated with (3). The true posterior distribution π(β, k|y, X) is then replaced with

π(β, k|y, X) ∝ f (y|X, β, k)π(β, k) ,
and used as such in all steps of the MCMC algorithm drafted above. The predictive distribution P(y n+1 = g|x n+1 , y, X) is correspondingly approximated by ( 7), based on the pseudo-sample thus produced.

While this replacement of the true distribution with the pseudo-likelihood approximation induces a bias in the estimation of (k, β) and in the predictive performance of the Bayes procedure, it has been intensively used in the past, if only because of its availability and simplicity. For instance, [START_REF] Holmes | Likelihood inference in nearest-neighbour classification models[END_REF] built their pseudo-joint distribution on such a product (with the difficulty that the components of the product were not true conditionals). As noted in [START_REF] Friel | Likelihood estimation and inference for the autologistic model[END_REF], pseudo-likelihood estimation can be very misleading and we will describe its performance in more detail in Section 3.5. (To the best of our knowledge, this Bayesian evaluation has not been conducted before.)

As illustrated on Figure 4 for Ripley's benchmark data, the random walk Metropolis-Hastings algorithm detailed above performs satisfactorily with the pseudo-likelihood approximation, even though the mixing is slow (cycles can be spotted on the bottom left graph). On that dataset, the pseudo-maximum-i.e., the maximum of (8)-is achieved for k = 53 and β = 2.28. If we use the last 10, 000 iterations of this MCMC run, the prediction performance of ( 7) is such that the error rate on the test set of 1000 points is 8.7%. Figure 4 also indicates how limited the information is about k. (Note that we settled on the value β max = 4 by trial-and-error.)

Path sampling

A now-standard approach to the estimation of normalising constants is path sampling, described in Gelman and Meng (1998) (see also [START_REF] Chen | Monte Carlo Methods in Bayesian Computation[END_REF], and derived from the Ogata (1989) method, in which the ratio of two normalising constants, Z(β , k)/Z(β, k), can be decomposed as an integral to be approximated by Monte Carlo techniques.

The basic derivation of the path sampling approximation is that, if

S(y) = i ∼ k i δ yi (y )/k , then Z(β, k) = y exp [βS(y)]
and

∂Z(β, k) ∂β = y S(y) exp[βS(y)] = Z(β, k) y S(y) exp(βS(y)) Z(β, k) = Z(β, k) E β [S(y)] .
Therefore, the ratio Z(β, k)/Z(β , k) can be derived from an integral, since for each value of (β, k)). On these graphs, the black curve is based on linear interpolation of the expectation and the red curve on second-order spline interpolation.

log {Z(β, k)/Z(β , k)} = β β E u,k [S(y)] du ,
which is easily evaluated by a numerical approximation.

The practical drawback with this approach is that each time a new ratio is to be computed, that is, at each step of a hybrid Gibbs scheme or of a Metropolis-Hastings proposal, an approximation of the above integral needs to be produced. A further step is thus necessary for path sampling to be used: we approximate the function Z(β, k) only once for each value of k and for a few selected values of β, and later we use numerical interpolation to extend the function to other values of β. Since the function Z(β, k) is very smooth, the degree of additional approximation is quite limited. Given that this approximation is only to be computed once, the resulting Metropolis-Hastings algorithm is very fast, as well as being efficient if enough care is taken with the approximation by checking that the slope of Z(β, k) is sufficiently smooth from one value of β to the next. (We stress however that the computational cost required to produce those approximations is fairly high, because of the joint approximation in (β, k).)

We illustrate this approximation using Ripley's benchmark dataset. Figure 5 provides the approximated expectations E β,k [S(y)] for a range of values of β and for two values of k. Within the expectation, the y's are simulated using a systematic scan Gibbs sampler, because using the perfect sampling scheme elaborated below in Section 3.4 makes little sense when only one expectation needs to be computed. As seen from this comparative graph, when β is small, the Gibbs sampler gives good mixing performance, while, for larger values, it has difficulty in converging, as illustrated by the poor fit on the right-hand plot when k = 125. The explanation is that the model is getting closer to the phase-transition boundary in that case.

For the approximation of Z(β, k), we use the fact that E β,k [S(y)] is known when β = 0, namely E 0,k [S(y)] = n/2. We can thus represent log{Z(β, k)} as

n log 2 + β 0 E u,k [S(y)] du
and use numerical integration to approximate the integral. As shown on Figure 6, which uses a bilinear interpolation based on a 50 × 12 grid of values of (β, k), the approximated constant Z(β, k) is mainly constant in k.

Once Z(β, k) has been approximated, we can use the genuine MCMC algorithm of Section 3.1 fairly easily, the main cost of this approach being thus in the approximation of Z(β, k). Figure 7 illustrates the output of the MCMC sampler for Ripley's benchmark, to be compared with Figure 4. A first item of interest is that the chain mixes much more rapidly(in terms of iterations) than its pseudo-likelihood counterpart. A more important point is that the range and shape of the approximations to both marginal posterior distributions differ widely between the two methods, a feature discussed in Section 3.5. When this output of the MCMC sampler is used for prediction purposes in (7), the error rate for Ripley's test set is equal to 8.5%.

Perfect sampling implementation and Gibbs approximation

A completely different approach to handling missing normalising constants has been developed recently by [START_REF] Møller | An efficient Markov chain Monte Carlo method for distributions with intractable normalising constants[END_REF] and is based on an auxiliary variable idea. If we introduce an auxiliary variable z on the same state space as y, with arbitrary conditional density g(z|β, k, y), and if we consider the joint posterior then simulating (β, k, z) from this posterior is equivalent to simulating (β, k) from the original posterior since z integrates out. If we now run a Metropolis-Hastings algorithm on this augmented scheme, with q 1 an arbitrary proposal density on (β, k) and with

π(β, k, z|y) ∝ π(β, k, z, y) = g(z|β, k, y) × f (y|β, k) × π(β, k) ,
q 2 (β , k , z |β, k, z) = q 1 (β , k |β, k, y)f (z |β , k ) ,
as the joint proposal on (β, k, z) (i.e., simulating z directly from the likelihood), the Metropolis-Hastings ratio associated with q 2 is

Z(β, k) Z(β , k) exp (β S(y)/k ) π(β , k ) exp (βS(y)/k) π(β, k) g(z |β , k , y) g(z|β, k, y) × q 1 (β, k|β , k, y) exp (βS(z)/k) q 1 (β , k |β, k, y) exp (β S(z)/k ) Z(β , k ) Z(β, k) ,
which means that the constants Z(β, k) and Z(β , k ) cancel out. The method of [START_REF] Møller | An efficient Markov chain Monte Carlo method for distributions with intractable normalising constants[END_REF] can thus be calibrated by the choice of the artificial target g(z|β, k, y) on the auxiliary variable z, and the authors advocate the choice

g(z|β, k, y) = exp βS(z)/ k Z( β, k),
as reasonable, where ( β, k) is a preliminary estimate, such as the maximum pseudo-likelihood estimate. While we follow this recommendation, we stress that the choice of ( β, k) is paramount for good performance of the algorithm, as explained below. The alternative of setting a target g(z|β, k, y) that truly depends on β and k is appealing but faces computational difficulties in that the most natural proposals involve normalising constants that cannot be computed.

Obviously, this approach of [START_REF] Møller | An efficient Markov chain Monte Carlo method for distributions with intractable normalising constants[END_REF] also has a major drawback, namely that the auxiliary variable z must be simulated from the distribution f (z|β, k) itself. However, there have been many developments in the simulation of Ising models, from [START_REF] Besag | Spatial interaction and the statistical analysis of lattice systems (with discussion)[END_REF] to [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF], and the particular case G = 2 allows for exact simulation of f (z|β, k) using perfect sampling. We refer the reader to [START_REF] Häggström | Finite Markov Chains and Algorithmic Applications[END_REF], [START_REF] Møller | Spatial Statistics and Computational Methods[END_REF], [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF] and Robert and Casella (2004, Chapter 13) for details of this simulation technique and for a discussion of its limitations. Without entering into technical details, we comment that, in the case of model ( 3) with G = 2, there also exists a monotone implementation of the Gibbs sampler that allows for a practical implementation of the perfect sampler [START_REF] Kendall | Perfect simulation using dominating processes on ordered spaces, with application to locally stable point processes[END_REF][START_REF] Berthelsen | Likelihood and non-parametric Bayesian MCMC inference for spatial point processes based on perfect simulation and path sampling[END_REF]. More precisely, we can use a coupling-from-the-past strategy [START_REF] Propp | Coupling from the past: a user's guide[END_REF]: in this setting, starting from the saturated situations in which the components of z are either all equal to 1 or all equal to 2, it is sufficient to monitor both associated chains further and further into the past until they coalesce by time 0. The sandwiching property of [START_REF] Kendall | Perfect simulation using dominating processes on ordered spaces, with application to locally stable point processes[END_REF] and the monotonicity of the Gibbs sampler ensure that all other chains associated with arbitrary starting values for z will also have coalesced by then. The only difficulty with this perfect sampler is the phase-transition phenomenon, which means that, for very large values of β, the convergence performance of the coupling from the past sampler deteriorates quite rapidly, a fact also noted in [START_REF] Møller | An efficient Markov chain Monte Carlo method for distributions with intractable normalising constants[END_REF] for the Ising model. We overcome this difficulty by using an additional accept-reject step based on smaller values of β that avoids this explosion in the computational time.

As shown on Figure 8, a poor choice for ( β, k) leads to very unsatisfactory performance with the algorithm. Starting from the pseudo-likelihood estimate and using this very value for the plug-in value ( β, k), we obtain a Markov chain with a very low energy and a very high rejection rate. However, use of the estimate ( k, β) = (13, 1.45) resulting from this poor run does improve considerably the performance of the algorithm, as shown by Figure 9. In this setting, the predictive error rate on the test dataset is equal to 0.084.

While this elegant solution based on an auxiliary variable completely removes the issue of the normalising constant, it faces several computational difficulties. First, as noted above, the choice of the artificial target g(z|β, k, y) is driving the algorithm and plug-in estimates need to be reassessed periodicaly. Second, perfect simulation from the distribution f (z|β, k) is extremely costly and may fail if β is close to the phase-transition boundary. Furthermore, the numerical value of this critical point is not known beforehand. Finally, the extension of the perfect sampling scheme to more than G = 2 classes has not yet been achieved.

For these different reasons, we advocate the substitution of a Gibbs sampler for the above perfect sampler in order to achieve manageable computing performance. If we replace the perfect sampling step with 500 (complete) iterations of the corresponding generic Gibbs sampler on z, the computing time is linear in the number n of observations and the results are virtually the same. One has to remember that the simulation of z is of second-order with respect to the original problem of simulating the posterior distribution of (β, k), since z is an auxiliary variable introduced to overcome the computation of the normalising constant. Therefore, the additional uncertainty induced by the use of the Gibbs sampler is far from severe. Figure 9 compares the Gibbs solution with the perfect sampling implementation and it shows how little loss is incurred by the use of the less expensive Gibbs sampler, while the gain in computing time is enormous. For 50, 000 iterations, the time required to run the Gibbs sampler is approximately 20 minutes, compared with more than a week for the corresponding perfect sampler (under the same C environment on the same machine).

Evaluation of the pseudo-likelihood approximation

Given that the above alternatives can all be implemented for small values of n, it is of direct interest to compare them in order to evaluate the effect of the pseudo-likelihood approximation. As demonstrated in the previous section, using Ripley's benchmark with a training set of 250 points, we are indeed able to run a perfect sampler over the range of possible β's, and this implementation gives a sampler in which the only approximation is due to running an MCMC sampler (a feature common to all three versions).

Histograms, for the same dataset, of simulated β's, conditional or unconditional, on k show gross misrepresentation of the samples produced by the pseudo-likelihood approximation; see Figures 10 and11. (The comparison for a fixed value of k was obtained directly by setting k to a fixed value in all three approaches and running the corresponding MCMC algorithms.) It could of course be argued that the defect lies with the path sampling evaluation of the constant, but this approach strongly coincides with the perfect sampling implementation, as showed on both figures. There is thus a fundamental discrepancy in using the pseudo-likelihood approximation; in other words, the pseudo-likelihood approximation defines a clearly different posterior distribution on (β, k).

As exhibited on Figure 10, the larger k is, the worse is this discrepancy, whereas Figure 11 shows that both β and k are significantly overestimated by the pseudo-likelihood approximation. (It is quite natural to find such a correlation between β and k when we realise that the likelihood depends mainly on β/k.) We can also note that the correspondence between path and perfect approximations is not absolute in the case of k, a difference that may be attributed to slower convergence in one or both samplers.

In order to assess the comparative predictive properties of both approaches, we also provide a comparison of the class probabilities P(y = 1|x, y, X) estimated at each point of the test sample. As shown by Figure 12, the predictions are quite different for values in the middle of the range, with no clear bias direction in using

Illustration on real datasets

In this Section, we illustrate the behaviour of the proposed methodology on some benchmark datasets.

We first describe the calibration of the algorithm used on each dataset. As starting value for the Gibbs approximation in the Møller scheme, we use the maximum pseudo-likelihood estimate. The Gibbs sampler is iterated 500 times as an approximation to the perfect sampling step. After 10,000 iterations, we modify the plug-in estimate using the current average and then we run 50,000 more iterations of the algorithm.

The first dataset is borrowed from the MASS library of R. It consists in the records of 532 Pima Indian women who were tested by the U.S. National Institute of Diabetes and Digestive and Kidney Diseases for diabetes. Seven quantitative covariates were recorded, along with the presence or absence of diabetes. The data are split at random into a training set of 200 women, including 68 diagnosed with diabetes, and a test set of the remaining 332 women, including 109 diagnosed with diabetes. The performance for various values of k on the test dataset is given in Table 2. If we use a standard leave-one-out cross-validation for selecting k (using only the training dataset), then there are 10 consecutive values of k leading to the same error rate, namely the range 57-66. The results are provided in Figure 13. Note that the simulated values of k tend to avoid the region found by the cross-validation procedure. One possible reason for this discrepancy is that, as noted in Section 2.2, the likelihood for our joint model is not directly equivalent to the k-nearest-neighbour objective function, since mutual neighbours are weighted twice as heavily as single neighbours in this likelihood. Over the final 20, 000 iterations, the prediction error is 0.209, quite in line with the k-nearest-neighbour solution in Table 2.

To illustrate the ability of our method to consider more than two classes, we also used the benchmark dataset forensic glass fragments, studied in [START_REF] Ripley | Neural networks and related methods for classification (with discussion)[END_REF]. This dataset involves nine covariates and six classes some of which are rather rare. Following the recommendation made in [START_REF] Ripley | Neural networks and related methods for classification (with discussion)[END_REF], we coalesced some classes to reduce the number of classes to four. We then randomly partitioned the dataset to obtain 89 individuals in the training dataset and 96 in the testing dataset. Leave-one-out cross-validation leads us to choose the value k = 17. The error rate of the 17-nearest-neighbour procedure on the test dataset is 0.35, whereas, using our procedure, we obtain an error rate of 0.29. The substantial gain from using our approach can be partially explained by the fact that the value of k chosen by the cross-validation procedure is much larger than those explored by our MCMC sampler.

Conclusions

While the probabilistic background to a Bayesian analysis of k-nearest-neighbour methods was initiated by [START_REF] Holmes | Likelihood inference in nearest-neighbour classification models[END_REF], the present paper straightens the connection between the original technique and a true probabilistic model by defining a coherent probabilistic model on the training dataset. This new model (3) then provides a sound setting for Bayesian inference and for evaluating not just the most likely allocations for the test dataset but also the uncertainty that goes with them. The advantages of using a probabilistic environment are clearly demonstrated: it is only within this setting that tools like predictive maps as in Figure 3 can be constructed. This obviously is a tremendous bonus for the experimenter, since boundaries between most likely classes can thus be estimated and regions can be established in which allocation to a specific class or to any class is uncertain. In addition, the probabilistic framework allows for a natural and integrated analysis of the number of neighbours involved in the class allocation, in a standard model-choice perspective. This perspective can be extended to the choice of the most significant components of the covariate x, even though this possibility is not explored in the current paper.

The present paper also addresses the computational difficulties related to this approach, namely the wellknown issue of the intractable normalising constant. While this has been thoroughly discussed in the literature, our comparison of three independent approximations leads to the strong conclusion that the pseudo-likelihood approximation is not to be trusted for training sets of moderate size. Furthermore, while the path sampling and perfect sampling approximations are useful in establishing this conclusion, they cannot be advocated at the operational level, but we also demonstrate that a Gibbs sampling alternative to the perfect sampling scheme of [START_REF] Møller | An efficient Markov chain Monte Carlo method for distributions with intractable normalising constants[END_REF] is both operational and practical.

Figure 2 :

 2 Figure 2: Misclassification leave-one-out error rate as a function of k for Ripley's training dataset.

Figure 3 :

 3 Figure3: (left) Level sets of the predictive probability to be in the black class, ranging from high (white) to low (red), and (right) consequences of the comparison with 0.5 of the 95% credibility intervals for the predictive probabilities. (These plots are based on an MCMC sample whose derivation is explained in Section 3.4.)

Figure 4 :

 4 Figure 4: Output of a random walk Metropolis-Hastings algorithm based on the pseudo-likelihood approximation of the normalising constant for 50, 000 iterations, with a 40, 000 iteration burn-in stage, and τ 2 = 0.05, r = 3. (top) sequence and marginal histogram for β when β max = 4 and (bottom) sequence and marginal barplot for k.

Figure 5 :

 5 Figure5: Approximation of the expectation E β,k [S(y)] for Ripley's benchmark, where the β's are equally spaced between 0 and β max = 4, and for k = 1 (left) and k = 125 (right) (10 4 iterations with 500 burn-in steps for each value of (β, k)). On these graphs, the black curve is based on linear interpolation of the expectation and the red curve on second-order spline interpolation.

Figure 6 :

 6 Figure 6: Approximation of the normalising constant Z(β, k) for Ripley's dataset where the β's are equally spaced between 0 and β max = 4, and k = 1, 10, 20, . . . , 110, 125 (based on 10 4 Monte Carlo iterations with 500 burn-in steps, and bilinear interpolation).

Figure 7 :

 7 Figure 7: Output of a random walk Metropolis-Hastings algorithm based on the path sampling approximation of the normalising constant for 50, 000 iterations, with a 40, 000 iteration burn-in stage and τ 2 = 0.05, r = 3. (top) sequence and marginal histogram for β when β max = 4 and (bottom) sequence and marginal barplot for k.

Figure 8 :Figure 9 :

 89 Figure 8: Output of a random walk Metropolis-Hastings algorithm based on the perfect sampling elimination of the normalising constant for a pseudo-likelihood plug-in estimate ( k, β) = (53, 2.28) and 20, 000 iterations, with a 10, 000 burn-in stage, β max = 4 and τ 2 = 0.05, r = 3: (top) sequence and marginal histogram for β and (bottom) sequence and marginal barplot for k.

Figure

  Figure Comparison of the approximations to the posterior distribution of β based on the pseudo (red), the path (green) and the perfect (yellow) schemes for Ripley's benchmark and k = 1, 10, 70, 125, for 20, 000 iterations and 10, 000 burn-in.

Figure 13 :

 13 Figure13: Pima Indian diabetes study based on 50, 000 iterations of the Gibbs-Møller sampling scheme with τ 2 = 0.05, r = 3, β max = 4, and K = 68.

  Figure1: Training (top) and test (bottom) groups for Ripley's benchmark: the points in red are those for which the label is equal to 1 and the points in black are those for which the label is equal to 2. the performance of the standard k-nearest-neighbour method on the test dataset for several values of k. The overall misclassification leave-one-out error rate on the training dataset as k varies is provided in Figure2and it shows that this criterion is not very discriminating for this dataset, with little variation for a wide range of values of k and with several values of k achieving the same overall minimum, namely17, 18, 35, 36, 45, 46, 51, 52, 53 and 54. There are therefore ten different values of k in competition. This range of values is an indicator of potential gains when averaging over k, and hence calls for a Bayesian perspective.

			k Misclassification
					error rate
			1		0.150
			3		0.134
			15		0.095
			17		0.087
			31		0.084
			54		0.081
	Table 1: k-nearest-neighbour performances on the Ripley test dataset
	Leave-one-out cross validation error rate	40 50 60
		30
		0	20	40	60	80	100	120
					k

Table 2 :

 2 Performance of k-nearest-neighbour methods on the Pima Indian test dataset.

	k Misclassification
		error rate
	1	0.316
	3	0.229
	15	0.226
	31	0.211
	57	0.205
	66	0.208

RR n°6173

http://www.inria.fr ISSN 0249-6399

Acknowledgements

The authors are grateful to Gilles Celeux for his numerous and insightful comments on the different perspectives offered by this probabilistic reassessment, as well as to the Associate Editor and to both referees for their constructive comments. Both JMM and CPR are also grateful to the Department of Statistics of the University of Glasgow for its warm welcome during various visits related to this work. This work was supported in part by the IST Programme of the European Community, under the PASCAL Network of Excellence, ST-2002-506778.

 4for the pseudo-likelihood approximation, in Figure 7 for the path sampling approximation and in Figure 9 for the perfect sampling approximation. pseudo-likelihood as an approximation. Note that the discrepancy may be substantial and may result in a large number of different classifications.