Keywords: Choix bayésien de modèles, algorithmes à sauts reversibles, méthodes de Monte-Carlo par Chaînes

Sur les difficultés associées à une méthode d'approximation de probabilités a posteriori de modèles Résumé : Dans le paradigme bayésien, Scott (2002) et Congdon (2006) ont proposé une méthode permettant d'approcher les probabilités a posteriori de modèles. Cette technique s'appuie sur des simulations a posteriori suivant les paramètres de chacun des modèles, ce qui est en pratique très appréciable. Bien qu'il soit possible d'approcher les probabilités a posteriori de modèles à partir de simulations suivant les paramètres de chaque modèle (Gelfand and Dey (1994) et Bartolucci et al. (2006)), nous montrons que la méthode de Scott (2002) et Congdon (2006) est biaisée et basée sur un argument fallacieux.

Introduction

Model selection is a fundamental statistical issue and a clear asset of the Bayesian methodology but it faces severe computational difficulties because of the requirement to explore simultaneously the parameter spaces of all models under comparison with enough of an accuracy to provide sufficient approximations to the posterior probabilities of all models. When [START_REF] Green | Reversible jump MCMC computation and Bayesian model determination[END_REF] introduced reversible jump techniques, it was perceived by the community as the second MCMC revolution in that it allowed for a valid and efficient exploration of the collection of models and the subsequent literature on the topic exploiting reversible jump MCMC is a testimony to the appeal of this method. Nonetheless, the implementation of reversible jump techniques in complex situations may face difficulties or at least inefficiencies of its own and, despite some recent advances in the devising of the jumps underlying reversible jump MCMC [START_REF] Brooks | Efficient construction of reversible jump Markov chain Monte Carlo proposal distributions (with discussion)[END_REF], the care required in the construction of those jumps often acts as a deterrent from its applications.

There are practical alternatives to reversible jump MCMC when the number of models under consideration is small enough to allow for a complete exploration of those models. Integral approximations using importance sampling techniques like those found in [START_REF] Gelfand | Bayesian model choice: asymptotics and exact calculations[END_REF], based on an harmonic mean representation of the marginal densities, and in [START_REF] Gelman | Simulating normalizing constants: From importance sampling to bridge sampling to path sampling[END_REF], focussing on the optimised selection of the importance function, are advocated as potential solutions, see [START_REF] Chen | Monte Carlo Methods in Bayesian Computation[END_REF] for a detailed entry. The reassessment of those methods by [START_REF] Bartolucci | Efficient Bayes factor estimation from the reversible jump output[END_REF] showed the connection between a virtual reversible jump MCMC and importance sampling (see also [START_REF] Chopin | Contemplating evidence: properties, extensions of, and alternatives to nested sampling[END_REF]. In particular, those papers demonstrated that the output of MCMC samplers on each single model could be used to produce approximations of posterior probabilities of those models, via some importance sampling methodologies also related to [START_REF] Newton | Approximate Bayesian inference by the weighted likelihood boostrap (with discussion)[END_REF].

In [START_REF] Scott | Bayesian methods for hidden Markov models: recursive computing in the 21st Century[END_REF] and [START_REF] Congdon | Bayesian model choice based on Monte Carlo estimates of posterior model probabilities[END_REF], a new and straightforward method is advanced to compute posterior probabilities of models under scrutinity based solely on MCMC outputs restricted to single models. While this simplicity is quite appealing for the approximation of those probabilities, we believe that both proposals of [START_REF] Scott | Bayesian methods for hidden Markov models: recursive computing in the 21st Century[END_REF] and [START_REF] Congdon | Bayesian model choice based on Monte Carlo estimates of posterior model probabilities[END_REF] are inherently biased and we advance in this note several arguments towards this thesis. In addition, we notice that, to overcome the bias we thus exhibited, a valid solution would call for the joint simulation of parameters under all models (using priors or pseudo-priors), and this step would thus loose the primary appeal of the methods against the one proposed by [START_REF] Carlin | Bayesian model choice through Markov chain Monte Carlo[END_REF], from which both [START_REF] Scott | Bayesian methods for hidden Markov models: recursive computing in the 21st Century[END_REF] and [START_REF] Congdon | Bayesian model choice based on Monte Carlo estimates of posterior model probabilities[END_REF] are inspired.

The methods

In a Bayesian framework of model comparison (see, e.g., [START_REF] Robert | The Bayesian Choice[END_REF], given D models in competition, M k , with densities f k (y|θ k ), and prior probabilities ̺ k = P (M = k) (k = 1, . . . , D), the posterior probabilities of the models M k conditional on the data y are given by

P (M = k|y) ∝ ̺ k f k (y|θ k )π k (θ k ) dθ k ,
the proportionality term being given by the sum of the above and M denoting the unknown model index.

In the specific setup of hidden Markov models, the solution of Scott (2002, Section 4.1) is to generate simultaneously and independently D MCMC chains (θ

(t) k ) t , 1 ≤ k ≤ D ,
with stationary distributions π k (θ k |y) and to approximate P (M = k|y) by

̺k (y) ∝ ̺ k T t=1 f k (y|θ (t) k ) D j=1 ̺ j f j (y|θ (t) j ) ,
as reported in formula (21) of [START_REF] Scott | Bayesian methods for hidden Markov models: recursive computing in the 21st Century[END_REF], with the mention that "(21) averages the D likelihoods corresponding to each θ j over the life of the Gibbs sampler" (p.347), the later being understood as "independently sampled D parallel Gibbs samplers" (p.347).

From a more general perspective, the proposal of [START_REF] Congdon | Bayesian model choice based on Monte Carlo estimates of posterior model probabilities[END_REF] for an approximation of the P (M = k|y)'s follows both from [START_REF] Scott | Bayesian methods for hidden Markov models: recursive computing in the 21st Century[END_REF] approximation and from the pseudo-prior construction of [START_REF] Carlin | Bayesian model choice through Markov chain Monte Carlo[END_REF] that predated reversible jump MCMC by saturating the parameter space with an artificial simulation of all parameters at each iteration. However, due to a very special (and, we believe, mistaken) choice of pseudo-priors discussed below, Congdon's (2006, p.349) approximation of P (M = k|y) eventually reduces to the estimator

̺k (y) ∝ ̺ k T t=1 f k (y|θ (t) k )π k (θ (t) k ) D j=1 ̺ j f j (y|θ (t) j )π j (θ (t) j ) ,
where the θ (t) k 's are samples from π k (θ k |y) (or approximate samples obtained by an MCMC algorithm). (A very similar proposal is found in [START_REF] Congdon | Model weights for model choice and averaging[END_REF] and, while some issues like the use of improper pseudo-priors discussed in Section 3.3 are corrected, the fundamental difficulty of simulating from the wrong target remains. We thus chose to address primarily the initial paper by [START_REF] Congdon | Bayesian model choice based on Monte Carlo estimates of posterior model probabilities[END_REF], especially since it generated follow-up papers like [START_REF] Chen | Bayesian model selection for heteroskedastic models[END_REF] and since its inherent simplicity is likely to appeal to unwary readers.)

Although both approximations ̺k (y) and ̺k (y) differ in their expressions, they fundamentally relate to the same notion that parameters from other models can be ignored when conditioning on the model index M . This approach is therefore bypassing the simultaneous exploration of several parameters spaces and restricts the simulation to marginal samplers on each separate model. This feature is very appealing since it cuts most of the complexity from the schemes both of [START_REF] Carlin | Bayesian model choice through Markov chain Monte Carlo[END_REF] and of [START_REF] Green | Reversible jump MCMC computation and Bayesian model determination[END_REF]. We however question the foundations of those approximations as presented in both [START_REF] Scott | Bayesian methods for hidden Markov models: recursive computing in the 21st Century[END_REF] and [START_REF] Congdon | Bayesian model choice based on Monte Carlo estimates of posterior model probabilities[END_REF] and advance below arguments that both authors are using incompatible versions of joint distributions on the collection of parameters that jeopardise the validity of the approximations.

Difficulties

The sections below expose the difficulties found with both methods, following the points made in [START_REF] Scott | Bayesian methods for hidden Markov models: recursive computing in the 21st Century[END_REF] and [START_REF] Congdon | Bayesian model choice based on Monte Carlo estimates of posterior model probabilities[END_REF], respectively. The fundamental difficulty with their approaches appears to us to be related to a confusion between the model dependent simulations and the joint simulations based on a pseudo-prior scheme as in [START_REF] Carlin | Bayesian model choice through Markov chain Monte Carlo[END_REF]. Once this difficulty is resolved, it appears that the approximation of P (M = k|y) by P (M = k|y) does require a joint simulation of all parameters and thus that the solutions proposed in [START_REF] Scott | Bayesian methods for hidden Markov models: recursive computing in the 21st Century[END_REF] and [START_REF] Congdon | Bayesian model choice based on Monte Carlo estimates of posterior model probabilities[END_REF] are of the same complexity as the proposal of [START_REF] Carlin | Bayesian model choice through Markov chain Monte Carlo[END_REF].

Incorrect marginals

We denote by θ = (θ 1 , . . . , θ D ) the collection of parameters for all models under consideration. Both [START_REF] Scott | Bayesian methods for hidden Markov models: recursive computing in the 21st Century[END_REF] and [START_REF] Congdon | Bayesian model choice based on Monte Carlo estimates of posterior model probabilities[END_REF] (1)

∝ D k=1 ̺ k f k (y|θ k ) j π j (θ j ) = D k=1 ̺ k m k (y) π k (θ k |y) j =k π j (θ j ) . (2)
In both papers, the θ (t) 's are instead simulated as independent outputs from the componentwise posteriors π k (θ k |y) and this divergence jeopardises the validity of the approximation. The error in their interpretations stems from the fact that, while the θ

(t)
k 's are (correctly) independent given the model index M , this independence does not hold once M is integrated out, which is the case in the above approximation P (M = k|y).

MCMC versus marginal MCMC

When [START_REF] Congdon | Bayesian model choice based on Monte Carlo estimates of posterior model probabilities[END_REF] defines a Markov chain (θ (t) ) at the top of page 349, he indicates that the components of θ (t) are made of independent Markov chains (θ (t) k ) simulated with MCMC samplers related to the respective marginal posteriors π k (θ k |y), following the approach of [START_REF] Scott | Bayesian methods for hidden Markov models: recursive computing in the 21st Century[END_REF]. The aggregated chain (θ (t) ) is thus stationary against the product of those marginals,

D k=1 π k (θ k |y) .
However, in the derivation of [START_REF] Carlin | Bayesian model choice through Markov chain Monte Carlo[END_REF], the model is defined in terms of (1) and the Markov chain should thus be constructed against (1), not against the product of the model marginals. Obviously, in the case of [START_REF] Congdon | Bayesian model choice based on Monte Carlo estimates of posterior model probabilities[END_REF], the fact that the pseudo-joint distribution does not exist because of the flat prior assumption (see Section 3.3) prevents this construction but, in the event the flat prior is replaced with a proper (pseudo-) prior (as in [START_REF] Congdon | Model weights for model choice and averaging[END_REF], the same statement holds: the probabilistic derivation of P (M = k|y) relies on the pseudo-prior construction and, to be valid, it does require the completion step at the core of [START_REF] Carlin | Bayesian model choice through Markov chain Monte Carlo[END_REF], where parameters need to be simulated from the pseudo-priors.

Similarly, in [START_REF] Scott | Bayesian methods for hidden Markov models: recursive computing in the 21st Century[END_REF], the target of the Markov chain (θ (t) , M (t) ) should be the distribution

P (θ, M = k|y) ∝ π k (θ k ) ̺ k f k (y|θ k ) j =k π j (θ j )
and the θ (t) j 's should thus be generated from the prior π j (θ j ) when M (t) = j-or equivalently from the corresponding marginal if one does not condition on M (t) , but simulating a Markov chain with stationary distribution ( 2) is certainly a challenge in many settings if the latent variable decomposing the sum is not to be used.

Since, in both [START_REF] Scott | Bayesian methods for hidden Markov models: recursive computing in the 21st Century[END_REF] and [START_REF] Congdon | Bayesian model choice based on Monte Carlo estimates of posterior model probabilities[END_REF], the (θ (t) )'s are not simulated against the correct target, the resulting averages of P (M = k|y, θ (t) ), ̺k (y) and ̺k (y), will both be biased, as demonstrated in the example of Section 3.4.

Improperty of the posterior

When resorting to the construction of pseudo-posteriors adopted by [START_REF] Carlin | Bayesian model choice through Markov chain Monte Carlo[END_REF], [START_REF] Congdon | Bayesian model choice based on Monte Carlo estimates of posterior model probabilities[END_REF] uses a flat prior as pseudo-prior on the parameters that are not in model M k . More precisely, the joint prior distribution on (θ, M ) is given by [START_REF] Congdon | Bayesian model choice based on Monte Carlo estimates of posterior model probabilities[END_REF] formula (2),

P (θ, M = k) = π k (θ k ) ̺ k j =k π(θ j |M = k) = π k (θ k ) ̺ k ,
which is indeed equivalent to assuming a flat prior as pseudo-prior on the parameters θ j that are not in model M k .

Unfortunately, this simplifying assumption has a dramatic consequence in that the corresponding joint posterior distribution of θ is never defined (as a probability distribution) since

π(θ|y) = D k=1 π k (θ k |y) P (M = k|y)
does not integrate to a finite value in any of the θ k 's (unless their support is compact). When [START_REF] Congdon | Bayesian model choice based on Monte Carlo estimates of posterior model probabilities[END_REF] points out "that it is not essential that the priors for P (θ j =k |M = k) are improper" (p.348), the whole issue is that they cannot be improper.

The fact that the posterior distribution on the saturated vector θ = (θ 1 , . . . , θ D ) does not exist obviously has dire consequences on the subsequent derivations, since a positive recurrent Markov chain with stationary distribution π(θ|y) cannot be constructed. Similarly, the fact that P (M = k|y) = P (θ, M = k|Y ) dθ does not hold any longer.

Note that [START_REF] Scott | Bayesian methods for hidden Markov models: recursive computing in the 21st Century[END_REF] does not follow the same track: when defining the pseudo-priors in his formula (20), he uses the product definition

1 P (θ, M = k) = π k (θ k ) ̺ k j =k π j (θ j ) ,
which (seemingly) means that the true priors are also used as pseudo-priors across all models. However, we stress that [START_REF] Scott | Bayesian methods for hidden Markov models: recursive computing in the 21st Century[END_REF] does not refer to the construction of [START_REF] Carlin | Bayesian model choice through Markov chain Monte Carlo[END_REF] in his proposal.

Illustration

We now proceed through a toy example where all posterior quantities can be computed in order to evaluate the bias brought by both approximations.

Example 1. Consider the case when a model M 1 : y|θ ∼ U(0, θ) with a prior θ ∼ Exp( 1) is opposed to a model M 2 : y|θ ∼ Exp(θ) with a prior θ ∼ Exp(1). We also assume equal prior weights on both models:

̺ 1 = ̺ 2 = 0.5.
The marginals are then

m 1 (y) = ∞ y θ -1 e -θ dθ = E 1 (y) ,
where E 1 denotes the exponential integral function tabulated both in Mathematica and in the GSL library, and

m 2 (y) = ∞ 0 θe -θ(y+1) dθ = 1 (1 + y) 2 .
For instance, when y = 0.2, the posterior probability of M 1 is thus equal to

P (M = 1|y) = m 1 (y)/{m 1 (y) + m 2 (y)} = E 1 (y)/{E 1 (y) + (1 + y) -2 } ≈ 0.6378 ,
while, for y = 0.9, it is approximately 0.4843. This means that, in the former case, the Bayes factor of M 1 against M 2 is B 12 ≈ 1.760, while for the later, it decreases to B 12 ≈ 0.939.

The posterior on θ in model M 2 is a gamma Ga(2, 1 + y) distribution and it can thus be simulated directly. For model M 1 , the posterior is proportional to θ -1 exp(-θ) for θ larger than y and it can be simulated using a standard accept-reject algorithm based on an exponential Exp(1) proposal translated by y.

Using simulations from the true (marginal) posteriors and the approximation of [START_REF] Congdon | Bayesian model choice based on Monte Carlo estimates of posterior model probabilities[END_REF], the numerical value of ̺1 (y) based on 10 6 simulations is 0.7919 when y = 0.2 and 0.5633 when y = 0.9, which translates into Bayes factors of 3.805 and of 1.288, respectively. For the approximation of [START_REF] Scott | Bayesian methods for hidden Markov models: recursive computing in the 21st Century[END_REF], the numerical value of ̺1 (y) is 0.6554 (corresponding to a Bayes factor of 1.898) when y = 0.2 and 0.6789 when y = 0.9 (corresponding to a Bayes factor of 2.11), based on the same simulations. Note that in the case y = 0.9, a selection based on either approximation of the Bayes factor would select the wrong model.

If we use instead a correct simulation from the joint posterior (2), which can be achieved by using a Gibbs scheme with target distribution P (θ, M = k|y), we then get a proper MCMC approximation to the posterior probabilities by the P (M = k|y)'s. For instance, based on 10 6 simulations, the numerical value of P (M = 1|y) when y = 0.2 is 0.6370, while, for y = 0.9, it is 0.4843. Note that, due to the impropriety difficulty exposed in Section 3.3, the equivalent correction for [START_REF] Congdon | Bayesian model choice based on Monte Carlo estimates of posterior model probabilities[END_REF] scheme cannot be implemented.

In Figure 1, the three approximations are compared to the exact value of P (M = 1|y) for a range of values of y. The correct simulation produces a graph that is indistinguishable from the true probability, while [START_REF] Congdon | Bayesian model choice based on Monte Carlo estimates of posterior model probabilities[END_REF] approximation stays within a reasonable range of the true value and [START_REF] Scott | Bayesian methods for hidden Markov models: recursive computing in the 21st Century[END_REF] surprisingly drifts apart for most values of y. ◭

The correspondence of what is essentially Carlin and Chib's (1995) scheme with the true numerical value of the posterior probabilities is obviously unsurprising in this toy example but more advanced setups see the approximation degenerate, since the simulations from the prior are most often inefficient, especially when the number of models under comparison is large. This is the reason why [START_REF] Carlin | Bayesian model choice through Markov chain Monte Carlo[END_REF] introduced pseudo-priors that were closer approximations to the true posteriors. [START_REF] Scott | Bayesian methods for hidden Markov models: recursive computing in the 21st Century[END_REF] approximation (in blue), [START_REF] Congdon | Bayesian model choice based on Monte Carlo estimates of posterior model probabilities[END_REF] approximation (in green), and correction of [START_REF] Scott | Bayesian methods for hidden Markov models: recursive computing in the 21st Century[END_REF] approximation (in red), indistinguishable from the true value (based on N = 10 6 simulations).

  start from the representation P (M = k|y) = P (M = k|y, θ)π(θ|y) dθ to justify the approximation P (M = k|y) = T t=1 P (M = k|y, θ (t) )/T . This is indeed an unbiased estimator of P (M = k|y) provided the θ (t) 's are generated from the correct (marginal) posterior π(θ|y) = D k=1 P (θ, M = k|y)

Figure 1 :

 1 Figure1: Comparison of three approximations of P (M = 1|y) with the true value (in black):[START_REF] Scott | Bayesian methods for hidden Markov models: recursive computing in the 21st Century[END_REF] approximation (in blue),[START_REF] Congdon | Bayesian model choice based on Monte Carlo estimates of posterior model probabilities[END_REF] approximation (in green), and correction of[START_REF] Scott | Bayesian methods for hidden Markov models: recursive computing in the 21st Century[END_REF] approximation (in red), indistinguishable from the true value (based on N = 10 6 simulations).

RR n➦ 0123456789

The indices on the priors have been added to make notations coherent, since[START_REF] Scott | Bayesian methods for hidden Markov models: recursive computing in the 21st Century[END_REF] denotes all priors with the same letter p.

Acknowledgements

Both authors are grateful to Brad Carlin for helpful discussions and to Antonietta Mira for providing a perfect setting for this work during the ISBA-IMS "MCMC'ski 2" conference in Bormio, Italy. The second author is also grateful to Kerrie Mengersen for her invitation to "Spring Bayes 2007" in Coolangatta, Australia, that started our reassessment of those papers. This work had been supported by the Agence Nationale de la Recherche (ANR, 212, rue de Bercy 75012 Paris) through the 2006-2008 project Adap'MC.


BUVARD online: a webtool using the R language (Shiny) to help end-users sizing their vegetative filter strips C. Catalogne, Claire Lauvernet, Nadia Carluer, Etienne Leblois

Method developed by Carluer et al. (2011) :

But handling the different modelling tools appeared to be very complex for end-users because of :

• The need of numerous and complex parameters, which are not always easily available in practical applications • The complex management of inputs/outputs between each step

Can we propose a simplified approach instead ?

1/ Describe contributive area geometry by using GIS features implemented in Hydrodem (Leblois and Creutin ,2013)

1/ Simplifying the approach by : • Generating a huge amount of scenarios for four french climatic areas.

• Summarizing the results through graphical abacus displaying efficiency for different sizes of VFS and different climatic conditions (summer or winter rains of different durations, initial moisture of soil).

2/ Implementation of an interactive web application to display abacus :

(using the R package « Shiny » features)

• No calculus is required, relevant abacus can be obtained with a parcimonious set of parameters (the most influent ones). • Readily available on the user web browser.

Introduction From a complex modeling toolkit… … to a user friendly interface presenting simplified results

BUVARD online: a webtool using the R language (Shiny) to help end-users sizing their vegetative filter strips Clotaire Catalogne, Claire Lauvernet * , Nadia Carluer, Etienne Leblois Irstea, UR MALY, 5 rue de la Doua, 69100 Villeurbanne, France (*corresponding author: claire.lauvernet@irstea.fr)

• Vegetative Filter Strips (VFS) are known as a relevant solution to limit pesticides transfers by runoff from fields to aquatic environment. • Their efficiency strongly depends on soil, agronomic and climatic conditions and they need to be optimized by considering appropriate sizing. • In this context, IRSTEA developed a complete toolkit, which allows designing site-specific VFS by simulating their efficiency to limit runoff transfers (Carluer et al., 2014). • Now, a supplementary step is carried out by implementing a web interface to help end-users to get simplified results. • Efforts to make research tools more accessible are very important to better disseminate them into the operational users community. This example illustrates a way to do so by using the scientific language R. • Contextual help is also needed to explain each parameter purpose and support their choice in difficult cases, especially here for the « Curve number », which is one of the most influent parameter.

Conclusion & perspectives

  https://cran.rproject.org/web/packages/shiny/index.html Leblois E. and Creutin J-D. (2013). Space-time simulation of intermittent rainfall with prescribed advection field : Adaptation of the turning band method. Water resources Research, 49(6):3375-3387. Carluer N., Noll D., Bernard K., Fontaine A., Lauvernet C. (2014). Dimensionner les zones tampons enherbées et boisées pour réduire le transfert hydrique des produits phytosanitaires. TSM, n°12, pp101-120. Muñoz-Carpena R., Lauvernet C., Carluer N. (2011). Development and testing of a mechanistic algorithm to calculate the influence of a shallow water table on flow dynamics through vegetative filter strips. Paper presented at the ASABE Annual International Meeting, Louisville, Kentucky (USA). USDA-SCS (1972). National Engineering Handbook, Part 630 Hydrology. Washington, D.C. The user fills in the parameter fields considering his study site characteristics At each step, (parameters choice, results interpretation) support is provided by contextual « help pages »

•

  Futur woks will improve the tool by implementation of a metamodel providing more accurate scenarios (seeLauvernet et al. oral presentation)

	Results panel
	displaying abacus
	and a numeric
	summary.
	The user can easily
	compare different
	climatic scenarios

BUVARD online BUffer strip for runoff Attenuation and pesticides Retention Design tool 5th

  international EcoSummit, 29 aout -1 septembre 2016, Montpellier

3/ Assess the VFS efficiency by using the mechanistic model VFS-mod (Munoz-Carpena et al., 2011)