Kerrie Mengersen

Christian P Robert

Judith Rousseau

Kate Lee

Reflecting about Selecting Noninformative Priors

Keywords: Mixture distribution, Non-informative prior, Bayesian analysis, Improper prior, Bayesian model choice, MCMC methods. Chapter Induced prior, Logistic model, Bayesian methods, Stability, Prior distribution Chapter 4 Noninformative prior, Mixture of distributions, Bayesian analysis, testing statistical hypotheses, Dirichlet prior, Posterior probability Chapter Noninformative prior, improper prior, Mixture of distributions, Bayesian analysis, Dirichlet prior, exchangeability, plane-sphere intersection, polar coordinates Mixture distribution, Non-informative prior, Bayesian analysis, Improper prior, Bayesian model choice, MCMC methods

the Metropolis-within-Gibbs algorithm, adaptive MCMC and the Parallel tempering algorithm. This part of the thesis is followed by the description of R package named Ultimixt which implements a generic reference Bayesian analysis of unidimensional mixtures of Gaussian distributions obtained by a location-scale parameterization of the model. This package can be applied to produce a Bayesian analysis of Gaussian mixtures with an arbitrary number of components, with no need to specify the prior distribution.

I have been very fortunate throughout my thesis to be supported by internationally renowned researchers. Now I have a great opportunity to say a big thank you to persons who contributed to my professional and human formation in the doctoral stage.

I would first like to thank my PhD supervisor, Christian P. Robert, for all the support he has given me throughout my time as a student. By his scientific personality and his rigor, he deeply influenced me and made me more and more passionate about my work. Numerous advices in my researches, always very judicious and confidence he has manifested me, encouraged me to improve myself and to progress every day.

My deep gratitude is for Hossein Bevrani who supports me as a professor for the license and Master and advisor for the PhD. I will always be grateful for his kindness, his confidence and all the time he gave me without ever counting.

During my PhD I had the chance to work with researchers like Kerrie Mengersen, Judith Rousseau and Kate Lee as my collaborators and I want to thank them for their support and for agreeing to work with me. I think in particular of working with Kate Lee which was really rewarding.

I would like to warmly thank the team of the Kurdish Institute of Paris, Campus France and French government for their financial support and for trusting me.

I am also very grateful to have been a PhD student at CEREMADE, université Paris-Dauphine, which always provided me with a lively and intellectually stimulating environment to work. I would like to thank all the faculty, Robin Ryder and particularly director Olivier Glass.

I am very grateful to my examiners, Gilles Celeux, Pierre Druilhet, Judith Rousseau and Anne Philippe, for taking the time to give me such useful feedback.

I would like to thank all the people I have met and worked with during my time at université Paris-Dauphine. It is not possible to thank everyone here but I would particularly like to mention Sofia Tsepletidou, Roxana Dumitrescu, Viviana Letizia and Nicolas Baradel ... for their scientific input, discussions, and friendship over the years.

Lastly, I would like to thank my husband Bewar, my parents and my sisters who always believed in me and always offered me unconditional love and also my friends Ayoub Moradi, Serwa Khoramdel and Fateme Mohamadipoor who have supported me all these years, even from thousands of miles away.

iii Résumé L'une des grandes applications de la statistique est la validation et la comparaison de modèles probabilistes au vu des données. Cette branche des statistiques a été développée depuis la formalisation de la fin du 19ième siècle par des pionniers comme Gosset, Pearson et Fisher. Dans le cas particulier de l'approche bayésienne, la solution à la comparaison de modèles est le facteur de Bayes, rapport des vraisemblances marginales, quelque soit le modèle évalué. Cette solution est obtenue par un raisonnement mathématique fondé sur une fonction de coût.

Ce facteur de Bayes pose cependant problème et ce pour deux raisons. D'une part, le facteur de Bayes est très peu utilisé du fait d'une forte dépendance à la loi a priori (ou de manière équivalente du fait d'une absence de calibration absolue). Néanmoins la sélection d'une loi a priori a un rôle vital dans la statistique bayésienne et par consequent l'une des difficultés avec la version traditionnelle de l'approche bayésienne est la discontinuité de l'utilisation des lois a priori impropres car ils ne sont pas justifiées dans la plupart des situations de test. La première partie de cette thèse traite d'un examen général sur les lois a priori non informatives, de leurs caractéristiques et montre la stabilité globale des distributions a posteriori en réévaluant les exemples de [Seaman III 2012].

Le second problème, indépendant, est que le facteur de Bayes est difficile à calculer à l'exception des cas les plus simples (lois conjuguées). Une branche des statistiques computationnelles s'est donc attachée à résoudre ce problème, avec des solutions empruntant à la physique statistique comme la méthode du path sampling de [Gelman 1998] et à la théorie du signal. Les solutions existantes ne sont cependant pas universelles et une réévaluation de ces méthodes suivie du développement de méthodes alternatives constitue une partie de la thèse. Nous considérons donc un nouveau paradigme pour les tests bayésiens d'hypothèses et la comparaison de modèles bayésiens en définissant une alternative à la construction traditionnelle de probabilités a posteriori qu'une hypothèse est vraie ou que les données proviennent d'un modèle spécifique. Cette méthode se fonde sur l'examen des modèles en compétition en tant que composants d'un modèle de mélange. En remplaçant le problème de test original avec une estimation qui se concentre sur le poids de probabilité d'un modèle donné dans un modèle de mélange, nous analysons la sensibilité sur la distribution a posteriori conséquente des poids pour divers modélisation préalables sur les poids et soulignons qu'un intérêt important de l'utilisation de cette perspective est que les lois a priori impropres génériques sont acceptables, tout en ne mettant pas en péril la convergence. Pour cela, les méthodes MCMC comme l'algorithme de Metropolis-Hastings et l'échantillonneur de Gibbs et des approximations de la probabilité par des méthodes empiriques sont utilisées. Une autre caractéristique de cette variante facilement mise en oeuvre est que les vitesses de convergence de la partie postérieure de la moyenne du poids et de probabilité a posteriori correspondant sont assez similaires à la solution bayésienne classique.

Dans la dernière partie de la thèse, nous sommes intéressés à la construction d'une analyse bayésienne de référence pour mélanges de gaussiennes par la création iv d'une nouvelle paramétrisation centrée sur la moyenne et la variance de ces modèles, ce qui nous permet de développer une loi a priori non-informative pour les mélanges avec un nombre arbitraire de composants. Nous démontrons que la distribution postérieure associée à ce préalable est propre et fournissons des implémentations MCMC qui exhibent l'échangeabilité attendu. L'analyse repose sur des méthodes MCMC comme l'algorithme de Metropolis-within-Gibbs, Adaptive MCMC et l'algorithme de "Parallel Tempering". Cette partie de la thèse est suivie par une package R nommée Ultimixt qui met en oeuvre une description de notre analyse bayésienne générique de mélanges de gaussiennes unidimensionnelles obtenues par une paramétrisation moyenne-variance du modèle. Ultimixt peut être appliqué à une analyse bayésienne des mélanges gaussiennes avec un nombre arbitraire de composants, sans avoir besoin de définir la loi a priori.

Mots clés: Distribution de mélange, Loi a priori non-informative, Analyse bayésienne, A priori impropre, Choix du modèle bayésien, Méthodes de MCMC. v Abstract One of the major applications of statistics is the validation and comparing probabilistic models given the data. This branch statistics has been developed since the formalization of the late 19th century by pioneers like Gosset, Pearson and Fisher. In the special case of the Bayesian approach, the comparison solution of models is the Bayes factor, ratio of marginal likelihoods, whatever the estimated model. This solution is obtained by a mathematical reasoning based on a loss function.

Despite a frequent use of Bayes factor and its equivalent, the posterior probability of models, by the Bayesian community, it is however problematic in some cases. First, this rule is highly dependent on the prior modeling even with large datasets and as the selection of a prior density has a vital role in Bayesian statistics, one of difficulties with the traditional handling of Bayesian tests is a discontinuity in the use of improper priors since they are not justified in most testing situations. The first part of this thesis deals with a general review on non-informative priors, their features and demonstrating the overall stability of posterior distributions by reassessing examples of [Seaman III 2012].

Beside that, Bayes factors are difficult to calculate except in the simplest cases (conjugate distributions). A branch of computational statistics has therefore emerged to resolve this problem with solutions borrowing from statistical physics as the path sampling method of [Gelman 1998] and from signal processing. The existing solutions are not, however, universal and a reassessment of the methods followed by alternative methods is a part of the thesis. We therefore consider a novel paradigm for Bayesian testing of hypotheses and Bayesian model comparison. The idea is to define an alternative to the traditional construction of posterior probabilities that a given hypothesis is true or that the data originates from a specific model which is based on considering the models under comparison as components of a mixture model. By replacing the original testing problem with an estimation version that focus on the probability weight of a given model within a mixture model, we analyze the sensitivity on the resulting posterior distribution of the weights for various prior modelings on the weights and stress that a major appeal in using this novel perspective is that generic improper priors are acceptable, while not putting convergence in jeopardy. MCMC methods like Metropolis-Hastings algorithm and the Gibbs sampler are used. From a computational viewpoint, another feature of this easily implemented alternative to the classical Bayesian solution is that the speeds of convergence of the posterior mean of the weight and of the corresponding posterior probability are quite similar.

In the last part of the thesis we construct a reference Bayesian analysis of mixtures of Gaussian distributions by creating a new parameterization centered on the mean and variance of those models itself. This enables us to develop a genuine non-informative prior for Gaussian mixtures with an arbitrary number of components. We demonstrate that the posterior distribution associated with this prior is almost surely proper and provide MCMC implementations that exhibit the expected component exchangeability. The analyses are based on MCMC methods as Chapter 1

General Introduction

Overview

In several areas of science, statistics is a powerful tool to analyze data both from controlled experiments such as natural sciences and from observational studies, mainly in the human sciences. Basically, a researcher expects to find methods which can provide means to judge a population from a subset of it, named a sample. Statistics has developed many different theories to be applied in different situations, and all of them have a characteristic in common and that is, given the uncertainty, they try to find the best strategy to answer scientists' queries.

In order to apply statistics to a problem, it is a common practice to start with a population or process to be studied. When the entire population is not available and only samples are studied, the inferential statistics is needed. These inferences can take the form of testing hypotheses, estimation, regression analysis, prediction and some other technics that have been recently developed such as spatial data and data mining. Furthermore, statistical inference defines random samples and describes the population being examined by a probability distribution that may have unknown parameters. Indeed, the main purpose of statistical theories is to infer properties about the probability distribution of the population of interest using observations. To do so, different paradigms of statistical inference have become established. Bayesian inference is considered as an important statistical technique especially in mathematical statistics because of its application in science beside a wide range of activities such as engineering, philosophy, medicine, sport, and law. In this thesis, we focus on the Bayesian inference. Although the original Bayesian theory was settled in the 18 th century, due to various previous computational difficulties, only in the last 30 years, the Bayesian method has grown substantially. This growth in research and applications of Bayesian methods refers to the 1980s which mostly attributed to the discovery of Markov Chain Monte Carlo methods which removed many of the computational problems.

This thesis consists of four general parts which are briefly introduced in the following sections.

Prior distribution

The Bayesian theory deals with probability statements which are conditional on the observed value and this conditional feature introduces the main difference between Bayesian and classical inferences. Despite the differences between these statistical methods, in many simple analysis we get superficially similar conclusions from the two approaches. A Bayesian statistical inference is based on a prior probability distribution of an uncertain quantity that expresses one's beliefs about this quantity before some evidence is taken into account. In other words, a prior distribution is the distribution of this uncertain quantity, named parameter, before any data is observed. Once this prior distribution is set, Bayesian inference is straightforward in terms of minimizing posterior losses, computing higher posterior density or finding the predictive distribution [Robert 2001]. But in general, a prior distribution is not easy to precisely find out and most of critics of the Bayesian analysis focussed on the choice of the prior distributions. Furthermore, different perspectives are available to choose a prior while the impact of this choice on the resulting posterior inference should not be omitted even in the case it is negligible. The main point here is about the existence of a prior or the determination of an exact or even a parametrized distribution for the prior on the parameter, which is never unique.

However, the prior plays a fundamental role in drawing Bayesian inference because of its exploitation combined with the probability distribution of data to yield the posterior distribution. Bayesian inference is fundamentally based on the posterior distribution which is used for future inference and decisions involving the parameter. [Gelman 2002] pointed out the assessment of the information that can be included in prior distributions and the properties of the resulting posterior distributions, as key issues in setting a prior. He also mentioned prior distributions as the key part of Bayesian inference and classified them to three categories: Noninformative priors, highly informative and moderately informative hierarchical prior distributions.

In fact, the existence of fairly precise scientific or lack of information about the parameter of interest leads to two classes of priors: Informative or subjective prior, and non-informative or objective priors. One method of determining the prior is a subjective evaluation of the prior probability that can be done by using past experiments of the same problem that is considered as an approximation to the real prior distribution [Robert 2001]. Another methods are based on the maximum entropy developed in [Jaynes 1980, Jaynes 1983] and as well as parametric approximations for priors resulting from restricting the choice of prior to a parametrized density and characterize the corresponding parameters using classical methods [Robert 2001]. Finally, other techniques such as empirical and hierarchical Bayes incorporate uncertainty about the prior distribution (for details see [Robert 2001]). All these methods depend on the availability of the information on the parameter of interest. In the case of limited prior input, conjugate priors can be used to construct the prior distribution, which originated in [START_REF] Raiffa | [END_REF]] and even if this choice may influence the resulting Bayesian inference, conjugate priors are not considered as part of the non-informative prior class [Robert 2001]. The most popular conjugate priors are related to the distributions associated with the exponential families which are called natural conjugate priors [Robert 2001]. This family of distributions is the only case where conjugate priors are guaranteed to exist. Despite the advantages such as being easy to deal with in both cases mathematically and computationally, the con-1.2. Prior distribution 3 jugate priors are not away from criticism. One reason is that these distributions are overly restrictive and also they are not necessarily considered as the most robust prior distributions.

The non-informative priors are requested when no information about the parameter is available. While informative priors are far from enough to allow hopes of achieving, the use of non-informative priors also underwent vary criticisms because of their influences on the relative posterior distribution. Laplace's prior is the simplest and oldest non-informative prior that is based on the principle of indifference by assigning equal probabilities to all possibilities. This prior was criticized because it results in improper resulting distributions in the case where the parameter space is infinite. This is not always a serious problem since it may lead to proper posteriors. However, the use of improper non-informative priors may also cause problem such as the marginalization paradox shown by [START_REF] Stone | UnBayesian implications of improper Bayes inference in routine statistical problems[END_REF]]. Some others are the possible inadmissibility of resulting Bayes estimators, Stein's paradox [START_REF] Syversveen | Noninformative Bayesian priors. Interpretation and problems with construction and applications[END_REF]] and in addition to the possibility of resulting improper posteriors [Kass 1996], considering equal probabilities for possible events is not coherent under partitioning as pointed out by [Robert 2001]. Another issue is the lack of invariance under the reparametrization of the parameter. The invariance of a prior is necessary when more than one inference about the parameter is needed. The best solution for obtaining invariant non-informative priors was represented by Jeffreys' distributions [Jeffreys 1939] where the information matrix of the sampling model is turned into a prior distribution. Jeffreys' prior is most often improper which means that it does not integrate to a finite value. Another method that was initially described by [Bernardo 1979] and further developed by [Berger 1979] is the reference prior. The advantages of this method compared with Jeffreys' method appear in the case of multidimensional problems [START_REF] Syversveen | Noninformative Bayesian priors. Interpretation and problems with construction and applications[END_REF]]. Some other methods have been also suggested by [START_REF] Box | [END_REF], Rissanen 2012, Welch 1963].

Since there is no best prior that one should use, research aims at acceding a prior so that posterior distribution is well behaved and proper while all available information about the parameter is taken into account. Recently, due to theoretical developments on sensitivity analysis, the dependence of posteriors on prior distributions can be checked by methods such as comparing posterior inferences under different reasonable choices of prior distribution. The first part of this thesis deals with selecting non-informative priors based on a critical review of [Seaman III 2012] and the main result of this work is to show that the Bayesian data analysis remains stable under different choices of non-informative prior distributions. A related paper was published in the journal of Applied and Computational Mathematics in July 21, 2014.

In the literature we can find a lot of theoretical and applied overviews of Bayesian statistics about the uses of non-informative priors (see [START_REF] Bernardo | [END_REF], Carlin 1996, Gelman 2013a]). A variety of methods of driving non-informative priors have been covered by [Yang 1996]. He also listed known properties of these prior distributions. Despite the wide application of non-informative priors by Bayesian community, the handling of non-informative Bayesian testing is mostly unresolved. In the following section, we briefly address hypotheses testing and related concepts.

Bayesian model choice

As mentioned at the beginning of this chapter, among many other types of statistical inference, hypotheses testing or equivalently model selection techniques are widely applied for data analysis. Statistical hypothesis tests define a procedure of controlling the probability of incorrectly deciding that a so-called null hypothesis is false.

Differences among statistical paradigms such as frequency-based or Bayesian methods are generally much more pronounced in model checking and selection than in fitting. In a Bayesian paradigm the typical method for comparing two models involves the Bayes factors or the posterior probability of the models which are based on a specification of both likelihood and prior distribution and both are compared together. Unlike standard frequency-based methods both Bayes factors and posterior probability treat the models under comparison essentially symmetrically. However, from both classical and Bayesian points of view, model selection is the problem in which we have to choose between some models on the basis of observed data but the Bayesian model comparison based on the Bayes factors does not depend on the parameters because of the integration over all parameters in each model. On the other hand, the use of Bayes factors has the advantage of automatically including a penalty for too much model structure [START_REF][END_REF].

The literature on Bayesian model choice is considerable by now and one of the earlier, reasonably thorough reviews, appears in [START_REF] Gelfand | [END_REF]]. The Bayes factors have also been the subject of much discussion in the literature in recent years and one of the comprehensive review of Bayes factors, their computation and usage in Bayesian hypothesis testing goes back to 1995 by [START_REF][END_REF]] who proposed this criterion as a solution for the comparison of models problem. However, the decision based on the Bayes factors requires a zero-one loss and [Kadane 1980] shows that these criterions are sufficient if and only if a zero-one loss obtains. Many other works on Bayesian model selection, Bayes factors and their features can be found in [START_REF] Good | Probability and the Weighting of Evidence[END_REF], Berger 1996].

Because of the difficulties caused by prior specification, the Bayesian approach to test hypotheses is not always straightforward especially in the case of an absolute lack of information. In fact, the use of non-informative prior distributions for testing hypotheses is delicate because of the sensitivity of Bayes factors to the choice of the prior. The typical strategy of using non-informative prior distributions with large variances clearly affects the Bayes factors [Robert 2001]. Furthermore, improper prior distributions result in improper prior predictive distributions and undefined Bayes factors. Among some other difficulties caused by Bayes factors that will be addressed in Chapter 4, a principal drawback from which both criterions, Bayes factors and posterior probability of models, suffer is that they can be difficult to compute. In all but the simplest cases, Bayes factors must be evaluated numerically using methods such as importance sampling, bridge sampling and reversible jump Markov Chain Monte Carlo [Green 1995]. Another method has also been recently produced by [O'Neill 2014] for computing Bayes factors that avoids the need to use reversible jump approaches. [O'Neill 2014] show that Bayes factors for the models can be expressed in terms of the posterior means of the mixture probabilities, and thus estimated from the MCMC output. In the other hand, one solution in the case that the likelihood is not available or too costly to evaluate numerically, is the approximate Bayesian computation. Some of related works can be found in [START_REF] Csilléry | Approximate Bayesian Computation (ABC) in practice[END_REF], Toni 2010, Rattan 2013] for instance. Other proposals have been made to solve particular problems with the ordinary Bayes factor such as intrinsic Bayes factors [Berger 1996] with further modifications such as the trimmed and median variants, fractional Bayes factors [O'Hagan 1995] and posterior Bayes factors [Aitkin 1991]. Consideration of Bayes factors also leads to two of the more common criteria used for model selection such as the Bayes Information Criterion (BIC) or Schwartz's criterion that provides a cursory first-order approximation to the Bayes factor [Robert 2001] and the Akaike Information Criterion (or AIC) [Akaike 1973]. A Bayesian alternative to both BIC and AIC based on the deviance has been developed by [Spiegelhalter 1998] which takes into account the prior information. Because the existing solutions are not, however, universal in the second part of this thesis our focus is towards addressing the difficulties with the traditional handling of Bayesian model selection using Bayes factors by proposing a method which goes some way to removing these complications. The key idea is to consider a mixture model whose components are the competing models of interest and the traditional method for the model choice is replaced by a kind of Bayesian estimation problem that focuses on the probability weight of the mixture model. The method includes a novel strategy of reparametrizing the competing models towards common meaning parameters in all models, that allows for using the non-informative priors at least on the common parameters. Two substantial advantages of our method are the usability of the non-informative priors for Bayesian model choice and the other is that due to the standard MCMC algorithms, the Bayesian estimation of the model is straightforward and there is no need to compute the marginal likelihoods. A related paper was submitted for publication.

The third part of this thesis focuses on the parametrization of the mixture distributions. In the following we briefly introduce the motivation of this work.

Mixture distributions

The earliest study about the mixture models was done by [START_REF] Pearson | Contributions to the mathematical theory of evolution[END_REF]] who investigated the estimation of parameters in the finite mixture model by the use of the method of moments. In 1894, [START_REF] Pearson | Contributions to the mathematical theory of evolution[END_REF]] studied the dissection of asymptotic and symmetric frequency curves into two components of normal distributions. Many other papers have appeared related to the problem of statistical inference about the parameters and probabilistic properties of these densities. Since this early work, finite mixture models have been widely used in many disciplines and there is a large body of literature on these distributions. For example in biology it is often desired to measure certain characteristics in natural populations of some particular species when the distribution of such characteristics may vary markedly with age of the individuals. Since age is difficult to ascertain in samples from populations, the biologist is dealing with a mixture of distributions and the mixing in this case is done over a parameter depending on the unobservable variate, age. Some other applications can be found in astronomy, ecology, genetics and so on due to the feature that they are easily applied to the data set in which two or more subpopulations are mixed together. In statistical applications, the mixture of densities can be used to approximate some parameters associated with a density.

The finite mixture models have also enjoyed intensive attentions over the recent years from both practical and theoretical viewpoints due to their flexibility in modeling. Some basic properties of mixtures were studied by [Robbins 1948] and [Robbins 1961] initiated the study of identifiability problem. Despite the popularity of mixtures, model estimation can be difficult when the number of components is unknown. In 1966, [Hasselblad 1966] first considered the estimation problem of mixtures by the method of maximum likelihood. [Rolph 1968] first considered Bayes estimation of the mixture parameters in the special case where the observations from the mixture population are restricted to the positive integers. In the framework of the Bayesian approach, one needs to assume that a prior distribution on component parameters is available. As summarized in [Frühwirth-Schnatter 2006], there are two main reasons why people may be interested in using the Bayesian method in finite mixture models. Firstly, including a suitable prior distribution or the parameters in the framework of the Bayesian approach may avoid spurious modes when maximizing the log-likelihood function. Secondly, when the posterior distribution for the unknown parameters is available, the Bayesian method can yield valid inference without relying on asymptotic normality. This is an advantage of the Bayesian method for estimating the parameters of a mixture distribution without the need of sample sizes very large. As mentioned before, the use of the conjugate prior produces the posterior distribution that may belong to the tractable distribution family. However, because of the complexity of mixtures, it is impossible to find a conjugate prior for the component parameters. While the posterior distributions derived from the mixture models are non standard, MCMC methods are used to generate samples from these complex distributions [START_REF] Editeurs | [END_REF], Frühwirth-Schnatter 2006]. Because the main idea of Bayesian estimation using MCMC methods followed by realizing a mixture model is considered as a special case of incomplete data problem with the missing component indicator variables, the problem with conjugate priors no longer poses serious obstacles to the application of Bayesian method.

In fact, the Bayesian estimators in mixture models are always well defined as long as priors are proper. Furthermore, the unidentifiability may be resolved by well defining the parameter space or using informative priors on parameters. However, in the case where no information is available for the component parameters, the choice of the prior is more delicate. [START_REF] Editeurs | [END_REF]] demonstrates that specifying improper prior to the component parameters results in improper posterior distribution that prohibits this kind of prior to be used for mixtures. In addition, non-informative priors assigned to the parameter of a specific component can also lead to identifiability problems. Because if each component has its own prior parameters and few observations are allocated to this component, there will be no information at all to estimate the parameter and in the case of Gibbs sampling, the sampler gets trapped in a local mode corresponding to this component.

This problem of non-identifiability in the posterior distribution can also be due to an overfitting phenomenon. Basically this happens when some components have weights equal to zero or merged together [Frühwirth-Schnatter 2006]. A full discussion about how over fitted mixtures behave can be found in [START_REF] Rousseau | [END_REF]] who proved that the posterior behavior of overfitted mixtures generally depends on both the choice of the prior on the weights and the number of free parameters.

[van Havre 2015] treated the issues such as non-identifiability due to overfitting, label switching and also the problem of lack of mixing caused by applying standard MCMC sampling techniques when the posterior contains multiple well separated modes.

Given the difficulty with non-informative priors, one solution is to use proper priors with the prior parameters chosen such that the prior is suitably weakly informative priors [START_REF] Richardson | [END_REF]]. This method is not always applicable because of the problem of multiple prior specifications. Another method proposed by [Diebolt 1994] is to use an improper prior under the condition of forcing each component to always have a minimal number of data points assigned to it. A related work has been recently developed by [Stoneking 2014] which does not result in any data dependence of the priors.

In the third part of this thesis we define a novel reparametrisation for the mixture of distributions based on the mean and standard deviation of the mixture itself, namely global parameters. The main feature of our method is that the noninformative prior distribution can be used on the global parameters of the mixture while the resulting posterior distribution is proper. A related paper was submitted for publication.

The reparametrized mixture model will be fitted with our R package named Ultimixt. Ultimixt provides the functionality for estimating reparametrized Gaussian mixture models with MCMC methods. The last part of this thesis pertains to the description of the implementation and the functions of Ultimixt. This package can accurately compute the posterior estimate of the parameters of reparametrized univariate Gaussian mixture distribution beside having the ability of graphically summarizing the posterior results.

Introduction

The choice of a particular prior for the Bayesian analysis of a statistical model is often seen more as an art than as a science. When the prior cannot be derived from the available information, it is generally constructed as a noninformative prior. This derivation is mostly mathematical and, even though the corresponding posterior distribution has to be proper and hence constitutes a correct probability density, it nonetheless leaves the door open to criticism. The focus of this note is the paper by [Seaman III 2012], where the authors consider using a particular noninformative distribution as a problem in itself, often bypassed by users of these priors: "if parameters with diffuse proper priors are subsequently transformed, the resulting induced priors can, of course, be far from diffuse, possibly resulting in unintended influence on the posterior of the transformed parameters" (p.77). Using the inexact argument that most problems rely on MCMC methods and hence require proper priors, the authors restrict the focus to those priors.

In their critical study, [Seaman III 2012] investigate the negative side effects of some specific prior choices related with specific examples. Our note aims at re-examining their investigation and at providing a more balanced discussion on these side effects. We first stress that a prior is considered as informative by [Seaman III 2012] "to the degree it renders some values of the quantity of interest more likely than others" (p.77), and with this definition, when comparing two priors, the prior that is more informative is deemed preferable. In contrast with this definition, we consider that an informative prior expresses specific, definite (prior) information about the parameter, providing quantitative information that is crucial to the estimation of a model through restrictions on the prior distribution [Robert 2007]. However, in most practical cases, a model parameter has no substance per se but instead calibrates the probability law of the random phenomenon observed therein. The prior is thus a tool that summarizes the information available on this phenomenon, as well as the uncertainty within the Bayesian structure. Many discussions can be found in the literature on how appropriate choices between the prior distributions can be decided. In this case, robustness considerations also have an important role to play [Lopes 2011, Stojanovski 2011]. This point of view will be obvious in this note as, e.g., in processing a logistic model in the following section. Within the sole setting of the examples first processed in [Seaman III 2012], we do exhibit a greater stability in the posterior distributions through various noninformative priors.

The plan of the note is as follows: we first provide a brief review of noninformative priors in Section 4.2. In Section 2.3, we propose a Bayesian analysis of a logistic model (Seaman III et al.'s (2012) first example) by choosing the normal distribution N (0, σ 2 ) as the regression coefficient prior. We then compare it with a g-prior, as well as flat and Jeffreys' priors, concluding to the stability of our results. The next sections cover the second to fourth examples of [Seaman III 2012], modeling covariance matrices, treatment effect in biomedical studies, and a multinomial distribution. When modeling covariance matrices, we compare two default priors for the standard deviations of the model coefficients. In the multinomial setting, we discuss the hyperparameters of a Dirichlet prior. Finally, we conclude with the argument that the use of noninformative priors is reasonable within a fair range and that they provide efficient Bayesian estimations when the information about the parameter is vague or very poor.

Noninformative priors

As mentioned above, when prior information is unavailable and if we stick to Bayesian analysis, we need to resort to one of the so-called noninformative priors. Since we aim at a prior with minimal impact on the final inference, we define a noninformative prior as a statistical distribution that expresses vague or general information about the parameter in which we are interested. In constructive terms, the first rule for determining a noninformative prior is the principle of indifference, using uniform distributions which assign equal probabilities to all possibilities [START_REF] Laplace | Théorie Analytique des Probabilités[END_REF]]. This distribution is however not invariant under reparametrization ,(see [Berger 1980, Robert 2007] for references). If the problem does not allow for an invariance structure, Jeffreys' priors [Jeffreys 1939], then reference priors, exploit the probabilistic structure of the problem under study in a more formalized way. Other methods have been advanced, like the little-known data-translated likelihood of [START_REF] Box | [END_REF]], maxent priors [Jaynes 2003], minimum description length priors [Rissanen 2012] and probability matching priors [Welch 1963]. [START_REF] Bernardo | [END_REF]] envision noninformative priors as a mere mathematical tool, while accepting their feature of minimizing the impact of the prior selection on inference: "Put bluntly, data cannot ever speak entirely for themselves, every prior specification has some informative posterior or predictive implications and vague is itself much too vague an idea to be useful. There is no "objective" prior that represents ignorance" (p.298). There is little to object against this quote since, indeed, prior distributions can never be quantified or elicited exactly, especially when no information is available on those parameters. Hence, the concept of "true" prior is meaningless and the quantification of prior beliefs operates under uncertainty. As stressed by [Berger 1994], noninformative priors enjoy the advantage that they can be considered to provide robust solutions to relevant problems even though "the user of these priors should be concerned with robustness with respect to the class of reasonable noninformative priors" (p.59).

Example 1: Bayesian analysis of the logistic model

The first example in [Seaman III 2012] is a standard logistic regression modeling the probability of coronary heart disease as dependent on the age x by

ρ(x) = exp(α + βx) 1 + exp(α + βx) . (2.1)
First we recall the original discussion in [Seaman III 2012] and then run our own analysis by selecting some normal priors as well as the g-prior, the flat prior and Jeffreys' prior.

Seaman et al.'s (2012) analysis

For both parameters of the model (2.1), [Seaman III 2012] chose a normal prior N (0, σ 2 ). A first surprising feature in this choice is to opt for an identical prior on both intercept and slope coefficients, instead of, e.g., a g-prior (discussed in the following) that would rescale each coefficient according to the variation of the corresponding covariate. Indeed, since x corresponds to age, the second term βx in the regression varies 50 times more than the intercept. When plotting logistic cdf's induced by a few thousands simulations from the prior, those cumulative functions mostly end up as constant functions with the extreme values 0 and 1. This behavior is obviously not particularly realistic since the predicted phenomenon is the occurrence of coronary heart disease. Under this minimal amount of information, the prior is thus using the wrong scale: the simulated cdfs should have a reasonable behavior over the range (20, 100) of the covariate x. For instance, it should focus on a -5 log-odds ratio at age 20 and a +5 log-odds ratio at 100, leading to the comparison pictured in Figure 2.1 (left versus right). Furthermore, the fact that the coefficient of x may be negative also bypasses a basic item of information about the model and answers the later self-criticism in [Seaman III 2012] that the prior probability that the ED50 is negative is 0.5. Using instead a flat prior would answer the authors' criticisms about the prior behavior, as we now demonstrate. We stress that [Seaman III 2012] produce no further justification for the choice of the prior variance σ 2 = 25 2 , other than there is no information about the model parameters. This is a completely arbitrary choice of prior, arbitrariness that does have a considerable impact on the resulting inference, as already discussed. [Seaman III 2012] further criticized the chosen prior by comparing both posterior mode and posterior mean derived from the normal prior assumption with the MLE. If the MLE is the golden standard there then one may wonder about the relevance of a Bayesian analysis! When the sample size N gets large, most simple Bayesian analyses based on noninformative prior distributions give results similar to standard non-Bayesian approaches [Gelman 2013a]. For instance, we can often interpret classical point estimates as exact or approximate posterior summaries based on some implicit full probability model. Therefore, as N increases, the influence of the prior on posterior inferences decreases and, when N goes to infinity, most priors lead to the same inference. However, for smaller sample sizes, it is inappropriate to summarize inference about the parameter by one value like the mode or the mean, especially when the posterior distribution of the parameter is more variable or even asymmetric.

The dataset used here to infer on (α, β) is the Swiss banknote benchmark (available in R). The response variable y indicates the state of the banknote, i.e. whether the bank note is genuine or counterfeit. The explanatory variable is the bill length. This data yields the maximum likelihood estimates α = 233.26 and β = -1.09. To check the impact of the normal prior variance, we used a random walk Metropolis-Hastings algorithm as in [START_REF][END_REF]] and derived the estimators reproduced in Table 2.1. We can spot definitive changes in the results that are caused by moves in the coefficient σ, hence concluding to the clear sensitivity of the posterior to the choice of hyperparameter σ (see also Figure 2.2).

Larger classes of priors

Normal priors are well-know for their lack of robustness (see e.g. [Berger 1994]) and the previous section demonstrates the long-term impact of σ. However, we can limit variations in the posteriors, using the g-priors of [Zellner 1986],

α, β | X ∼ N 2 (0, g(X T X) -1 ).
(2.2)

where the prior variance-covariance matrix is a scalar multiple of the information matrix for the linear regression. This coefficient g plays a decisive role in the analysis, however large values of g imply a more diffuse prior and, as shown e.g. in [START_REF][END_REF]], if the value of g is large enough, the Bayes estimate stabilizes. We will select g as equal to the sample size 200, following [Liang 2008], as it means that the amount of information about the parameter is equal to the amount of information contained in one single observation.

A second reference prior is the flat prior π(α, β) = 1. And Jeffreys' prior constitutes our third prior as in [START_REF][END_REF]]. In the logistic case, Fisher's information matrix is I(α, β, X) = X T W X, where X = {x ir } is the design matrix,

W = diag{m i π i (1 -π i )
} and m i is the binomial index for the ith count [Firth 1993]. This leads to Jeffreys' prior {det(I(α, β, X))}

1 2 , proportional to   n � i=1 exp(α + βx i ) {1 + exp(α + βx i )} 2 n � i=1 x 2 i exp(α + βx i ) {1 + exp(α + βx i )} 2 - � n � i=1 x i exp(α + βx i ) {1 + exp(α + βx i )} 2 � 2   1 2
This is a nonstandard distribution on (α, β) but it can be easily approximated by a Metropolis-Hastings algorithm whose proposal is the normal Fisher approximation of the likelihood, as in [START_REF][END_REF]]. Bayesian estimates of the regression coefficients associated with the above three noninformative priors are summarized in Table 2.2. Those estimates vary quite moderately from one choice to the next, as well as relatively to the MLEs and to the results shown in Table 2.1 when σ = 900. Figure 2.3 is even more definitive about this stability of Bayesian inferences under different noninformative prior choices.

Example 2: Modeling covariance matrices

The second choice of prior criticized by [Seaman III 2012], was proposed by [Barnard 2000] for the modeling of covariance matrices. However the paper falls short of demonstrating a clear impact of this prior modeling on posterior inference. Furthermore the adopted solution of using another proper prior resulting in a "wider" dispersion requires a prior knowledge of how wide is wide enough. We thus run Bayesian analyses considering prior beliefs specified by both [Seaman III 2012] and [Barnard 2000].

Setting

The multivariate regression model of [Barnard 2000] is

Y j | X j , β j , τ j ∼ N (X j β j , τ 2 j I n j ), j = 1, 2, . . . , m. (2.3)
where Y j is a vector of n j dependent variables, X j is an n j ×k matrix of covariate variables, and β j is a k-dimensional parameter vector. For this model, [Barnard 2000] considered an iid normal distribution as the prior

β j ∼ N ( β, Σ)
conditional on β, Σ where β, τ 2 j for j = 1, 2, . . . , m are independent and follow a normal and inverse-gamma priors, respectively. Assuming that β, τ 2 j 's and Σ are a priori independent, [Barnard 2000] firstly provide a full discussion on how to choose a prior for Σ because it determines the nature of the shrinkage of the posterior of the individual β j is towards a common target. The covariance matrix Σ is defined as a diagonal matrix with diagonal elements S, multiplied by a k × k correlation matrix R,

Σ = diag(S)Rdiag(S) .
Note that S is the k × 1 vector of standard deviations of the β j s, (S 1 , . . . , S k ). [Barnard 2000] propose lognormal distributions as priors on S j . The correlation matrix could have (1) a joint uniform prior p(R) ∝ 1, or (2) a marginal prior obtained from the inverse-Wishart distribution for Σ which means p(R) is derived from the integral over S 1 , . . . , S k of a standard inverse-Wishart distribution. In the second case, all the marginal densities for r ij are uniform when i � = j [ Barnard 2000].

Considering the case of a single regressor, i.e. k = 2, [Seaman III 2012] chose a different prior structure, with a flat prior on the correlations and a lognormal prior with means 1 and -1, and standard deviations 1 and 0.5 on the standard deviations of the intercept and slope, respectively. Simulating from this prior, they concluded at a high concentration near zero. They then suggested that the lognormal distribution should be replaced by a gamma distribution G(4, 1) as it implies a more diffuse prior. The main question here is whether or not the induced prior is more diffuse should make us prefer gamma to lognormal as a prior for S j , as discussed below.

Prior beliefs

First, Barnard et al.'s (2000) basic modeling intuition is "that each regression is a particular instance of the same type of relationship" (p.1292). This means an exchangeable prior belief on the regression parameters. As an example, they suppose that m regressions are similar models where each regression corresponds to a different firm in the same industry branch. Exploiting this assumption, when β j has a normal prior like β ij ∼ N ( βi , σ 2 i ), j = 1, 2, . . . , m, the standard deviation of β ij (S i = σ i ) should be small as well so "that the coefficient for the ith explanatory variable is similar in the different regressions" (p.1293). In other words, S i concentrated on small values implies little variation in the ith coefficient. Toward this goal, [Barnard 2000] chose a prior concentrated close to zero for the standard deviation of the slope so that the posterior of this coefficient would be shrunken together across the regressions. Based on this basic idea and taking tight priors on Σ for β j , j = 1, . . . , m, they investigated the shrinkage of the posterior on β j as well as the degree of similarity of the slopes. Their analysis showed that a standard deviation prior that is more concentrated on small values results in substantial shrinkage in the coefficients relative to other prior choices.

Consider for instance the variation between the choices of lognormal and gamma distributions as priors of S 2 , standard deviation of the regression slope. In this case, most of the mass of the lognormal prior is concentrated on values close to zero whereas the gamma prior is more diffuse. The 10, 50, 90 percentiles of LN (-1, 0.5) and G(4, 1) are 0.19, 0.37, 0.7 and 1.74, 3.67, 6.68, respectively. Thus, choosing LN (-1, 0.5) as the prior of S 2 is equivalent to believe that values of β 2 in the m regressions are much closer together than the situation where we assume S 2 ∼ G(4, 1). To assess the difference between both prior choices on S 2 and their impact on the degree of similarity of the regression coefficients, we resort to a simulated example, similar to [Barnard 2000], except that m = 4 and n j = 36.

The explanatory variables are simulated standard normal variates. We also take τ j ∼ IG(3, 1) and β ∼ N (0, 1000I). The prior for Σ is such that π(R) ∝ 1 and we run Seaman et al.'s (2012) analyses under S 2 ∼ LN (-1, 0.5) and S 2 ∼ G(4, 1).

S i ∼ LN (-

Comparison of posterior outputs

As seen in Tables 2.3 and 2.4, respectively. The differences between the regression estimates are quite limited from one prior to the next, while the estimates of the standard deviations vary much more. In the lognormal case, the posterior of S i is concentrated on smaller values relative to the gamma prior. Figure 2.5 displays the posterior distributions of those parameters. The impact of the prior choice is quite clear on the standard deviations. Therefore, since the posteriors of both intercepts and slopes for all four regressions are centered in (16.5, 17) and (-10, -9), respectively, we can conclude at the stability of Bayesian inferences on β j when selecting two different prior distributions on S j . That the posteriors on the S i 's differ is in fine natural since those are hyperparameters that are poorly informed by the data, thus reflecting more the modeling choices of the experimenter.

Examples 3 and 4: Prior choices for a proportion and the multinomial coefficients

This section considers more briefly the third and fourth examples of [Seaman III 2012].

The third example relates to a treatment effect analyzed by [Cowles 2002] and the fourth one covers a standard multinomial setting.

Proportion of treatment effect captured

In [Cowles 2002] two models are compared for surrogate endpoints, using a link function g that either includes the surrogate marker or not. The quantity of interest is a proportion of treatment effect captured: it is defined as

PTE ≡ 1 -β 1 /β R,1 ,
where β 1 , β R,1 are the coefficients of an indicator variable for treatment in the first and second regression models under comparison, respectively. [Seaman III 2012] restricted this proportion to the interval (0, 1) and under this assumption they proposed to use a generalized beta distribution on β 1 , β R,1 so that PTE stayed within (0, 1). We find this example most intriguing in that, even if PTE could be turned into a meaningful quantity (given that it depends on parameters from different models), the criticism that it may take values outside (0, 1) is rather dead-born since it suffices to impose a joint prior that ensures the ratio stays within (0, 1). This actually is the solution eventually proposed by the authors. If we have prior beliefs about the parameter space (which depends on β 1 /β R,1 in this example) the prior specified on the quantity of interest should integrate these beliefs. In the current setting, there is seemingly no prior information about (β 1 , β R,1 ) and hence imposing a prior restriction to (0, 1) is not a logical specification. For instance, using normal priors on β 1 and β R,1 lead to a Cauchy prior on β 1 /β R,1 , which support is not limited to (0, 1). We will not discuss this rather artificial example any further.

Multinomial model and evenness index

The final example in [Seaman III 2012] deals with a measure called evenness index

H(θ) = - � θ i log(θ i ) � log(K)
that is a function of a vector θ of proportions θ i , i = 1, . . . , K. The authors assume a Dirichlet prior on θ with hyperparameters first equal to 1 then to 0.25. For the transform H(θ), Figure 2.6 shows that the first prior concentrates on (0.5, 1) whereas the second does not. Since there is nothing special about the uniform prior, re-running the evaluation with the Jeffreys prior reduces this feature, which anyway is a characteristic of the prior distribution, not of a posterior distribution that would account for the data. The authors actually propose to use the Dir(1/4, 1/4, . . . , 1/4) prior, presumably on the basis that the induced prior on the evenness is then centered close to 0.5. If we consider the more generic Dir(γ 1 , . . . , γ K ) prior, we can investigate the impact of the γ i 's when they move from 0.1 to 1. In Figure 2.6, the induced priors on H(θ) indeed show a decreasing concentration of the posterior on (0.5, 1) as γ i decreases towards zero. To further the comparison, we generated datasets of size N = 50, 100, 250, 1000, 10, 000. Figure 2.7 shows the posteriors associated with each of the four Dirichlet priors for these samples, including modes that are all close to 0.4 when N = 10 4 . Even for moderate sample sizes like 50, the induced posteriors are almost similar. When the sample size is 50, Table 2.5 shows there is some degree of variation between the posterior means, even though, as expected, this difference vanishes when the sample size increases.

Note that, while Dirichlet distributions are conjugate priors, hence potentially lacking in robustness, Jeffreys's prior is a special case corresponding to γ i = 1/K (here K is equal to 8). Figure 2.8 reproduces the transform of Jeffreys' prior for the evenness index (left) and the induced posterior densities for the same values of N . Since it is a special case of the above, the same features appear. A potential alternative we did not explore is to set a non-informative prior on the hyperparameters of the Dirichlet distribution.

Conclusion

In this note, we have reassessed the examples supporting the critical review of [Seaman III 2012], mostly showing that off-the-shelf noninformative priors are not suffering from the shortcomings pointed out by those authors. Indeed, according to the outcomes produced therein, those noninformative priors result in stable posterior inferences and reasonable Bayesian estimations for the parameters at hand. We thus consider the level of criticism found in the original paper rather unfounded, as it either relies on a highly specific choice of a proper prior distribution or on bypassing basic prior information later used for criticism. The paper of [Seaman III 2012] available on some quantities of interest. Our sole point of contention is the repeated and recommended reference to MLE, if only because it implies assessing or building the prior from the data. The most specific (if related to the above) recommendation is to use conditional mean priors as exposed by [Christensen 2011]. For instance, in the first (logistic) example, this meant putting a prior on the cdfs at age 40 and age 60. The authors picked a uniform in both cases, which sounds inconsistent with the presupposed shape of the probability function.

In conclusion, we find there is nothing pathologically wrong with either the paper of [Seaman III 2012] or the use of "noninformative" priors! Looking at induced priors on more intuitive transforms of the original parameters is a commendable suggestion, provided some intuition or prior information is already available on those. Using a collection of priors including reference or invariant priors helps as well towards building a feeling about the appropriate choice or range of priors and looking at the dataset induced by simulating from the corresponding predictive cannot hurt.

Chapter 3

Supplementary material: Reflecting about Selecting Noninformative Priors

This chapter contains the statistical tools, computational details and some more data analyses related to the examples studied in the chapter 2.

Example 1

The first example of Chapter 2 is about the Bayesian analysis of the logistic model. The non standard posterior distributions resulted by assigning different non-informative priors to the parameters of the model are given as follows � for a flat prior π(α, β) = 1:

f (α, β|ρ, x) = exp( � n i=1 ρ i (α+βx i )) / � n i=1 (1+exp(α+βx i ))
� for g-prior α, β|X ∼ N (0, g(X T X) -1 ):

f (α, β|ρ, x) = |X T X| 1/2 exp � -g/2( α β ) T X T X( α β )+ � n i=1 ρ i (α+βx i ) � /2π √ g � n i=1 (1+exp(α+βx i ))
� and for the Jeffrey's prior, the log-likelihood of the logistic model is given by

�(α, β) = n � i=1 (ρ i (α + βxi) -ln(1 + exp(α + βx i )))
The second derivate of � with respect to α and β is

∂ 2 �(α, β) ∂α 2 = - n � i=1 exp(α+βx i ) /(1+exp(α+βx i )) 2 ∂ 2 �(α, β) ∂β 2 = - n � i=1 x 2 i exp(α+βx i ) /(1+exp(α+βx i )) 2 ∂ 2 �(α, β) ∂α∂β = - n � i=1 x i exp(α+βx i ) /(1+exp(α+βx i )) 2
Chapter 3. Supplementary material: Reflecting about Selecting Noninformative Priors

and the matrix of Fisher information can be written as following:

I(α, β, x) = � - � n i=1 exp(α+βx i ) (1+exp(α+βx i )) 2 - � n i=1 x i exp(α+βx i ) (1+exp(α+βx i )) 2 - � n i=1 x i exp(α+βx i ) (1+exp(α+βx i )) 2 - � n i=1 x 2 i exp(α+βx i ) (1+exp(α+βx i )) 2 � �
The invariant Jeffreys prior computed from � det(I(α, β, x)) yields the following posterior distribution for α and β

f (α, β|ρ, x) = � � � � n � i=1 exp(α + βx i ) (1 + exp(α + βx i )) 2 n � i=1 x 2 i exp(α + βx i ) (1 + exp(α + βx i )) 2 -{ n � i=1 x i exp(α + βx i ) (1 + exp(α + βx i )) 2 } 2 × exp ( � n i=1 ρ i (α + βx i )) � n i=1 (1 + exp (α + βx i ))
.

We can sample from the posterior distributions above using the Metropolis-Hastings algorithm for each prior specification in which the proposal distribution is a random walk multivariate normal distribution based on the maximum likelihood estimate as starting value and the asymptotic covariance matrix of the maximum likelihood estimate as the covariance matrix of the proposal. The implementation in R can be found in [Kamary 2016a].

We run the Metropolis-Hastings algorithm with 10 4 iterations for the bank dataset by considering the bill length as the explanatory variable, and we test three different proposal scales, τ = 0.1, 1, 5. As shown in Figures 3. 1, 3.2 and 3.3, for all values of τ , the chains simulated for α and β are convergent to the target distribution and able to move around the normal range with decreasing autocorrelations. However, in the case where τ = 5, the acceptance rate is low and the histograms of the output are far from the target distribution even after 10 4 iterations. The autocorrelation graph for τ = 1 decreases quicker than the cases where τ = 0.1, 5. By comparing the raw sequences and the autocorrelation graphs provided by three algorithms above and also the corresponding acceptance rates, the best mixing behavior is related to τ = 1.

By comparing the plots shown in Figures 3.1, 3.2 and 3.3, despite the fact that three different non informative priors were assigned to the parameters of the logistic model, there is no visible difference between the posterior draws.

Example 2

Bayesian inference of multivariate regression model 2.3 using the prior modeling defined in section 2.4 derives the following joint posterior probability

(β j , τ 2 j , β, S, R|Y j , X j ) ∝ �(β j , τ 2 j , β|Y j , X j ) × π(β j | β, S, R)π(τ 2 j )π( β)π(S, R) ∝ 1 /(τ 2 j ) n j/2-a-1 exp � -(Y j -X j β j ) T (Y j -X j β j ) /2τ 2 j � × |diag(S)Rdiag(S)| -1/2 exp � -(β j -β) T (diag(S)Rdiag(S)) -1 (β j -β) /2 � × exp( -b /τ 2 j -βT β/2000)π(S)
Chapter 3. Supplementary material: Reflecting about Selecting Noninformative Priors that is resulted under the assumption of independence between the parameter S and R and π(R) ∝ 1. The conditional posterior distribution of the parameter β j can be obtained as following

π(β j |τ 2 j , β, S, R, Y j , X j ) ∝ exp � -(Y j -X j β j ) T (Y j -X j β j ) /2τ 2 j -(β j -β) T (diag(S)Rdiag(S)) -1 (β j -β) /2 � ∝ exp � -(Y j -X j β j ) T (τ 2 j I k ) -1 (Y j -X j β j )-(β j -β) T (diag(S)Rdiag(S)) -1 (β j -β) /2
� where I k is the identity matrix of size k. If we replace diag(S)Rdiag(S) by Σ, dropping multiplicative terms that do not involve β j gives

π(β j |τ 2 j , β, S, R, Y j , X j ) ∝ exp(-1 /2 � X j β j -Y j β j - β � T � τ 2 j I k 0 0 Σ � -1 � X j β j -Y j β j - β � ) or π(β j |τ 2 j , β, S, R, Y j , X j ) ∝ exp(-1 /2( � X j I k � β j - � Y j β � ) T � τ 2 j I k 0 0 Σ � -1 ( � X j I k � β j - � Y j β � ))
As shown in [Christensen 2011], if we define

βj = � τ -2 j X T j X j + Σ -1 � -1 (τ -2 j X T j Y j + Σ -1 β)
we can rewrite the posterior density as

3.2. Example 2 29 π(β j |τ 2 j , β, S, R, Y j , X j ) ∝ exp � -1 /2(β j -βj ) T � τ -2 j X T j X j + Σ -1 � (β j -βj )
� which implies a multivariate Gaussian distribution as following

β j |τ 2 j , β, S, R, Y j , X j ∼ N � βj , � τ -2 j X T j X j + Σ -1 � -1 � .
The posterior density of the parameter τ 2 j given β j , β, S, R, Y j , X j is inverse gamma distribution

τ 2 j |β j , β, S, R, Y j , X j ∼ IG � n j +2a /2, (Y j -X j β j ) T (Y j -X j β j )+2b /2 �
when the prior distribution is supposed to be IG(a, b). To obtain the full conditional for β, we can write

β|β j , τ 2 j , S, R, Y j , X j ∝ exp � -(β j -β) T (diag(S)Rdiag(S)) -1 (β j -β) /2 -βT β/2 � ∝ exp(-1 /2 � β -β j β -0 � T � Σ 0 0 1000I k � -1 � β -β j β -0 � )
which implies a multivariate Gaussian distribution with the following form

β|β j , τ 2 j , S, R, Y j , X j ∼ N � β 0 , (Σ -1 + 1 /1000I k ) -1 � (3.1)
where β 0 = (Σ -1 + 1 /1000I k ) -1 Σ -1 β j . For the correlation matrix R, the full conditional density will therefore be

R|β j , β, τ 2 j , S, Y j , X j ∝ |diag(S)Rdiag(S)| -1/2 exp � -(β j -β) T (diag(S)Rdiag(S)) -1 (β j -β) /2 � .
When the prior of the standard deviations, S, is log normal LN (µ, σ), the conditional posterior is given by

S|β j , β, τ 2 j , R, Y j , X j ∝ |diag(S)Rdiag(S)| -1/2 exp � -(β j -β) T (diag(S)Rdiag(S)) -1 (β j -β) /2 � × k � j=1 1 /s j exp(-(ln(s j )-µ) 2 /2σ 2 )
and in the case where gamma G(δ, ζ) prior is placed on the elements of S, we will have

S|β j , β, τ 2 j , R, Y j , X j ∝ |diag(S)Rdiag(S)| -1/2 exp � -(β j -β) T (diag(S)Rdiag(S)) -1 (β j -β) /2 � × k � j=1 s δ-1 j exp(-s j/ζ)
Chapter 3. Supplementary material: Reflecting about Selecting Noninformative Priors The analyses of this example in Chapter 2 are based on the Metropolis-within-Gibbs algorithm in which the correlations r ij ; i � = j are independently simulated from uniform proposal distributions. For both prior specifications of S, log normal and gamma distributions, the implementation in R can be seen in [Kamary 2016a].

In the algorithm, the parameter β j and Σ = diag(S)Rdiag(S) are initialized by the maximum likelihood estimate and the asymptotic covariance matrix of the maximum likelihood estimate and the other parameters are started from a random value simulated from their prior distributions. The proposal distribution of logarithm function of S is a random walk multivariate normal distribution with the possibility of calibrating the proposal scale. As discussed in Chapter 2, replacing log normal distribution assigned to the hyper parameter S by gamma distribution does not influence the conditional posterior distribution of the regression intercept and slope. Here, we deal with the convergence of the simulated samples obtained by implementing the Metropolis-within-Gibbs algorithm and also the impact of both prior choices on Bayesian inference of the other parameters. To do so, we simulate n 1 = 36 data points from normal regression model with an explanatory variable and run the program over 10 4 iterations when the scale of the proposal distribution of s j is equal to 1. the output of the simulation run in the case of log normal prior) the chains simulated for β1j have a lot of density over a narrow interval near the true values while the distribution of β1j spreads out over a wide range in the case of gamma prior. This impact of the prior choices for the standard deviations on the posterior distribution of β becomes more visible when comparing the empirical density plotted in Figure 3.6. In both cases, the empirical densities are centered over the true values of the parameter β while the dispersion of two cases is impacted by the change in the prior choice for s j 's. The main reason for this effect is that as shown in (3.1), the conditional posterior distribution of β depends on Σ and so on S. A stretched or squeezed prior choice allocated to s j 's influences the posterior results of standard deviation S and therefore impacts those of β.

Introduction

While a if not the central problem of statistical inference and a dramatically differentiating feature between classical and Bayesian paradigms [Neyman 1933[START_REF] Berger | [END_REF][START_REF] Casella | [END_REF], Gigerenzer 1991, Berger 2003a[START_REF] Mayo | Frequentist statistics as a theory of inductive inference[END_REF], Gelman 2008], the handling of hypothesis testing by Bayesian theory is wide open to controversy and divergent opinions, even within the Bayesian community [Jeffreys 1939, Bernardo 1980, Berger 1985, Aitkin 1991, Berger 1992, De Santis 1997, Bayarri 2007, Christensen 2011[START_REF][END_REF], Gelman 2013a, Robert 2014]. In particular, the handling of the noninformative Bayesian testing case is mostly unresolved and has produced much debate, witness the specific case of the Lindley or Jeffreys-Lindley paradox [Lindley 1957, Shafer 1982, DeGroot 1982, Robert 1993, Lad 2003, Spanos 2013, Sprenger 2013, Robert 2014].

Bayesian model selection is understood here as the comparison of several potential statistical models towards the selection of the model that fits the current data the "best". For instance, [Christensen 2011] consider this is a decision issue that pertains to testing, while [Robert 2001] expressed it as a model index estimation setting and [Gelman 2013a] do not agree about the decisional aspect. A mostly accepted perspective is however that Bayesian model selection does not primarily seek to identify which model is "true" (if any), but rather to indicate which model fits the data better given all the available information. As discussed in the Bayesian literature (see, e.g. [Berger 1992[START_REF] Madigan | [END_REF], Balasubramanian 1997, MacKay 2002[START_REF] Consonni | [END_REF]), tools like the Bayes factor [Jeffreys 1939] naturally include a penalization factor addressing model complexity, penalization mimicked by approximations like the Bayes Information (BIC) and the Deviance Information (DIC) criteria [Schwarz 1978, Csiszár 2000, Spiegelhalter 2002, Plummer 2008]. Posterior predictive tools have been successfully advocated in [Gelman 2013a], even though they can be criticized for multiple uses of the (same) data.

Let us recall very briefly ( referring to [Berger 1985, Robert 2001]) that the standard Bayesian approach to testing is to consider two families of models, one for each of the hypotheses under comparison,

M 1 : x ∼ f 1 (x|θ 1 ) , θ 1 ∈ Θ 1 and M 2 : x ∼ f 2 (x|θ 2 ) , θ 2 ∈ Θ 2 ,
and to associate with each of those models a prior distribution,

θ 1 ∼ π 1 (θ 1 ) and θ 2 ∼ π 2 (θ 2 ) ,
in order to compare the marginal likelihoods

m 1 (x) = � Θ 1 f 1 (x|θ 1 ) π 1 (θ 1 ) dθ 1 and m 2 (x) = � Θ 2 f 2 (x|θ 2 ) π 1 (θ 2 ) dθ 2
either through the Bayes factor or through the posterior probability, respectively:

B 12 = m 1 (x) m 2 (x) , P(M 1 |x) = ω 1 m 1 (x) ω 1 m 1 (x) + ω 2 m 2 (x)
;

the latter depends on the prior weights ω i of both models. Both testing and model selection are thus expressed as a comparison of models. The Bayesian decision step proceeds by comparing the Bayes factor B 12 to the threshold value of one or comparing the posterior probability P(M 1 |x) to a bound derived from a 0-1 loss function (or a "golden" bound like α = 0.05 inspired from frequentist practice [START_REF] Berger | [END_REF], Berger 1997, Berger 1999, Berger 2003a[START_REF] Ziliak | [END_REF]]. As a general rule, when comparing more than two models, the model with the highest posterior probability is the one selected, but this rule is highly dependent on the prior modeling, even with large datasets, which makes it hard to promote as the default solution in practical studies. Some well-documented difficulties with this traditional handling of Bayesian tests and Bayesian model choices via posterior probabilities are, among others [START_REF] Vehtari | [END_REF], Vehtari 2012]:

� a tension between using posterior probabilities as justified by a binary loss function but depending on unnatural prior weights and using Bayes factors [Jeffreys 1939] that eliminate this dependence but escape as well the direct connection with the posterior distribution, unless the prior weights are integrated within the loss function [Berger 1985, Robert 2001];;

� a subsequent and delicate interpretation (or calibration) of the strength of the Bayes factor [Jeffreys 1939, Dickey 1978[START_REF][END_REF][START_REF] Lavine | [END_REF]] towards supporting a given hypothesis or model, mostly due to the fact that it is not a Bayesian decision rule (once more, unless the loss function is artificially modified to incorporate the prior weights);

� a similar difficulty with posterior probabilities, with the correlated tendency to interpret them as p-values (rather than the opposite) when they only report through a marginal likelihood ratio the respective strengths of fitting the data to both models (and nothing about the "truth" of either model);

� a long-lasting impact of the prior modeling, meaning the choice of the prior distributions on the parameter spaces of both models under comparison, despite the existence of an overall consistency proof for the Bayes factor [Berger 2003b, Rousseau 2007, McVinish 2009];

� a discontinuity in the use of improper priors since they are not justified in most testing situations [DeGroot 1970, DeGroot 1973, Robert 2001, Robert 2014], leading to many alternative if ad hoc solutions, where the data is either used twice [Aitkin 1991, Aitkin 2010, Gelman 2013b]; or split in artifial ways [O'Hagan 1995, Berger 1996, Berger 1998, Berger 2001];

� a binary (accept vs.reject) outcome more suited for immediate decision (if any) than for model evaluation, in connection with the use of a rudimentary binary loss function that many deem unnatural [Gelman 2013a];

� a related impossibility to ascertain simultaneous misfit (i.e., a lack of fit for both models under comparison) or to detect the presence of outliers;

� a lack of assessment of the uncertainty associated with the decision itself;

� a difficult computation of marginal likelihoods in most settings [Chen 2000[START_REF][END_REF] with further controversies about which solution to adopt [START_REF][END_REF], Neal 1994, Green 1995, Chib 1995, Neal 1999, Skilling 2006, Steele 2006, Chopin 2010];

� a strong dependence of the values of posterior probabilities on conditioning statistics, which in turn undermines their validity for model assessment, as exhibited in Approximate Bayesian computation (ABC) settings by [Robert 2011] and [START_REF] Editeurs | [END_REF]];

� a temptation to create pseudo-frequentist equivalents such as q-values [START_REF][END_REF], Johnson 2013b, Johnson 2013a] with even less Bayesian justifications.

Rather than vainly attempting to solve those numerous issues in the light of the many attempts listed above, which clearly failed to produce a consensus, we therefore propose a paradigm shift in the Bayesian processing of hypothesis testing and of model selection, namely to adopt a completely novel perspective on this issue, perspective that provides a convergent and naturally interpretable solution, while allowing for a more extended use of improper priors. This approach relies on the simple representation of (or embedding into) the problem as a two-component mixture estimation problem where the weights are formally equal to 0 or 1. The mixture model [Frühwirth-Schnatter 2006] thus contains both models under comparison as extreme cases. This approach is inspired from the consistency result of [START_REF] Rousseau | [END_REF]] on estimated overfitting mixtures, where the authors established that over-parameterised mixtures can be consistently estimated, despite the parameter standing on a (or several) boundary(ies) of the parameter space. While this mixture representation is not directly equivalent to the use of a posterior probability, i.e., the posterior estimator of the mixture weight cannot be considered as a proxy to the posterior probability value, we do not perceive this as a negative feature but rather as a new tool having the potential of a better approach to testing, with a further valuable property of not expanding the number of parameters in the model (and hence keeping in line with Occam's razor, see, e.g., [START_REF] Adams | William ockham[END_REF][START_REF] Jefferys | [END_REF][START_REF] Rasmussen | [END_REF]). Our new paradigm to Bayesian testing requires a calibration of the posterior distribution of the weight of a model, while moving from the admittedly artificial and rarely understood notion of the posterior probability of a model.

The plan of the paper is as follows: Section 4.2 provides a description of the mixture model specifically created for this setting, while Section 6.2 details the implementation issues with estimating the parameters of the mixture. Section 4.3 details at great length how the mixture approach performs in the most standard i.i.d. models. Section 4.4 demonstrates its application on a survival dataset. Section 4.5 expands [START_REF] Rousseau | [END_REF]] to provide conditions on the hyperameters of the mixture model that are sufficient to achieve convergence. Section 4.6 concludes on the generic applicability of the above principle.

Testing problems as estimating mixture models

A new paradigm for testing

Given two classes of statistical models,

M 1 : x ∼ f 1 (x|θ 1 ) , θ 1 ∈ Θ 1 and M 2 : x ∼ f 2 (x|θ 2 ) , θ 2 ∈ Θ 2 ,
which may correspond to an hypothesis to be tested and its alternative, respectively, it is always possible to embed both models within an encompassing mixture model

M α : x ∼ αf 1 (x|θ 1 ) + (1 -α)f 2 (x|θ 2 ) , 0 ≤ α ≤ 1 . (4.1)
Indeed, both models correspond to very special cases of the mixture model, one for α = 1 and the other for α = 0 (with a slight notational inconsistency in the indices). 1 When considering a sample (x 1 , . . . , x n ) from one of the two models, the mixture representation still holds at the likelihood level, namely the likelihood for each model is a special case of the weighted sum of both likelihoods. However, this is not directly appealing for estimation purposes since it corresponds to a mixture with a single observation. See however [O'Neill 2014] for a computational solution based upon this representation.

What we propose in this paper is to draw inference on the individual mixture representation (6.1), acting as if each observation was individually and independently 2 produced by the mixture model. While this apparently constitutes an approximation to the real (unknown) model, except in the cases when α = 0, 1, we see several definitive advantages to this paradigm shift:

� relying on a Bayesian estimate of the weight α rather than on the posterior probability of model M 1 does produce an equally convergent indicator of which model is "true" (see Section 4.5), while removing the need of overwhelmingly artificial prior probabilities on model indices, ω 1 and ω 2 ;

� the interpretation of this estimator of α is at least as natural as handling the posterior probability, while avoiding the caricatural zero-one loss setting [DeGroot 1970, DeGroot 1973, Berger 1985]. The quantity α and its posterior distribution provide a measure of proximity to both models for the data at 1 The choice of possible encompassing models is obviously unlimited: for instance, a Geometric mixture

x ∼ fα(x) ∝ f1(x|θ1) α f2(x|θ2) 1-α
is a conceivable alternative. However, such alternatives are less practical to manage, starting with the issue of the intractable normalizing constant. Note also that when f1 and f2 are Gaussian densities, the Geometric mixture remains Gaussian for all values of α. Similar drawbacks can be found with harmonic mixtures.

2 An extension to the iid case will be considered in Example 4.3.6 for linear models. Dependent observations like Markov chains can be modeled by a straightforward extension of (6.1) where both terms in the mixture are conditional on the relevant past observations. hand, while being also interpretable as a propensity of the data to stand with (or to stem from) one of the two models. This representation further allows for alternative perspectives on testing and model choice, through the notions of predictive tools [Gelman 2013a], cross-validation [START_REF] Vehtari | [END_REF]], and information indices like WAIC [Vehtari 2012];

� the highly problematic computation [Chen 2000[START_REF][END_REF]] of the marginal likelihoods is bypassed, standard algorithms being available for Bayesian mixture estimation [START_REF] Richardson | [END_REF], Berkhof 2003, Frühwirth-Schnatter 2006, Lee 2009];

� the extension to a finite collection of models to be compared is straightforward, as this simply involves a larger number of components. This approach further allows to consider all models at once rather than engaging in pairwise costly comparisons and thus to eliminate the least likely models by simulation, those being not explored by the corresponding algorithm [START_REF][END_REF][START_REF] Richardson | [END_REF]];

� the (simultaneously conceptual and computational) difficulty of "label switching" [Celeux 2000, Stephens 2000[START_REF] Jasra | [END_REF]] that plagues both Bayesian estimation and Bayesian computation for most mixture models completely vanishes in this particular context, since components are no longer exchangeable.

In particular, we compute neither a Bayes factor3 nor a posterior probability related with the substitute mixture model and we hence avoid the difficulty of recovering the modes of the posterior distribution [Berkhof 2003, Lee 2009, Rodriguez 2014]. Our perspective is solely centered on estimating the parameters of a mixture model where both components are always identifiable;

� the posterior distribution of α evaluates more thoroughly the strength of the support for a given model than the single figure outcome of a Bayes factor or of a posterior probability. The variability of the posterior distribution on α allows for a more thorough assessment of the strength of the support of one model against the other;

� an additional feature missing from traditional Bayesian answers is that a mixture model also acknowledges the possibility that, for a finite dataset, both models or none could be acceptable. This possibility will be seen in some illustrations below (Section 4.3)

� while standard (proper and informative) prior modeling can be painlessly reproduced in this novel setting, non-informative (improper) priors now are manageable therein, provided both models under comparison are first reparameterised towards common-meaning and shared parameters, as for instance with location and scale parameters. In the special case when all parameters can be made common to both models 4 , the mixture model (6.1) can read as

M α : x ∼ αf 1 (x|θ) + (1 -α)f 2 (x|θ) , 0 ≤ α ≤ 1 .
For instance, if θ is a location parameter, a flat prior π(θ) ∝ 1 can be used with no foundational difficulty, in opposition to the testing case [DeGroot 1973, Berger 1998];

� continuing from the previous argument, using the same parameters or some identical parameters on both components is an essential feature of this reformulation of Bayesian testing, as it highlights the fact that the opposition between the two components of the mixture is not an issue of enjoying different parameters, but quite the opposite. As further stressed below, this or even those common parameter(s) is (are) nuisance parameters that need be integrated out (as they also are in the traditional Bayesian approach through the computation of the marginal likelihoods);

� even in the setting when the parameters of the mixture components, θ 1 and θ 2 , differ, they can be integrated out by mere Monte Carlo methods;

� the choice of the prior model probabilities is rarely discussed in a classical Bayesian approach, even though those probabilities linearly impact the posterior probabilities and can be argued to promote the alternative of using the Bayes factor instead. In the mixture estimation setting, prior modeling only involves selecting a prior on α, for instance a Beta B(a 0 , a 0 ) distribution, with a wide range of acceptable values for the hyperparameter a 0 , as demonstrated in Section 4.5. While the value of a 0 impacts the posterior distribution of α, it can be argued that (a) it nonetheless leads to an accumulation of the mass near 1 or 0, i.e. to favor the most likely or the true model over the other one, and (b) a sensitivity analysis on the impact of a 0 is straightforward to carry on;

� in most settings, this approach can furthermore be easily calibrated by a parametric boostrap experiment providing a posterior distribution of α under each of the models under comparison. The prior predictive error can therefore be directly estimated and can drive the choice of the hyperparameter a 0 , if need be.

4 While this may sound like an extremely restrictive requirement in a traditional mixture model, let us stress here that the presence of common parameters becomes quite natural within a testing setting. To wit, when comparing two different models for the same data, moments like E[X γ ] are defined in terms of the observed data and hence should be the same for both models. Reparametrising the models in terms of those common meaning moments does lead to a mixture model with some and maybe all common parameters. We thus advise the use of a common parameterisation, whenever possible.

Mixture estimation

Before studying the application of the above principle to some standard examples in Section 4.3, we point out a few specificities of mixture estimation in such a particular setting. While the likelihood is a regular mixture likelihood, the fact that the weights are a priori close to the boundaries means that the usual completion approach of [Diebolt 1994] is bound to be quite inefficient as soon as the sample size grows to moderate values. More precisely, if we consider a sample x = (x 1 , x 2 , . . . , x n ) from (6.1) (or assumed to be from (6.1)), the completion of the sample by the latent component indicators

ζ i (i = 1, . . . , n) leads to the completed likelihood L(θ, α 1 , α 2 | x, ζ) = n � i=1 α ζ i f (x i | θ ζ i ) = α n 1 (1 -α) n 2 n � i=1 f (x i | θ ζ i ) , (4.2)
where

(n 1 , n 2 ) = ( � n i=1 I ζ i =1 , � n i=1 I ζ i =2 ) under the constratin n = � 2 j=1 � n i=1 I ζ i =j .
This decomposition leads to a natural Gibbs implementation [Diebolt 1994] where the latent variables ζ i and the parameters are generated from their respective conditional distributions. For instance, under a Beta Be(a 1 , a 2 ) prior, α is generated from a Beta Be(a

1 + n 1 , a 2 + n 2 ).
However, while this Gibbs sampling scheme is valid from a theoretical point of view, it faces convergence difficulties in the current setting, especially with large samples, due to the prior concentration on the boundaries of (0, 1) for the mixture weight α. This feature is illustrated by Figure 4.1: as the sample size n grows, the Gibbs sample of the α's shows less and less switches between the vicinity of zero and the vicinity of one. The lack of label switching for regular mixture models is well-know, see, e.g., [Celeux 2000] and [Lee 2009]. It is due to the low probability of switching all component labels ζ i at once. This issue is simply exacerbated on the current case due to extreme values for α.

Therefore, an alternative to the Gibbs sampler is needed [Lee 2009] and we resort to a simple Metropolis-Hastings algorithm where the model parameters θ i are generated from the respective posteriors of both models (that is, based on the entire sample) and where the mixture weight α is generated either from the prior distribution or from a random walk proposal on (0, 1). It is indeed a quite rare occurrence for mixtures when we can use independent proposals. In the testing setting, the parameter θ i can be considered independently within each model and its posterior can be based on the whole dataset. (In cases when a common parameter is used in both components, one of the two available posteriors is chosen at random at each iteration, either uniformly or based on the current value of α.) The equivalent of Figure 4.1 for this Metropolis-Hastings implementation, Figure 4.2 exhibits a clear difference in the exploration abilities of the resulting chain.

We also point out that, due to the specific pattern of the posterior distribution on α accumulating most of its weight on the endpoints of (0, 1), the use of the posterior mean is highly inefficient and thus we advocate that the posterior median be instead used as the relevant estimator of α.

Illustrations

In this Section, we proceed through a series of experiments in highly classical statistical settings in order to assess the performances of the mixture estimation approach for separating the models under comparison. As we will see throughout those examples, this experimentation brings a decisive confirmation of the consistency results obtained in Section 4.5. The first two examples are direct applications of Theorem 1 while the third is an application of Theorem 2.

Example 4.3.1 For a model choice test between a Poisson P(λ) and a Geometric Geo(p) (defined as a number of failures, hence also starting at zero) distribution, we can model the mixture (6.1) towards using the same parameter λ in the Poisson P(λ) and in the Geometric Geo(p) distribution if we set p = 1 /1+λ. The resulting mixture, to be estimated, is then defined as

M α : αP(λ) + (1 -α)Geo( 1 /1+λ)
This common parameterisation allows for the call to Jeffreys' (1939) improper prior π(λ) = 1 /λ since the resulting posterior is then proper. Indeed, in a Gibbs sampling implementation, the full posterior distribution on λ, conditional on the allocation vector ζ is given by

π(λ | x, ζ) ∝ exp(-n 1 (ζ)λ + log{λ} (nx n -1) (λ + 1) -{n 2 (ζ)+s 2 (ζ)} , (4.3)
where

n 1 (ζ) = n -n 2 (ζ)
is the number of observations allocated to the Poisson component, while s 2 (ζ) is the sum of the observations that are allocated to the Geometric component. This conditional posterior is well-defined for every ζ when n > 0, which implies that the marginal posterior is similarly well-defined since ζ takes its values in a finite set. The distribution (4.3) can easily be simulated via a independent Metropolis-within-Gibbs step where the proposal distribution on λ is the Gamma distribution corresponding to the Poisson posterior. (The motivation for this choice is that, since both distributions share the same mean parameter, using the posterior distribution associated with either one of the components and all the observations should be realistic enough to produce high acceptance rates, even when the data is Geometric rather than Poisson. This is what happens in practice with acceptance rates higher than 75% in the Geometric case and close to 1 in the Poisson case. This strategy of relying on a model-based posterior as a proposal will be used throughout the examples. It obviously would not work in a regular mixture model.) Under a Be(a 0 , a 0 ) prior on α, the full conditional posterior density on α is a Be(n 1 (ζ) + a 0 , n 2 (ζ) + a 0 ) distribution and the exact Bayes factor opposing the Poisson to the Geometric models is given by

B 12 = n nxn n � i=1 x i ! Γ � n + 2 + n � i=1 x i � � Γ(n + 2).
This Bayes factor is however undefined from a purely mathematical viewpoint, since it is associated with an improper prior on the parameter [Jeffreys 1939, DeGroot 1973, Berger 1998, Robert 2009b]. The posterior probability of the Poisson model is then derived as

P(M 1 |x) = B 12 1 + B 12
when adopting (without much of a justification) identical prior weights on both models.

A first experiment in assessing our approach is based on 100 datasets simulated from a Poisson P(4) distribution. As shown in Figure 4.3, not only is the parameter λ properly estimated, but the estimation of α is very close to 1 for a sample size equal to n = 1000. In this case, the smaller the value of a 0 , the better in terms of proximity to 1 of the posterior distribution on α. Note that the choice of a 0 does not significantly impact the posterior distribution of λ. The sensitivity of the posterior distribution of α on the hyperparameter a 0 is clearly expressed by that graph. While all posterior means and medians converge to 1 in this simulation, the impact of small values of a 0 on the estimates is such that we consider values a 0 ≤ .1 as having too strong and too lengthy an influence on the posterior distribution to be acceptable.

We can also compare the outcome of a traditional (albeit invalid, since relying on improper priors) Bayesian analysis with our estimates of α. difference reinforces our earlier recommendation that smaller values of a 0 should be avoided, as they overwhelm the information contained in the data for small sample sizes.

A symmetric experiment is to study the behavior of the posterior distribution on α for data from the alternative model, i.e., a Geometric distribution. Based on 100 datasets from a Geometric G(0.1) distribution, Figure 4.7 displays the very quick convergence of the posterior median to 0 for all values of a 0 considered, even though the impact of this hyperprior is noticeable.

Example 4.3.2 For the model comparison of a normal N (θ 1 , 1) with a normal N (θ 2 , 2) distribution, we again model the mixture so that the same location parameter θ is used in both the normal N (θ, 1) and the normal N (θ, 2) distribution. Therefore, Jeffreys' (1939) noninformative prior π(θ) = 1 can be used, in contrast with the corresponding Bayes factor. Indeed, when considering the mixture of normal models, αN (θ, 1) + (1α)N (θ, 2), and a Beta B(a 0 , a 0 ) prior on α, considering the posterior distribution on (α, θ), conditional on the allocation vector ζ, leads to conditional independence between θ and α:

θ|x, ζ ∼ N � n 1 x1 + .5n 2 x2 n 1 + .5n 2 , 1 n 1 + .5n 2 � , α|ζ ∼ Be(a 0 + n 1 , a 0 + n 2 ) ,
where n i and xi denote the number of observations and the empirical mean of the observations allocated to component i, respectively (with the convention that n i xi = 0 when n i = 0. Since this conditional posterior distribution is well-defined for every possible value of ζ and since the distribution ζ has a finite support, π(θ|x) is proper.5 For the same purpose of evaluating the convergence rates of the estimates of the mixture weights, we simulated 100 N (0, 1) datasets. same concentration effect (if a lingering impact of a o ) when n increases. We also included the posterior probability of M 1 in the comparison, derived from the Bayes factor

B 12 = 2 n-1 /2 � exp 1 /4 n � i=1 (x i -x) 2 ,
with equal prior weights, even though it formally is not well-defined since based on an improper prior. The shrinkage of the posterior expectations towards 0.5 confirms our recommendation to use the posterior median instead. The same concentration phenomenon occurs for the N (0, 2) case, as illustrated on Figure 4.10 for a single N (0, 2) dataset. In order to better understand the nature of the convergence of the posterior distribution of α towards the proper limiting value, we plotted in Figure 4.9 a zoomed version of this convergence, by comparing log(n) log (1-E[α|x]) with log(1p(M 1 |x)) as the sample size n grows. Most interestingly, the variation range is of the same magnitude for both procedures, even though the choice of the hyperparameter a 0 impacts the variability of the mixture solution. This is due to the fact that the asymptotic regime is not quite reached for those sample sizes, as 1-P(M 1 |x) ≤ e -cn for some positive c with high probability, while

E[α|x] = O(n -1 /2 ), leading to log(n) log(1 -E[α|x]) � -(log n) 2 .
Furthermore, the alternative of considering the posterior probability of having the entire sample being generated from a single component is not relevant for the comparison as this estimate is always very close to zero. This means that, while α captures the model preferred by the data, the mixture modelling itself never favours a completely homogeneous sample that would come from one and only one component. By comparison, note that if we had instead called the algorithm of [van Havre 2014], we would have obtained mostly homogeneous samples for very small values of a 0 . Their algorithm is a special type of tempering MCMC, where tempering is obtained by choosing successive values of a 0 , ranging from large to very small. Example 4.3.3 We now consider a setting where we oppose a N (0, 1) model against a N (µ, 1) model, hence testing whether or not µ = 0. This being an embedded case, we cannot use an improper prior on µ and thus settle for a µ ∼ N (0, 1) prior. As discussed above in Section 6.2, Gibbs sampling applied to this mixture posterior model shows poor performances and should be replaced with a Metropolis-Hastings algorithm.

The resulting inference on the weight of the N (µ, 1) component, α, is unsurprisingly contrasted between the case when the data is distributed as N (0, 1) and when it is not from this null distribution. In the former case, obtaining values of α close to one requires larger sample sizes than in the latter case. behavior of the posterior distribution of α when the sample comes from a normal distribution N (1, 1). For a sample of size 10 2 , the accumulation of α on (.8, 1) illustrates the strength of the support for the model N (µ, 1) which is reduced with the increase of a 0 . The impact of the small sample size on the posterior distributions of α is shown in the right side of Figure 4.11 for the case where a 0 = .1 such that for n = 5 we can not recognize which model is fitter to the data.

Example 4.3.4 Inspired from [START_REF] Editeurs | [END_REF]], we oppose the normal N (µ, 1) model to the double-exponential L(µ, √ 2) model. The scale √ 2 is intentionally chosen to make both distributions share the same variance. As in the normal case in Example 4.3.2, the location parameter µ can be shared by both models and allows for the use of the flat Jeffreys' prior. As in all previous examples, Beta distributions B(a 0 , a 0 ) are compared wrt their hyperparameter a 0 .

While, in those previous examples, we illustrated that the posterior distribution of the weight of the true model converged to 1, we now consider the setting of a dataset produce by another model than those in competition, using, e.g., N (0, .7 2 ) to simulate the data. In this specific case, both posterior means and medians of α fail to concentrate near 0 and 1 as the sample size increases, as shown on Figure 4.12. So in a majority of cases in this experiment, the outcome indicates that neither of both models is favored by the data. This example does not exactly follow the assumptions of Theorem 1 since the Laplace distribution is not differentiable everywhere, however it is almost surely differentiable and it is differentiable in quadratic mean and so we expect to see the same types of behavior as predicted by Theorem 1. In this example. the Bayes factor associated with Jeffreys' prior is defined as

B 12 = exp { - � n i=1 (x i -x) 2 /2} ( √ 2π) n-1 √ n � � ∞ -∞ exp { - � n i=1 |x i -µ| / √ 2} (2 √ 2) n dµ
where the denominator is available in closed form (see 5.5). Since the prior is improper, it is formally undefined. Using nonetheless the above expression, we can compare Bayes estimators of α with the posterior probability of the model being a N (µ, 1) distribution. Based on a Monte Carlo experiment involving 100 replicas of a N (0, .7 2 ) dataset, Figure 4.13 demonstrates how the mixture estimate mostly stay away from 0 and 1 while P(M 1 |x) varies all over between 0 and 1 for all sample sizes considered here. While this is a weakly informative indication, the right hand side of Figure 4.13 shows that, on average, the posterior estimates of α converge toward a value between .1 and .4 for all a 0 while the posterior probabilities converge to .6.

In that respect, both criteria offer a similar interpretation about the data because neither α nor P (M 1 |x) confirm that either of the models is true. In the following two examples we consider regression models. Although they do not strictly speaking follow the identically setup, the methodology can be extended to this case and so can the theory, asuming that the design is random. Hence example 4.3.5 is an application of Theorem 1 while example 4.3.6 is an application of Theorem 2.

Example 4.3.5 In this example, we apply our testing strategy to a binary dataset, using the R dataset about diabetes in Pima Indian women [R Development Core Team 2006] as a benchmark [START_REF][END_REF]]. This dataset contains a randomly selected table of 200 women tested for diabetes according to WHO criteria. The response variable y is "Yes" or "No", for presence or absence of diabetes and the explanatory variable x is restricted here to the bmi, body mass index weight in kg/(height in m) 2 . For this binary dataset, either logistic or probit regression models could be suitable. We are thus comparing both fits via our method. If y = (y 1 y 2 . . . y n ) is the vector of binary responses and X = [I n x 1 ] is the n × 2 matrix of corresponding explanatory variables, the models in competition can be defined as (i = 1, . . . , n)

M 1 : y i | x i , θ 1 ∼ B(1, p i ) where p i = exp(x i θ 1 ) 1 + exp(x i θ 1 ) M 2 : y i | x i , θ 2 ∼ B(1, q i ) where q i = Φ(x i θ 2 ) (4.4)
where x i = (1 x i1 ) is the vector of explanatory variables and where θ j , j = 1, 2, is a 2 × 1 vector made of the intercept and of the regression coefficient under either M 1 or M 2 . We once again consider the case where both models share the same parameter.

However, the model is a generalized linear model and there is no moment equation that relates θ 1 and θ 2 . We thus adopt a local reparameterisation strategy by rescaling the parameters of the probit model M 2 so that the MLE's of both models coincide. This strategy follows from [Choudhury 2007] remark on the connection between the normal cdf and a logistic function

Φ(x i θ 2 ) ≈ exp(kx i θ 2 ) 1 + exp(kx i θ 2 )
and we attempt to find the best estimate of k to bring both parameters into coherency. Given

(k 0 , k 1 ) = ( � θ 01/ � θ 02 , � θ 11/ � θ 12 ) ,
ratios of the maximum likelihood estimates of the logistic model parameters to those for the probit model, we reparameterise M 1 and M 2 defined in (4.4) as

M 1 : y i | x i , θ ∼ B(1, p i ) where p i = exp(x i θ) 1 + exp(x i θ) M 2 : y i | x i , θ ∼ B(1, q i ) where q i = Φ(x i (κ -1 θ)) , (4.5) where κ -1 θ = ( θ 0/k 0 , θ 1/k 1 ).
Once the mixture model is thus parameterised, we set our now standard Beta B(a 0 , a 0 ) on the weight of M 1 , α, and choose the default g-prior on the regression parameter (see, e.g., Chapter 4. [START_REF][END_REF]]),

θ ∼ N 2 (0, n(X T X) -1 )
In a Gibbs representation (not implemented here), the full conditional posterior distributions given the allocation vector ζ are that α ∼ B(a 0 + n 1 , a 0 + n 2 ) and that

π(θ | y, X, ζ) ∝ exp �� i I ζ i =1 y i x i θ � � i;ζ i =1 [1 + exp(x i θ)] exp � -θ T (X T X)θ � 2n � × � i;ζ i =2 Φ(x i (κ -1 θ)) y i (1 -Φ(x i (κ -1 θ))) (1-y i ) (4.6)
where n 1 and n 2 are the number of observations allocated to the logistic and probit models, respectively. This conditional representation shows that the posterior distribution is then clearly defined, which is obvious when considering that for once the chosen prior is proper.

For the Pima dataset, the maximum likelihood estimates of the GLMs are θ1 = (-4.11, 0.10) and θ2 = (-2.54, 0.065), respectively, and so k = (1.616, 1.617). We compare the outcomes of this Bayesian analysis when a 0 = .1, .2, .3, .4, .5 in Table 4.1. As clearly shown by the Table, the estimates of α are close to 0.5, no matter what the value of a 0 while the estimates of θ 0 and θ 1 are very stable (and quite similar to the MLEs). We note a slight increase of α towards 0.5 as a 0 increases, but do not want to over-interpret the phenomenon. This behavior leads us to conclude that (a) none or both of the models are appropriate for the Pima Indian data; (b) the sample size may be insufficiently large for allowing a discrimination between the logit and the probit models.

Logistic model parameters Probit model parameters

a 0 α θ 0 θ 1 θ 0 k 0 θ 1 k 1 .1 .
M 1 α M 2 α True model: logistic with θ 1 = (5, 1.5) probit with θ 2 = (3.5, .8) a 0 α θ 0 θ 1 θ 0 k 0 θ 1 k 1 α θ 0 θ 1 θ 0 k 0 θ 1 k 1 .1 .
Since the benchmark dataset is apparently too small to reach the asymptotic regime, we ran a second experiment with simulated logit and probit datasets and a larger sample size n = 10, 000. For the logit model, we used the regression coefficients (5, 1.5) and for the probit model the regression coefficients (3.5, .8). The estimates of the parameters of both M α 1 and M α 2 and for both datasets are produced in Table 4.2. For every a 0 , the estimates in the true model are quite close to the true values and the posterior estimates of α are either close to 1 in the logit case and to 0 in the probit case. For this large setting, there is thus consistency in the selection of the proper model. In addition, Figure 4.14 shows that when the sample size is large enough, the posterior distribution of α concentrates its mass near 1 and 0 when the data is simulated from a logit and a probit model, respectively.

Example 4.3.6 We now examine the classical issue of variable selection in a Gaussian linear regression model. Given a vector of outcomes (y 1 , y 2 , . . . , y n ) and the corresponding explanatory variables represented by the n ×

(k + 1) matrix X = [1 n X 1 . . . X k ] (including 1 n , a first column of 1's), we assume that y | X, β, σ 2 ∼ N n (Xβ, σ 2 I n ) (4.7)
where β = (β 0 , β 1 , . . . , β k ) T is a k + 1-vector of k + 1 elements with β 0 the intercept. If we consider the generic case where any covariate could be removed from the model, we are facing the comparison of 2 k+1 -1 models, corresponding to every possible subset of explanatory variables. In our framework, this means evaluating a mixture model (6.1) with γ = 2 k+1 -1 components. For j = 1, . . . , γ, M j will denote the corresponding model, υ j the number of explanatory variables used in M j , β j the vector of the υ j regression coefficients and X j the sub-matrix of X derived from the covariate variables included in M j .

The corresponding mixture model used for testing is therefore given by

M α : y ∼ γ � j=1 α j N (X j β j , σ 2 I n ) γ � j=1 α j = 1 . (4.8)
When introducing a missing variable representation, each observation y i is associated with a missing variable ζ i taking values in 1, 2, . . . , γ. The weights of the mixture (4.8) are associated with a symmetric Dirichlet prior (α 1 , . . . , α γ ) ∼ D γ (a 0 , . . . , a 0 ). Contrary to the previous examples of this section, we now consider two different settings, corresponding to the separate versus common parameterisations of the different models M j .

Case 1. If M f denotes the full regression model, including all k explanatory variables, we impose that β j is a subvector of β f for all j's. Therefore the models M j and therefore the mixture model (4.8) all are parameterised in terms of the same β f . To simplify the notation, we will denote this common parameter vector by β = (β 0 , β 1 , . . . , β k ) T . Therefore, conditional on ζ i = j, we have

y i ∼ N � X(i) • j 2 β, σ 2 � ,
where X(i) denotes the i-th row of X and j 2 is the binary (base 2) representation of the integer j, with the convention that X(i) • j 2 means a term-by-term multiplication, i.e., that this vector contains zero entries for the components of j 2 that are equal to zero:

X(i) • j 2 = (X(i) 1 j 2 [1], . . . , X(i) k j 2 [k]) .
Assuming υ j > 0 and gathering all observations such that ζ i = j under the notation y i;ζ i =j and the corresponding covariates by X i;ζ i =j , we then have

y i;ζ i =j ∼ N υ j � X i;ζ i =j • j 2 β, σ 2 I υ j � ,
with the same convention about the term-by-term multiplication. The overall model conditional on ζ = (ζ 1 , . . . , ζ n ), the conditional distribution of the dataset is therefore

     y i;ζ i =1 y i;ζ i =2 . . . y i;ζ i =γ      n×1 =      1 i;ζ i =1 1 2 [1][X 1 ] i;ζ i =1 . . . 1 2 [k][X k ] i;ζ i =1 1 i;ζ i =2 2 2 [1][X 1 ] i;ζ i =2 . . . 2 2 [k][X k ] i;ζ i =2 . . . . . . . . . 1 i;ζ i =γ γ 2 [1][X 1 ] i;ζ i =γ . . . γ 2 [k][X k ] i;ζ i =γ      n×(k+1)      β 0 β 1 . . . β k      (k+1)×1 +ε n×1
where 1 i;ζ i =j is a υ j -dimensional vector of 1's. By convention, any value of j such that υ j = 0 does not appear in the above. If we summarize the above equation as y ζ = X ζ β + ε and use a Zellner's [Zellner 1986] G-prior,

β|σ ∼ N k+1 � M k+1 , cσ 2 (X T X) -1 � , π(σ 2 ) ∝ 1/σ 2 ,
the full conditional posterior distribution on the parameters is defined as

(α 1 , . . . , α γ )|ζ ∼ D γ (υ 1 +a 0 , . . . , υ γ +a 0 ) , β|y, ζ, σ ∼ N k+1 ( β, Σ) , σ 2 |y, β ∼ IG(a, b) ,
where

β = Σ � X T XM /cσ 2 + X ζ T y ζ/σ 2 � Σ = � X T X /cσ 2 + X ζ T X ζ/σ 2 � -1 a = (n + k + 1)/2 b = (y ζ -X ζ β) T (y ζ -X ζ β) /2 + (β-M ) T (X T X)(β-M ) /2c .
The MCMC implementation of this version of the model then leads to a straightforward Gibbs sampler.

Case 2. The alternative parameterisation of the mixture (6.1) is to consider all regression coefficients as independent between models. This means that, for j = 1, . . . , γ, the regression model M j is written as y = X j β M j + ε and that the β M j 's are independent. We still assume σ is common to all components. In this representation, we allocate a Zellner's G-prior to each parameter vector,

β M j ∼ N υ j (M j , cσ 2 ({X j } T X j ) -1 )
and, conditional on the allocation vector ζ, the full conditional posterior distributions are easily derived:

(α 1 , . . . , α γ )|ζ ∼ D γ (υ 1 +a 0 , . . . , υ γ +a 0 ) , β M j |y, σ, ζ ∼ N υ j (η j , ϕ j ) , σ 2 |y, β ∼ IG(a, b) , 4.3. Illustrations 57 where η = ϕ � {X j } T X j M j/cσ 2 + X j i;ζ i =j y i;ζ i =j/σ 2 � ϕ = � {X j } T X j /cσ 2 + {X j i;ζ i =j } T X j i;ζ i =j/σ 2 � -1 a = (n + s)/2 b = 1 2 n � i=1 γ � j=1 c(y i;ζ i =j -X j i;ζ i =j β M j ) T (y i;ζ i =j -X j i;ζ i =j β M j ) + c -1 (β M j -M j ) T ({X j } T X j )(β M j -M j )
where s is the total number of the regression coefficients of all models under comparison and where the indexing conventions are the same as in Case 1.

The comparison of the performances of the mixture approach in both cases is conducted via simulated data with k = 3 covariates, meaning that (1

≤ i ≤ n) E[y i | β, X] = β 0 + β 1 x i1 + β 2 x i2 + β 3 x i3 .
This setting thus involves 15 models to be compared (since the model where the mean of the observations is zero is of no interest). The parameters used for the data simulation are (β 0 , β 1 , β 2 , β 3 ) = (2, -3, 0, 0), σ = 1, with X 1 , X 2 and X 3 simulated from N (0, 1), B(1, .5) and U(10, 11), respectively. We are seeking to identify the true regression model M 2 :

y i = 2 -3X i1 + ε i ,
by running (Gibbs) mixture estimations algorithms.

Based on a single simulated dataset, Figure 4.15 summarizes the results of those simulations by representing the convergence of the posterior medians of the true model weight in both cases as the sample size n increases. Comments that stem from these results are that � all posterior medians of the true model weight α 2 converge to 1 when the sample size increases to n = 10, 000, which means that the mixture procedure eventually supports M 2 against the other models;

� in those graphs, the impact of the prior modeling, i.e., of the value of a 0 is such that the convergence is faster when a 0 is smaller;

� even for small sample sizes, the posterior medians of α 2 are close to 1;

� the difference between both mixture parameterisations, i.e., Case 1 and Case 2, are negligible;

� for almost every sample size and prior hyperparameter, the method concludes that M 2 is likely to be more appropriate than the others. The most interesting conclusion is therefore that using completely independent parameterisation between the components of the mixture does not induce a strong degradation in the performances of the method, although the convergence to 1 is slightly slower on the right hand side of Figure 4.15. Table 4.3 produces the posterior means of α 2 under different Dirichlet hyperparameters a 0 , which shows a stronger difference only for a 0 = 0.5, which then appears as a less reliable upper bound.

In order to assess the difference with the classical Bayesian analysis of this model, we compare our posterior means of 1α 2 with the posterior probability of M 2 computed using G-prior for the regression parameters in Figure 4.16. This picture shows that the convergence of log(1 -E(α 2 |y, X)) is faster than for log(1 -P(M 2 |y, X)). It also exhibits a difference between Cases 1 and 2 for the larger sample size, with log(1 -E(α 2 |y, X)) concentrated between -6.5 and -5 in Case 1 method, about -4 in Case 2, and about -2 for log(1 -P(M 2 |y, X)). Although those pictures are As a second check on the performances of the mixture approach for linear models, for the same set of three regressors, we simulated 50 datasets with 500 observations from each of the models M 1 , . . . , M 15 and looked at the respective averages of the Bayes estimates and of the posterior probabilities. In all cases reported in Table 4.4, the posterior means and medians support much more strongly the correct model than the posterior probability, which may sometimes get close to zero.

Case study : a survival analysis

Survival and reliability models are employed in a large number of disciplines ranging from engineering to health. An important modeling decision in these problems is the choice of the survival function. Among the many parametric alternatives, common choices include the Weibull, log-Normal, logistic, log-logistic, exponential, hypo-and hyper-exponential extensions, Gompertz, Birnhaum-Saunders, Erlang, Coxian, and Pareto distributions. The Weibull distribution is also a representative of the class of models used for extreme value modelling. Other models in this class include the extreme value, Stable, Gumbel and Frechet.

We apply here our testing paradigm to choosing between three potential survival models. Given data (x 1 , ..., x n ) with corresponding censoring indicators (c 1 , .., c n ), we wish to test the hypothesis that the data are drawn from a log-Normal(φ, κ 2 ), a Weibull(α, λ), or a log-Logistic(γ, δ) distribution. The corresponding mixture is thus given by the density

α 1 exp{-(log x -φ) 2 /2κ 2 }/ √ 2πxκ + α 2 α λ exp{-(x/λ) α }((x/λ) α-1 + α 3 (δ/γ)(x/γ) δ-1 /(1 + (x/γ) δ ) 2
where α 3 = 1α 1α 2 . A more amenable version can be obtained by working on the scale Y =log(X), which then provides a comparison between the N(θ, σ 2 ), Gumbel(µ, β), and Logistic(ξ, ζ) distributions. This gives rise to the mixture density

f θ,α (y) = α 1 exp{-(y -φ) 2 /2σ 2 } � ( √ 2πσ)+ α 2/β exp{-(y -µ)/β)} exp � -e -(y-µ)/β � + α 3 exp{-(y -ξ)/ζ} � {ζ(1 + exp{-(y -ξ)/ζ})) 2 } .
If we opt for a common parameterisation of those different models, we have the 4.4. Case study : a survival analysis 61 following moments matching equations

φ = µ + γβ = ξ σ 2 = π 2 β 2 /6 = ζ 2 π 2 /3
where γ ≈ 0.5772 is EulerâĂŞMascheroni constant. As above, this choice allows the use of a noninformative prior on the common location scale parameter, π(φ, σ 2 ) = 1/σ 2 . Once more, we use a Dirichlet prior D(a 0 , a 0 , a 0 ) on (α 1 , α 2 , α 3 ). Appendix 2 establishes that the corresponding posterior is proper provided the observations are not all equal.

A common feature in survival data is the presence of censoring. In this case, the mixture equation becomes

f θ,α (y, c) = α 1 � e -(y-φ) 2 /2σ 2 / √ 2πσ � c Φ[(y -φ)/σ] 1-c + α 2 � 1 /βe -(y-µ)/β exp � -e -(y-µ)/β �� c � exp � -e -(y-µ)/β �� 1-c + α 3 � e -(y-ξ)/ζ � � ζ(1 + e -(y-ξ)/ζ ) 2 �� c � 1/ � (1 + e -(y-ξ)/ζ ) �� 1-c
Three experiments were performed. First, the performance of the model selection approach in distinguishing between the Weibull, lognormal and log-logistic distributions was assessed by simulating 1000 observations from a Normal(0, 1) density (with no censoring), and testing a Normal versus Gumbel and Logistic distributions as described above. The experiment was then repeated using 1000 simulations from a Gumbel and then from a Logistic distribution. For illustration, the momentmatched Normal, Gumbel and Logistic densities are depicted Figure 4.17 by solid, dashed and dotted lines, respectively.

The Gibbs sampler was run for 10,000 iterations using a prior value of a 0 = 1.0 for the hyperparameter on the mixture weights. The resultant probabilities of selecting the various distributions are shown in Figure 4.17, left panel. It can be seen that in all cases, the correct model was overwhelmingly identified. As expected, the probabilities of a correct selection increase with the sample size; in an analogous experiment with n = 10 5 , all probabilities were larger than 0.90 (figures not shown).

A second experiment was undertaken to assess the influence of the hyperparameter, a 0 . The above Monte Carlo experiment was repeated with n = 1000 simulated observations, comparing the impact of four values of a 0 , namely a 0 = 0.01, 0.1, 1.0, 10.0, for all three distributions and for each pair of distributions. As illustrated on Figure 4.18 and 4.19, the probabilities of selecting a (true) Normal or (true) Gumbel model, and in agreement with earlier comparisons, the value of a 0 impacts the probability of a correct model selection, although in all cases covered by this Figure, the correct model was overwhelmingly identified (as the most likely one). Note further that in this experiment the values of a 0 were higher than those we recommender above, namely a + 0 ≤ 0.5. As before, increasing the sample size from n = 1, 000 to n = 10, 000 pushes the posterior probabilities toward the boundaries. In addition to this Monte Carlo evaluation of the mixture approach, we considered a real case study involving modeling survival times for breast cancer in Queensland, Australia. A sample of 25125 individuals with breast cancer was provided by Cancer Council Queensland. Among the subjects, 83.5% were recorded as censored and the remainder (n = 4155) were recorded as deaths from any cause. The median survival times were 4.35 and 2.02 years for each of these groups, respectively.

The response variable used in the following analyses is the hazard function, defined as the probability of death at t + Δ years given survival to t years, adjusted for age, sex and the expected mortality rate, that is, the age-and sex-adjusted background population risk of death. Of interest is whether or not this distribution is best fitted by a log Normal, Weibull, or log-Logistic distribution, or, equivalently, whether or not the log hazard is best fitted by a Normal, Gumbel, or Logistic distribution. Alternatively, it may be preferable to follow a model averaged approach, which is intrinsically part of the mixture model approach since the MCMC outcome provides in addition a posterior approximation of the overall mixture. (It could actually be argued that this approach is even better than standard model averaging as each observation in the sample selects the best fitted component of the mixture.) The choice of an appropriate model or of a combination of models is important for the prediction of survival for cancer patients, which then impacts on decisions about personalized management and treatment options.

Figure 4.20 provides histograms of both hazard and log hazard for all data and for deaths only (i.e., excluding censored observations). The corresponding q-q plots associated with fitting the three distributions, ignoring censoring, are also shown. The result of the modeling is that these distributions have different fit characteristics: whereas the Normal (and hence the log Normal) distribution fits the centre of the distribution more closely, the Weibull (and hence the Gumbel) distribution captures the tail behavior more accurately. The logistic distribution appears to have a similar fit to the Normal, but it accommodates slightly more diffuse tails. Based on a choice of hyperparameter a 0 = 1.0, the mixture test for the breast cancer data resulted in the choice of the logistic distribution with probability 0.996 (s.d. 1.4 10 -3 ), with the remaining probability mass almost equally split between the Normal and Gumbel distributions.

Asymptotic consistency

In this section we prove posterior consistency for our mixture testing procedure. More precisely we study the asymptotic behavior of the posterior distribution of α. We consider two different cases. In the first case, the two models, M 1 and M 2 , are well separated while, in the second case, model M 1 is a submodel of M 2 . We denote (right) for the non-censored data; (bottom panel) q-q plots for fitting the three distributions by π the prior distribution on (α, θ 1 , θ 2 ) and assume that θ j ∈ Θ j ⊂ R d j . We first prove that, under weak regularity conditions on each model, we can obtain posterior concentration rates for the marginal density

f θ,α (•) = αf 1,θ 1 (•) + (1 -α)f 2,θ 2 (•). Let x n = (x 1 , • • • , x n ) a n sample with true density f * . Proposition 1 Assume that, for all C 1 > 0, there exist Θ n a subset of Θ 1 × Θ 2 and B > 0 such that π [Θ c n ] ≤ n -C 1 , Θ n ⊂ {�θ 1 � + �θ 2 � ≤ n B } (4.9)
and that there exist H ≥ 0 and L, δ > 0 such that, for j = 1, 2, sup

θ,θ � ∈Θn �f j,θ j -f j,θ � j � 1 ≤ Ln H �θ j -θ � j �, θ = (θ 1 , θ 2 ), θ � = (θ � 1 , θ � 2 ) , ∀�θ j -θ * j � ≤ δ; KL(f j,θ j , f j,θ * j ) � �θ j -θ * j � .
(4.10)

We then have that, when

f * = f θ * ,α * , with α * ∈ [0, 1], there exists M > 0 such that π � (α, θ); �f θ,α -f * � 1 > M � log n/n|x n � = o p (1)
.

The proof of Proposition 1 is a direct consequence of Theorem 2.1 of [START_REF] Ghosal | [END_REF]] and is omitted for the sake of conciseness. Condition (4.10) is a weak regularity condition on each of the candidate models. Combined with condition (4.9) it allows to consider noncompact parameter sets in the usual way, see for instance [START_REF] Ghosal | [END_REF]].

It is satisfied in all examples considered in Section 4.3. We build on Proposition 1 to describe the asymptotic behavior of the posterior distribution on the parameters.

The case of separated models

Assume that both models are separated in the sense that there is identifiability:

∀α, α � ∈ [0, 1], ∀θ j , θ � j , j = 1, 2 P θ,α = P θ � ,α � ⇒ α = α � , θ = θ � , (4.11)
where P θ,α denotes the distribution associated with f θ,α . We assume that (4.11) also holds on the boundary of Θ 1 × Θ 2 . In other words, the following

inf θ 1 ∈Θ 1 inf θ 2 ∈Θ 2 �f 1,θ 1 -f 2,θ 2 � 1 > 0
holds. We also assume that, for all θ * j ∈ Θ j , j = 1, 2, if P θ j converges in the weak topology to P θ * j , then θ j converges in the Euclidean topology to θ * j . The following result then holds: Theorem 1 Assume that (4.11) is satisfied, together with (4.9) and (4.10), then for all ε > 0

π [|α -α * | > ε|x n ] = o p (1).
In addition, assume that the mapping θ j → f j,θ j is twice continuously differentiable in a neighborhood of θ * j , j = 1, 2, and that

f 1,θ * 1 -f 2,θ * 2 , ∇f 1,θ * 1 , ∇f 2,θ * 2
are linearly independent as functions of y and that there exists δ > 0 such that

∇f 1,θ * 1 , ∇f 2,θ * 2 , sup |θ 1 -θ * 1 |<δ |D 2 f 1,θ 1 |, sup |θ 2 -θ * 2 |<δ |D 2 f 2,θ 2 | ∈ L 1 . Then π � |α -α * | > M � log n/n|x n � = o p (1).
(4.12)

Theorem 1 allows for the interpretation of the quantity α under the posterior distribution. In particular, if the data x n is generated from model M 1 (resp. M 2 ), then the posterior distribution on α concentrates around α = 1 (resp. around α = 0), which establishes the consistency of our mixture approach.

Proof: Using Proposition 1, we have that

π (A n |x n ) = 1 + o p (1)
with

A n = {(α, θ); �f θ,α -f θ * ,α * � 1 ≤ δ n } and δ n = M � log n/n.
Consider a subsequence α n , P 1,θ 1n , P 2,θ 2n which converges to α, µ 1 , µ 2 where convergence holds in the sense that α n → α and P j,θ jn converges weakly to µ j . Note that µ j (X ) ≤ 1 by precompacity of the unit ball under the weak topology. At the limit

αµ 1 + (1 -α)µ 2 = α * P 1,θ * 1 + (1 -α * )P 2,θ * 2
The above equality implies that µ 1 and µ 2 are probabilities. Using (4.11), we obtain that α = α * , µ j = P j,θ * j , which implies posterior consistency for α. The proof of (4.12) follows the same line as in [START_REF] Rousseau | [END_REF]]. Consider first the case where α * ∈ (0, 1). Then the posterior distribution on θ concentrates around θ * .

Writing

L � = (f 1,θ * 1 -f 2,θ * 2 , α * ∇f 1,θ * 1 , (1 -α * )∇f 2,θ * 2 ) := (L α , L 1 , L 2 ) L " = diag(0, α * D 2 f 1,θ * 1 , (1-α * )D 2 f 2,θ * 2 ) and η = (α-α * , θ 1 -θ * 1 , θ 2 -θ * 2 ), ω = η/|η| , we then have �f θ,α -f θ * ,α * � 1 = |η| � � �ω T L � + |η| /2ω T L " ω + |η|ω 1 � ω 2 L � 2 + ω 3 L � 3 � + o(|η|) � � � (4.13) For all (α, θ) ∈ A n , set η = (α -α * , θ 1 -θ * 1 , θ 2 -θ *
2 ) goes to 0 and for n large enough there exists ε > 0 such that |αα * | + |θθ * | ≤ ε. We now prove that there exists c > 0 such that for all (α, θ)

∈ A n v(ω) = � � � � ω T L � + |η| 2 ω T L " ω + |η|ω 1 � ω T 2 L � 2 + ω T 3 L � 3 � + o(|η|) � � � � > c,
where ω is defined with respect to α, θ. Were it not the case, there would exist a sequence

(α n , θ n ) ∈ A n such that the associated v(ω n ) ≤ c n with c n = o(1)
. As ω n belongs to a compact set we could find a subsequence converging to a point ω. At the limit we would obtain ωT L � = 0 and by linear independence ω = 0 which is not possible. Thus for all (α, θ)

∈ A n |α -α * | + |θ -θ * | � δ n .
Assume now instead that α * = 0. Then define L � = (L α , L 2 ) and

L " = diag(0, D 2 f 2,θ * 2 ) and η = (α -α * , θ 2 -θ * 2 ), ω = η/|η|
and consider a Taylor expansion with θ 1 fixed, θ * 1 = θ 1 and |η| going to 0. This leads to

�f θ,α -f α * ,θ * � 1 = |η| � � � � ω T L � + |η| 2 ω T L " ω + |η|ω 1 ω 3 L � 3 � � � � + o(|η|) (4.14)
in place of (4.13) and the posterior concentration rate δ n is obtained in the same way. � We now consider the embedded case.

Embedded case

In this section we assume that M 1 is a submodel of M 2 , in the sense that θ 2 = (θ 1 , ψ) with ψ ∈ S ⊂ R d and that f 2,θ 2 ∈ M 1 when θ 2 = (θ 1 , ψ 0 ) for some given value ψ 0 , say ψ 0 = 0. Condition (4.11) is no longer verified for all α's: we assume however that it is verified for all α, α * ∈ (0, 1] and that θ * 2 = (θ * 1 , ψ * ) satisfies ψ * � = 0. In this case, under the same conditions as in Theorem 1, we immediately obtain the posterior concentration rate � log n/n for estimating α when α * ∈ (0, 1) and ψ * � = 0. We now treat the case where ψ * = 0; in other words, f * is in model M 1 .

As in [START_REF] Rousseau | [END_REF]], we consider both possible paths to approximate f * : either α goes to 1 or ψ goes to ψ 0 = 0. In the first case, called path 1,

(α * , θ * ) = (1, θ * 1 , θ * 1 , ψ) with ψ ∈ S, in the second, called path 2, (α * , θ * ) = (α, θ * 1 , θ * 1 , 0) with α ∈ [0, 1].
In either case, we write P * the distribution. We also denote F * g = � f * (x)g(x)dµ(x) for any integrable function g. For sparsity reasons, we consider the following structure for the prior on (α, θ):

π(α, θ) = π α (α)π 1 (θ 1 )π ψ (ψ), θ 2 = (θ 1 , ψ).
This means that the parameter θ 1 is common to both models, i.e., that θ 2 shares the parameter θ 1 with f 1,θ 1 .

Condition (4.11) is replaced by

P θ,α = P * ⇒ α = 1, θ 1 = θ * 1 , θ 2 = (θ * 1 , ψ) or α ≤ 1, θ 1 = θ * 1 , θ 2 = (θ * 1 , 0) (4.15) Let Θ * the above parameter set.
As in the case of separated models, the posterior distribution concentrates on Θ * . We now describe more precisely the asymptotic behavior of the posterior distribution, using [START_REF] Rousseau | [END_REF]]. We cannot apply directly Theorem 1 of [START_REF] Rousseau | [END_REF]], hence the following result is an adaptation of it. We require the following assumptions with

f * = f 1,θ * 1 .
For the sake of simplicity, we assume that Θ 1 and S are compact. Extension to non compact sets can be handled similarly to [START_REF] Rousseau | [END_REF]].

B1 Regularity: Assume that θ 1 → f 1,θ 1 and θ 2 → f 2,θ 2 are 3 times continuously differentiable and that

F *   f 3 1,θ * 1 f 3 1,θ * 1   < +∞, f1,θ * 1 = sup |θ 1 -θ * 1 |<δ f 1,θ 1 , f 1,θ * 1 = inf |θ 1 -θ * 1 |<δ f 1,θ 1 F *   sup |θ 1 -θ * 1 |<δ |∇f 1,θ * 1 | 3 f 3 1,θ * 1   < +∞, F * � |∇f 1,θ * 1 | 4 f 4 1,θ * 1 � < +∞, F *   sup |θ 1 -θ * 1 |<δ |D 2 f 1,θ * 1 | 2 f 2 1,θ * 1   < +∞, F * � sup |θ 1 -θ * 1 |<δ |D 3 f 1,θ * 1 | f 1,θ * 1 �

< +∞

B2 Integrability: There exists S 0 ⊂ S ∩ {|ψ| > δ 0 }, for some positive δ 0 and satisfying Leb(S 0 ) > 0, and such that for all ψ ∈ S 0 ,

F * � sup |θ 1 -θ * 1 |<δ f 2,θ 1 ,ψ f 4 1,θ * 1 � < +∞, F * � sup |θ 1 -θ * 1 |<δ f 3 2,θ 1 ,ψ f 3 1,θ1 * � < +∞, B3 Stronger identifiability : Set ∇f 2,θ * 1 ,ψ * (x) = � ∇ θ 1 f 2,θ * 1 ,ψ * (x) T , ∇ ψ f 2,θ * 1 ,ψ * (x) T � T .
Then for all ψ ∈ S with

ψ � = 0, if η 0 ∈ R, η 1 ∈ R d 1 η 0 (f 1,θ * 1 -f 2,θ * 1 ,ψ )+η T 1 [∇ θ 1 f 1,θ * 1 -∇ θ 1 f 2,θ * 1 ,ψ (x)] = 0 ⇔ η 1 = 0, η 2 = 0 (4.16)
We can now state the main theorem:

Theorem 2 Given the model

f θ 1 ,ψ,α = αf 1,θ 1 + (1 -α)f 2,θ 1 ,ψ ,
assume that the data is made of the n sample

x n = (x 1 , • • • , x n ) issued from f 1,θ * 1
for some θ * 1 ∈ Θ 1 , that assumptions B1 -B3 are satisfied, and that there exists

M > 0 such that π � (α, θ); �f θ,α -f * � 1 > M � log n/n|x n � = o p (1).
If the prior π α on α is a Beta B(a 1 , a 2 ) distribution, with a 2 < d 2 , and if the prior π θ 1 ,ψ is absolutely continuous with positive and continuous density at (θ * 1 , 0), then for all M n going to infinity,

π � |α -α * | > M n (log n) γ / √ n|x n � = o p (1), γ = max((d 1 + a 2 )/(d 2 -a 2 ), 1)/2, ( 4
.17)

Proof: We must find a precise lower bound on

D n := � α � Θ e ln(f θ,α )-ln(f * ) dπ θ (θ)dπ α (α)
Consider the approximating set

S n (ε) = {(θ, α), α > 1 -1/ √ n, |θ 1 -θ * 1 | ≤ 1/ √ n, |ψ -ψ| ≤ ε}
with | ψ| > 2ε some fixed parameter in S. Using the same computations as in [START_REF] Rousseau | [END_REF]], it holds that for all δ > 0 there exists C δ > 0 such that

P * � D n < e -C δ π(S n (ε))/2 � < δ. (4.18)
So that with probability greater than 1δ,

D n � n -(b+d 1 )/2 . Denote B n = {(θ, α); �f θ,α -f * � 1 ≤ M � log n/n} and A n = {(θ, α) ∈ B n ; 1 -α > z n / √ n} with z n = M n (log n) γ /
√ n and M n a sequence increasing to infinity. We split B n into

B n,1 (ε) = B n ∩ {(θ, α), θ = (θ 1 , ψ); |ψ| < ε}, B n,2 (ε) = B n ∩ B n,1 (ε) c .
To prove Theorem 2 it is enough to verify that

π(A n ) = o(n -(a 2 +d 1 )/2 ).
To simplify notations we also write δ n = M � log n/n. First we prove that for all ε > 0, A n ∩ B n,2 (ε) = ∅, when n is large enough. Let ε > 0, then for any

(θ, α) ∈ A n ∩ B n,2 (ε), We thus have |ψ| � = o(1), α = 1 + o(1) and |θ 1 -θ * 1 | = o(1)
. Consider a Taylor expansion of f θ,α around α = 1 and θ 1 = θ * 1 , with ψ fixed. This leads to

f θ,α -f * = (α -1)[f 1,θ * 1 -f 2,θ * 1 ,ψ ] + (θ 1 -θ * 1 )[∇ θ 1 f 1,θ * 1 -∇ θ 1 f 2,θ * 1 ,ψ (x)] + 1 2 (θ 1 -θ * 1 ) T � ᾱD 2 θ 1 f 1, θ1 + (1 -ᾱ)D 2 θ 1 f 2, θ1 ,ψ � (θ 1 -θ * 1 ) + (α -1)(θ 1 -θ * 1 ) T [∇ θ 1 f 1, θ1 -∇ θ 1 f 2, θ1 ,ψ ] = (α -1)[f 1,θ * 1 -f 2,θ * 1 ,ψ ] + (θ 1 -θ * 1 )∇ θ 1 f 1,θ * 1 + o(|α -1| + |θ 1 -θ * 1 |)
with ᾱ ∈ (0, 1) and θ1 ∈ (θ 1 , θ * 1 ) and the o( 1) is uniform over

A n ∩ B n,2 (ε). Set η = (α -1, θ 1 -θ * 1 ) and x = η/|η| if |η| > 0. Then �f θ,α -f * � 1 = |η| � x T L 1 (ψ) + o(1) � , L 1 = (f 1,θ * 1 -f 2,θ * 1 ,ψ , ∇ θ 1 f 1,θ * 1 )
We now prove that on

A n ∩ B n,2 (ε), �f θ,α -f * � 1 � |η|.
Assume that it is not the case then there exist c n > 0 going to 0 and a sequence (θ n , α n ) such that along that subsequence

|x T n L 1 (ψ n ) + o(1)| ≤ c n with x n = η n /|η n |.
Since it belongs to a compact, together with ψ n , any converging subsequence satisfies at the limit (x, ψ),

xT L 1 ( ψ) = 0 , which is not possible. Hence |α -1| � M √ log n/ √ n = o(M n (log n) γ / √ n), which
is not possible so that A n ∩ B n,2 (ε) = ∅ when n is large enough. We now bound π(A n ∩ B n,1 (ε)) for ε > 0 small enough but fixed. We consider a Taylor expansion around θ * = (θ * 1 , 0), leaving α fixed. Note that

∇ θ 1 f 2,θ * = ∇ θ 1 f 1,θ * 1 . We have f θ,α -f * = (θ 1 -θ * 1 ) T ∇ θ 1 f 2,θ * + (1 -α)ψ T ∇ ψ f 2,θ * 1 2 (θ -θ * ) T H α, θ(θ -θ * )
where H α, θ is the bloc matrix

H α, θ = � αD 2 θ 1 f 1, θ1 + (1 -α)D 2 θ 1 ,θ 1 f 2, θ (1 -α)D 2 θ 1 ,ψ f 2, θ (1 -α)D 2 ψ,θ 1 f 2, θ (1 -α)D 2 ψ,ψ f 2, θ �
Since H α, θ is bounded in L 1 (in the sense that each of its components is bounded as functions in L 1 ), uniformly in neighborhoods of θ * , we have writing η

= (θ 1 - θ * 1 , (1 -α)ψ) and x = η/|η|, that |η| = o(1) on A n ∩ B n,1 (ε) and �f θ,α -f * � 1 � |η| � x T ∇f 2,θ * + o(1) � ,
if ε is small enough. Using a similar argument to before, this leads to |η| � δ n on A n ∩ B n,1 (ε), so that

π (A n ∩ B n,1 (ε)) � δ d 1 n � 1 zn/ √ n (δ n /u) d 2 u b-1 du � δ d 1 +b n z b-d 2 n � n -(d 1 +a 2 )/2 M a 2 -d 2 n ,
which terminates the proof. �

Conclusion

Bayesian inference has been used in a very wide range over the past twenty years, mostly thanks to enhanced computing abilities, and many of those applications of the Bayesian paradigm have concentrated on the comparison of scientific theories and on testing of null hypotheses. Due to the ever increasing complexity of the statistical models handled in such applications, the natural and understandable tendency of practitioners has been to rely on the default solution of the posterior probability (or equivalently of the Bayes factor) without ever questioning its validity. It is only in rare cases that warnings were heeded [Robert 2011] about the poorly understood sensitivity of such tools to both prior modeling and posterior calibration. In this area, objective Bayes solutions remain tentative and do not meet with consensus. We thus believe Bayesian analysis has reached the time for a paradigm shift in the matter of hypothesis testing and model selection, albeit the solution does not have to be found outside the Bayesian paradigm, as for instance the frequentist priors of [Johnson 2013b, Johnson 2013a] and the integrated likelihood setting of [Aitkin 2010]. The novel paradigm we proposed here for Bayesian testing of hypotheses and Bayesian model comparison offers many incentives while answering some of the classical attacks against posterior probabilities and Bayes factors. Our alternative to the construction of traditional posterior probabilities that a given hypothesis is true or that the data originates from a specific model is therefore to rely on the encompassing mixture model. Not only do we replace the original testing problem with a better controlled estimation target that focus on the frequency of a given model within the mixture model, but we also allow for posterior variability over this frequency as opposed to the deterministic characteristics of the standard Bayesian approach. The posterior distribution on the weights of both components in the mixture offers a setting for deciding about which model is most favored by the data that is at least as intuitive as the sole number corresponding to either the posterior probability or the Bayes factor. The range of acceptance, rejection and indecision conclusions can easily be calibrated by simulation under both models, as well as by deciding on the values of the weights that are extreme enough in favor of one model. The examples provided in this paper have showed that the posterior medians of such weights are very quickly settling near the boundary values of 0 and 1, depending on which model is right. Even though we do not advocate such practice, it is even possible to derive a Bayesian p-value by looking at the posterior area under the tail of the distribution of the weight.

Besides decision making, another issue of potential concern about this new approach is the impact of the prior modeling. We demonstrated through all our examples that a partly common parameterisation is always feasible and hence allows for reference priors, at least on the common parameters. This proposal thus allows for a removal of the absolute prohibition of using improper priors in hypothesis testing [DeGroot 1973], a problem which has plagued the objective Bayes literature for decades. Concerning the prior on the weight parameter, we analyzed the sensitivity on the resulting posterior distribution of various prior Beta modelings on those weights. While the sensitivity is clearly present, it naturally vanishes as the sample size increases, in agreement with our consistency results, and remains of a moderate magnitude, which leads us to suggest the default value of a 0 = 0.5 in the Beta prior, in connection with both the earlier result of [START_REF] Rousseau | [END_REF]] and Jeffreys' prior in the simplest mixture setting.

A last point about our proposal is that it does not induce additional computational strain on the analysis. Provided algorithmic solutions exist for both models under comparison, such solutions can be recycled towards estimating the encompassing mixture model. As demonstrated through the various examples in the paper, the setting is actually easier than with a standard mixture estimation problem [Diebolt 1994, Marin 2005] because of the existence of common parameters that allow for the original MCMC samplers to be turned into proposals. Gibbs sampling completions are useful for assessing the potential outliers in a model but altogether not essential to achieve a conclusion about the overall problem.

Chapter 5

Supplementary material: Testing hypotheses as a mixture estimation model

In Chapter 4, we expressed how Bayesian model choice via posterior probabilities of models can be replaced by an estimation based on the probability weight of a model within a mixture model and the reasonable performance of this transformation illustrated by several examples. This chapter deals with some more Bayesian inferences, MCMC algorithms, the behavior of resulting Markov chains and also statistical tools in more details.

Mixture weight distribution

Bayesian inference of the mixture model weights is based on a beta distribution as a conjugate prior probability distribution with shape parameters that take values a 0 = .1, .2, .3, .4, .5 while a 0 = .5 yields Jeffreys prior. The identical shape parameters lead a symmetric density function about .5 that looks like a basin in a one-dimensional curve and tends to infinity in the boundaries of unit interval as shown in Figure 5.1.

If n 1 denotes the number of observations associated with the mixture component M 1 , the evolution of the posterior probability over n 1 based on the Jeffreys prior, which is Beta(.5 + n 1 , .5 + n 2 ), is shown on the right side of Figure 5.1. This evolution implies that the smaller the value of n 1 is the more the posterior density of M 1 tightens up near zero. This means that only in the case where a very large number of the observations is allocated to the model M 1 , the Bayesian estimate of corresponding mixture weight is very close to 1. This behavior plays a fundamental role in replacing the posterior probabilities of the models by the posterior estimations of the mixture model weights.

Poisson versus geometric

In the first example of Chapter 4, Poisson distribution is compared with Geometric distribution under the assumption of using the same parameter in both models. This allows us to consider a non-informative prior for the common parameter λ. Here, we firstly study this case with more details in section 5.2.1 and then in section 5.2.2, we proceed with comparing these distributions when they have deferent parameters in estimation model order to assess the behavior of the mixture model weights in the case of informative prior modeling.

Non-informative prior modeling

Let us to consider an i.i.d sample x 1 , . . . , x n from model M α . With the assumption of the same parameter in both components of the mixture model, the likelihood is such that

�(λ, α|x) = n � i=1 α exp(-λ)λ x i /x i ! + (1-α)λ x i /(1+λ) x i +1
and under the condition of the missing variable ζ i associated with each x i , we will have

�(λ, α|x, ζ) = α n 1 (1 -α) n 2 � i;ζ i =1 exp(-λ)λ x i /x i ! � i;ζ i =2 ( 1 /1+λ)(1 -1 /1+λ) x i .
When λ and α are independent and π(λ) = 1 /λ, the posterior distribution of λ is given by

π(λ|x, ζ) ∝ exp(-n 1 λ) λ � i I ζ i =1 x i + � i I ζ i =2 x i -1 /(1+λ) n 2 + � i I ζ i =2 x i
In order to simulate parameter λ from this non-standard posterior distribution, we apply the Metropolis-within-Gibbs algorithm. The related R code is available in [Kamary 2016b] in which α is simulated from the conditional posterior density Beta(a 0 + n 1 , a 0 + n 2 ) and λ from gamma proposal with parameter (

� n i=1 x i/2, n /2).
The proposal distribution of λ results in an unbiased estimate for λ because the expected value is

� n i=1 x i/n which converges to E( � n i=1 x i/n) = λ.
This algorithm works well in terms of accurately estimating the parameters λ and θ by comparing them with the true values and the resulting distribution of α supports the true model when the sample size is high enough. However, in the case where the sample size is small, this algorithm results in poor estimation of α due to the problem of label switching as shown in Figure 5.2. Another R code can be written by considering the acceptance ratio of the Metropolis-Hastings step based on the target distributions without associating any indicator variable ζ to the observations. In other words, we can sample from the following posterior densities

π(λ|x) = � n � i=1 α exp(-λ)λ x i /x i ! + (1-α)λ x i /(1+λ) x i +1 � 1 /λ π(α|x) ∝ � n � i=1 α exp(-λ)λ x i /x i ! + (1-α)λ x i /(1+λ) x i +1 � α a 0 -1 (1 -α) a 0 -1
which are both non-standard and require MCMC algorithm to sample. The advantages of the implementation of corresponding algorithm in R is that it is almost twice faster than the Gibbs sampler method and the label switching does not happen anymore in the output of α even for small sample sizes. The R code is available in [Kamary 2016b]. Our first check on convergence of Markov chains provided by this algorithm is to consider four samples of size 50, 1000, 400, 600 simulated from P(.64), P(10), Geo(.4) and P(.59), respectively, and to plot histories of λ[t], p[t], α[t], as shown in Figures 5.3 and 5.4. The trace plots indicate that the Markov chains have stabilized and appear constant over the graphs. Moreover, the chains have good mixing and are dense in the sense that they quickly traverse the support of the distribution. They are also able to explore both the tails and the mode areas efficiently. The autocorrelation plots show a very small degree of autocorrelation among the posterior samples, and the histograms estimate the posterior marginal distributions for the parameters. Note that we observe the same behavior for the simulated samples when a 0 = .2, .3, .4, .5.

In order to evaluate the accuracy of the parameter estimations, we simulate 50 datasets of sizes from 10 to 1000 once from Poisson distribution with parameter estimation model λ varied from .5 to 10 and another from geometric distribution with parameter p diverse from .06 to 1. By running MCMC algorithm for all 100 datasets, the mean absolute errors of the resulting estimates of λ and p based on the median of the simulated samples are summarized in Table 5.1 . Very small values of MAE lead us to conclude that the parameters of both models are accurately estimated in both cases. In addition to that, the evolution of the corresponding posterior estimations of α over the sample size shown in Figure 5.5 indicates the degree of the support of α toward the true model.

Informative prior modeling

Suppose that P(λ) is tested against Geo(p). We consider the following conjugate priors for the parameters λ and p,

λ ∼ G(β 1 , β 2 ); p ∼ Beta(δ 1 , δ 2 ).
The joint posterior distribution of λ and p is therefore given by The posterior samples from this non-standard posterior density are produced by implementing the Metropolis-within-Gibbs in R [Kamary 2016b]. The candidate values for λ and p are independently proposed by gamma and beta distributions. Once again, we simulate two datasets of sizes 50, 500 from P(.64) and Geo(.49), respectively, and for each dataset, we execute the R code 50 times by choosing different initial values for the chains in order to assess the convergence of the simulated samples. For both datasets, the marginal posterior distributions of λ and p are stable under 50 repetitions of the Metropolis-within-Gibbs algorithm with 10, 000 iterations, as shown in Figure 5.6. They are also identical and centered on the true values whatever the value of a 0 is. The estimated densities for α shown on the right side of Figure 5.6 have the same behavior as the case where both models under comparison share the same parameter. This means that the distributions are concentrated over 1 for the true model and the smaller the value of a 0 is, the more the densities tighten up over the boundaries of unit interval. In addition to this, for each value of a 0 , the generated samples of α resulted by 50 repetitions of MCMC algorithm have the same behavior and this guarantees the convergence of the chains toward stationarity in this case.

π(λ, p|x) ∝ � n � i=1 α exp(-λ)λ x i /x i ! + (1 -α)(1 -p) x i p � λ β 1 -1 exp(-β 2 λ)p δ 1 -1 (1-p) δ 2 -1

N (θ, 1) versus N (θ, 2)

When comparing two normal distributions with the same location parameter θ, both Gibbs sampler based on allocating missing variable ζ to the observations and the Metropolis-Hastings algorithm based on following posterior π(θ|x, α) ∝ �(θ, α|x)π(θ)

(5.1)

= n � i=1 � α exp(-(x i -θ) 2 /2) + (1 -α) exp(-(x i -θ) 2 /4)/ √ 2 � (5.2)
yield the same Bayesian inference for the mixture model. The label switching does not happen using Gibbs sampler in this case, as shown at the bottom of Figures 5.7 and 5.8. The implementations of both algorithms using programming language R are given in [Kamary 2016b].

In the Metropolis-Hastings algorithm, the proposals of θ are drawn from a normal distribution centered on the empirical mean of the observations and its standard deviation can be calibrated by an argument in the input of the corresponding function. By simulating 3 datasets of sizes 10, 510, 1000 once from N (θ, 1) and another estimation model Another data analysis is summarized in Tables 5.2 and 5.3. The tables report a series of posterior summaries related to 11 datasets simulated from normal distributions, N (θ, 1) and N (θ, 2) such that the posterior median and standard deviation of posterior draws for θ and α. The number of the observations N and the true value of θ are also listed. Included in the tables is also a convergence test using [Gelman 1992]'s criterion that is done by gelman.diag function in R, named gd.. This test is based on four MCMC chains produced in parallel starting from an arbitrary position for each parameter. Tables 5.2 and 5.3 display that for all datasets, the parameter θ is accurately estimated by the median of the posterior draws with a standard deviation less than .5 while the point estimate of α is always very close to 1 for the true model. The results of gd. test shows a clear stabilization around the target value 1 from 10,000 iterations which indicates that the four chains have converged on the same region, resulting in a perfect fit to the target.

N (

Standard normal distribution versus N (µ, 1)

The mixture of the standard normal distribution and N (µ, 1) is defined by M α : αN (0, 1) + (1α)N (µ, 1) and the conditional posterior distributions of µ and α are given by

π(µ|x) ∝ � n � i=1 α exp(-x 2 i/2) + (1 -α) exp(-(x i -µ) 2 /2) � exp(-µ 2 /2) π(α|x) ∝ � n � i=1 α exp(-x 2 i/2) + (1 -α) exp(-(x i -µ) 2 /2) � α a 0 -1 (1 -α) a 0 -1 .
The Metropolis-within-Gibbs algorithm will be applied to simulate from the conditional posteriors above and the corresponding R code can be seen in [Kamary 2016b]. Convergence verification of the chains produced by this algorithm is done by plotting the resulting posterior draws for some datasets simulated from both competing models as shown in Figure 5.9. Small autocorrelations indicate very low degree of correlation between the draws. The trace plots illustrate good mixing of the chains which are moving around the parameter space. The marginal posterior distribution of each parameter is also shown by the histograms in Figure 5.9. The plots are related to a 0 = .1 and we get the same results for the cases where a 0 = .2, .3, .4, .5. Another experiment is to run the algorithm several times for different datasets of different sizes simulated from each model under comparison. Table 5.4 lists the information about datasets such as the number of observations and the true value of the parameter µ in the case of simulating the data points from N (µ, 1). The table reports also the posterior mean and standard deviation of µ which are accurately estimated for all datasets.

Figure 5.10 displays that when N (µ, 1) is the model from which the dataset is simulated, the posterior estimate of α strongly supports this model for the data 5.4 lists more details about the datasets.

Normal versus double-exponential distribution

The mixture of a normal and a double-exponential distributions can be defined as αN (µ, 1) + (1α)L(µ, √ 2) and the conditional posterior distributions of µ and α are therefore given by

π(µ|x, α) ∝ � n � i=1 α exp(-(x i -µ) 2 /2) + (1 -α) exp(-|x i -µ| / √ 2) � π(α|x, µ) ∝ � n � i=1 α exp(-(x i -µ) 2 /2) + (1 -α) exp(-|x i -µ| / √ 2) � (α(1 -α)) a 0 -1 .
An implementation of the Metropolis-within-Gibbs can be used to simulate from these non-standard posteriors in which the parameter µ is simulated from a normal distribution centered on the empirical mean of the dataset while the standard deviation is calibrated by the user. The code in R is available in [Kamary 2016b]. In Chapter 4, we illustrated that the posterior estimates of the mixture model weights fail to concentrate near 0 or 1 even for high sample sizes when the analyzed dataset is produced by another model than those in competition. Here, we proceed by analyzing the behavior of α for the datasets simulated from one of the models under comparison. To do so, 21 samples of sizes from 5 to 1000 are simulated once from standard normal distribution another time from L(0, √ 2). The trace plots of µ are shown in Figure 5.11 and indicate the stabilization of the Markov chains over the true value 0 for any sample sizes. Figure 5.11 also shows that the higher the sample size is the more the concentration of the Markov chain is over the true value while the marginal posterior densities of α are always in favor of the true model even for the samples of size 5. In other words, the marginal posterior densities of α tightenup near 1 when the data comes from N (0, 1). We emphasize here that these results are preliminary for the convergence check. In order to check the convergence of the Markov chains, [Gelman 2003] suggests to compute the statistic gelman.diag for each scalar estimate of interest, and to continue running the chains until the statistics are all less than 1.1. We therefore compute this criterion based on four chains produced for µ and for α and list the values in Table 5.5. This table also displays a 97.5% upper limit of this diagnostic for both parameters α and µ. The values of [Gelman 2003]'s statistics shown in this table are all less than 1.1, that illustrate satisfactory convergence has been achieved.

When we test a normal

M 1 = N (µ, 1) against a double-exponential distribution M 2 = L(µ, √
2) under the flat prior by using the Bayes factor, we need to compute the marginal likelihood under both models. The marginal distribution under M 1 is and under M 2 , we have

π 1 (x) = � ∞ -∞ (2π) -n /2 exp (- � n i=1 (x i -µ) 2 /2) dµ = exp(- � n i=1 (x i -x) 2 /2) /(2π) -n /2 � ∞ -∞ exp(-n(µ-x) 2 /2)dµ = exp(- � n i=1 (x i -x) 2 /2) /(2π) -(n-
π 2 (x) = � ∞ -∞ (2 √ 2) -n exp (- � n i=1 |x i -µ| / √ 2) dµ = (2 √ 2) -n � x (1)
-∞ exp (-

� n i=1 |x i -µ| / √ 2) dµ + (2 √ 2) -n n-1 � i=1 � x (i+1)
x (i) exp (-

� n j=1 |x j -µ| / √ 2) dµ + (2 √ 2) -n � ∞ x (n)
exp (-

� n i=1 |x i -µ| / √ 2) dµ
where x (1) < . . . < x (n) . From µ < x (1) , we obtain |x i -µ| = x iµ for i = 1, . . . , n and we can rewrite the first integral as following

� x (1) -∞ exp (- � n i=1 |x i -µ| / √ 2) dµ = exp � - n � i=1 x i/ √ 2 � � x (1) -∞ exp (n µ / √ 2) dµ = √ 2 /n exp � - n � i=1 x i/ √ 2 + nx (1)/ √ 2 � Since µ > x (n) , |x i -µ| = µ -x i
, the third integral can be rewritten as

� ∞ x (n)
exp (-

� n i=1 |x i -µ| / √ 2) dµ = exp � n � i=1 x i/ √ 2 � � ∞ x (n) exp � - n � i=1 µ / √ 2 � dµ = √ 2 /n exp � n � i=1 x i/ √ 2 -nx (n)/
√ 2

�

For i = 1, . . . , n -1, we also have x (i) < µ < x (i+1) from which we deduce

|x (j) -µ| = � µ -x (j) for j < i + 1 x (j) -µ for j ≥ i + 1
and we will therefore have

� x (i+1)
x (i)

exp (-

� n j=1 |x i -µ| / √ 2) dµ = � x (i+1) x (i) exp   - i � j=1 µ-x (j)/ √ 2 - n � j=i+1 x (j) -µ / √ 2   dµ = exp �� i j=1 x (j) - � n j=i+1 x (j)/ √ 2 � � x (i+1) x (i) exp   - i � j=1 µ / √ 2 + n � j=i+1 µ / √ 2   dµ = exp �� i j=1 x (j) - � n j=i+1 x (j)/ √ 2 � � x (i+1) x (i) exp ( (n-2i)µ / √ 2) dµ estimation model
The last integral is equal to xn /2+1xn /2 when i = n /2 and we can therefore write

n-1 � i=1 � x (i+1) x (i) exp (- � n j=1 |x j -µ| / √ 2) dµ = n-1 � i=1;i� = n /2 √ 2 /(n-2i) exp �� i j=1 x (j) - � n j=i+1 x (j)/ √ 2 � exp ( (n-2i)x (i+1)/ √ 2) -exp ( (n-2i)x (i)/ √ 2) + exp   n /2 � j=1 x (j)/ √ 2 - n � j= n /2+1 x (j)/ √ 2   (x ( n /2+1) -x ( n /2) ) = n-1 � i=1;i� = n /2 √ 2 /n-2i exp � - � n j=i+1 x (j) - � i j=1 x (j) -(n-2i)x (i+1)/ √ 2 � -exp � - � n j=i+1 x (j) - � i-1 j=1 x (j) -(n-2i+1)x (i)/ √ 2 � + exp   n /2 � j=1 x (j)/ √ 2 - n � j= n /2+1 x (j)/ √ 2   (x ( n /2+1) -x ( n /2) )
The Bayes factor in Example 4.3.4 can be derivate from the marginal likelihood of the double-exponential L(µ, √ 2) model under a flat prior, that is:

� ∞ -∞ exp � -1 / √ 2 n � i=1 |x i -µ| � dµ = √ 2 /n exp    -1 / √ 2   n � j=1 x (j) -nx (1)      + √ 2 n-1 � i=1 i� =n/2 1 /n-2i exp    -1 / √ 2   n � j=i+1 x (j) - i � j=1 x (j) -(n -2i)x (i+1)      - √ 2 n-1 � i=1 i� =n/2 1 /n-2i exp    -1 / √ 2   n � j=i+1 x (j) - i-1 � j=1 x (j) -(n -2i + 1)x (i)      + (xn /2+1 -xn /2 ) exp    -1 / √ 2   n � j= n 2 +1 x (j) - n 2 � j=1 x (j)      + √ 2 /n exp    -1 / √ 2   nx (n) - n � j=1 x (j)      .
We computed the Bayes factor in Chapter 4 with the intention of comparing it with the results of our approach. However, the use of the improper prior avoid considering the Bayes factor as a validate criterion in this case.

Logistic versus probit regression model 91

Logistic versus probit regression model

We return to the problem of testing a logistic against a probit model for the binary outcomes. In chapter 4, we illustrated that when estimating a mixture of these regression models, the posterior estimate of α strongly supports the true model. This means that we can easily distinguish the true model when the sample size is large enough. The analyses are based on the mixture model defined as

M α : α( exp(y i x i θ) /(1+exp(x i θ))) + (1 -α)Φ(x i (κ -1 θ)) y i (1 -Φ(x i (κ -1 θ))) 1-y i
The likelihood and the conditional posterior distributions of θ and α derived from the g-prior and beta prior without considering missing variable ζ can be written as

�(θ, α|y, X) = � n � i=1 α( exp(y i x i θ) /(1+exp(x i θ))) + (1 -α)Φ(x i (κ -1 θ)) y i (1 -Φ(x i (κ -1 θ))) 1-y i � π(θ|y, X, α) ∝ n � i=1 α( exp(y i x i θ) /(1+exp(x i θ))) + (1 -α)Φ(x i (κ -1 θ)) y i (1 -Φ(x i (κ -1 θ))) 1-y i exp(-θ T (X T X)θ /2n) π(α|y, X, θ) ∝ n � i=1 α( exp(y i x i θ) /(1+exp(x i θ))) + (1 -α)Φ(x i (κ -1 θ)) y i (1 -Φ(x i (κ -1 θ))) 1-y i (α(1 -α)) a 0 -1
We can use the Metropolis-within-Gibbs algorithm to simulate from the conditional posteriors above in which the parameter θ is simulated according to a random walk multivariate normal distribution N (θ[t], τ Σ) starting from the maximum likelihood estimates. Σ is the asymptotic covariance matrix of the maximum likelihood estimates of the coefficients. For the R code, see [Kamary 2016b]. In chapter 4, two models were analyzed with an explanatory variable. Here, we test the regression models with intercept once over 3 another over 6 explanatory variables. The datasets are simulated once from logistic another from probit model. The explanatory variables are simulated from N (0, 1), U(0, 1),U (-1, 1), N (1, 4), U(2, 3) and U(-1, 1). The approximate Bayes estimates of θ are obtained by running the Metropoliswithin-Gibbs with scale τ = 1 over 10 4 iterations. The results summarized in Table 5.6 are slightly close to the true values.

Figure 5.12 gives an assessment of the convergence of the chains obtained for 10 4 data points simulated from the logistic model which is summarized in Table 5.6 in the case where a 0 = .5. The row sequences and the autocorrelation graphs illustrate the good mixing behavior of the chains. The posterior distributions of α displayed in Figure 5.13 are related to the four datasets used to estimate the regression coefficients in Table 5.6. Once again, the accumulation of the posterior draws of α is over 1 for the true model even if the number of the explanatory variables is more than 2.

Variable selection

The use of the mixture model for the variable selection in a Gaussian regression model is considered as a decision problem in that all potential models have to be considered as the mixture components against the mixture weights that ranks them in this context. If k is the number of predictor variables to explain the output y, every subset of explanatory variables can constitute a proper set of explanatory variables for the regression of y and the related model should be considered as a mixture model component. As an example, when k = 3, all possible models for y are shown in Table 5.7.

In the case where the mixture model is parametrized in terms of the same potential parameter β, the regression model is denoted by

y ζ = X ζ β + ε. The likelihood function is �(β, σ 2 , α|y ζ , X ζ , ζ) = γ � j=1 α υ j j ( √ 2πσ) -n exp � -(y ζ -X ζ β) T (y ζ -X ζ β) /2σ 2 �
where υ j = � n i=1 I ζ i =j ; j = 1, . . . , γ. The joint prior for β, σ 2 and α is the improper prior Table 5.7: Variable selection 5.7: All potential models for regression y when k = 3 and i = 1, . . . , n.

π(β, σ 2 , α) ∝ (σ 2 ) -k+1 /2-1 exp � -(β-M k+1 ) T (X T X)(β-M k+1 ) /2cσ 2 � γ � j=1 α a 0 -1 j
where α = (α 1 , . . . , α k ) is the vector of γ component weights. The conditional posterior distributions of β, σ 2 and α given latent variable ζ can be computed as follows

π(β|σ 2 , α, y ζ , X ζ , ζ) ∝ exp � -(y ζ -X ζ β) T (y ζ -X ζ β) /2σ 2 -(β-M k+1 ) T (X T X)(β-M k+1 ) /2cσ 2 � ∝ exp � -( X ζ β-y ζ β-M k+1 ) T   σ 2 I n 0 0 cσ 2 (X T X) -1   -1 ( X ζ β-y ζ β-M k+1 ) /2 � ∝ exp � -{( X ζ I k+1 )β-( y ζ M k+1 )} T   σ 2 I n 0 0 cσ 2 (X T X) -1   -1 {( X ζ I k+1 )β-( y ζ M k+1 )} /2

�

The last term of the equation above yields the following posterior density for

β π(β|σ 2 , α, y ζ , X ζ , ζ) ∝ exp � -(β-β) T (σ -2 X ζ T X ζ +(X T X))(β-β) /2 � β = {σ -2 X ζ T X ζ + (cσ 2 ) -1 X T X} -1 {σ -2 X ζ T y ζ + (cσ 2 ) -1 X T XM k+1 }
which implies the multivariate Gaussian distribution obtained in Chapter 4. By dropping the term that does not involve σ 2 from the multiplication of likelihood and the joint prior of β, σ 2 and α, we will obtain

π(σ 2 |β, α, y ζ , X ζ , ζ) ∝ (σ 2 ) -(n+k+1) /2-1 exp � -(y ζ -X ζ β) T (y ζ -X ζ β) /2-(β-M k+1 ) T (X T X)(β-M k+1 ) /2c /σ 2 �
from which we can easily deduce inverse-gamma distribution for σ 2 . For mixture weights, we have

π(α|ζ) = γ � j=1 α υ j +a 0 -1 j
which results in a Dirichlet distribution with the concentration parameter υ j + a 0 ; j = 1, . . . , γ. In Chapter 4, three stages Gibbs sampler algorithm is applied in order to obtain samples from the conditional posterior distributions above and the corresponding R code is available in [Kamary 2016b].

The example 4.3.6 illustrates the efficient performance of the Gibbs sampler in this context because the weight of the potential model from which the outputs are simulated converges to 1 with the sample size. However, in the case where the sample size is small, Gibbs sampler needs a high number of MCMC iterations to get the convergence of the chains and that extremely increases the system time. As an example, we analyze caterpillar dataset extracted from a 1973 study on pine processionary caterpillars [START_REF][END_REF]] when the response variable is the logarithmic transform of the average number of nests of caterpillars per tree. Three explanatory variables are considered for the regression model, which are supposed to be x 1 : the altitude, x 2 : the slope and x 3 : the number of the pines in the area. According to the classical analysis, the coefficient β 3 is not significant. This means that the appropriate model to y would have the form as M 5 in Table 5.7 in this case and the maximum likelihood estimate of the components, β0 , β1 , β2 is 4.94, -0.002, -0.035.

After running the Gibbs sampler algorithm by considering c equal to the sample size and M k+1 = 0 4 , the convergence of the chains is achieved when the number of the iterations is 10 5 as shown in Figure 5.15. In this case, the posterior distributions of β 0 , β 1 , β 2 , β 3 are centered on the same values obtained by maximum likelihood method. α 5 is concentrated over 1 which indicates that the posterior draws support the model M 5 for the output y. The result that is in agreement with the classical conclusion.

We obtain the same results when we implement the Metropolis-within-Gibbs algorithm to sample from the posterior distributions of the mixture parameters, except that we do not need to consider a large number of iterations and the convergence is achieved by producing 10 4 MCMC iterations. The code in R is also shown in [Kamary 2016b] in which the parameters β, σ 2 and α are simulated from multivariate normal, inverse-gamma and Dirichlet proposal distributions when the acceptance probabilities are based on the following posteriors

π(β|σ 2 , α, y, X) ∝ γ � j=1 α j exp � -(y-X j β j ) T (y-X j β j ) /2σ 2 � exp � -(β-M k+1 ) T (X T X)(β-M k+1 ) /2cσ 2 � π(σ 2 |β, α, y, X) ∝ (σ 2 ) -n+k+1 /2-1 γ � j=1 α j exp � -(y-X j β j ) T (y-X j β j ) /2σ 2 � exp � -(β-M k+1 ) T (X T X)(β-M k+1 ) /2cσ 2 � π(α|β, σ 2 , y, X) ∝ γ � j=1 α j exp � -(y-X j β j ) T (y-X j β j ) /2σ 2 � γ � j=1 α a 0 -1 j .
The second case studied in variable selection problem is related to the condition that the possible regression models are independent in the sense that each model estimation model has the regression coefficients that should be independently estimated from those of the other models. In this case, the number of the parameters rises very quickly by increasing the number of the explanatory variables. It means that using a large number of explanatory variables requires a huge number of parameters to be estimated. Consequently, the time system of the MCMC programs increases a lot. When we have three explanatory variables for the response y, 32 regression coefficients should be estimated for 15 potential models shown in Table 5.8.

The regression model is defined as y = X j β M j + ε for each model, M j and the likelihood conditional on missing variable ζ can therefore be written by Table 5.8: Variable selection 5.7: All potential models for regression y when k = 3 and i = 1, . . . , n.

�(β, σ 2 , α|y, X, ζ) = (2πσ 2 ) -n /2 γ � j=1 α υ j j exp � - � n i=1 (y i;ζ i =j -X j i;ζ i =j β M j ) T (y i;ζ i =j -X j i;ζ i =j β M j ) /2σ 2 �
where � γ j=1 υ j = n and the joint prior of β, α and σ 2 is

π(β, σ 2 , α) ∝ (σ 2 ) -s /2-1 exp   - γ � j=1 (β M j -M j ) T ({X j } T X j )(β M j -M j ) /2cσ 2   γ � j=1 α a 0 -1 j .
The conditional posterior distribution of β M j ; j = 1, . . . , γ is given by

π(β M j |σ 2 , α, y, X, ζ) ∝ exp � - � n i=1 (y i;ζ i =j -X j i;ζ i =j β M j ) T (y i;ζ i =j -X j i;ζ i =j β M j ) /2σ 2 � exp � -(β M j -M j ) T ({X j } T X j )(β M j -M j ) /2cσ 2 � ∝ exp � -(β M j -η) T ( {X j } T X j /cσ 2 + {X j i;ζ i =j } T X i;ζ i =j/σ 2 )(β M j -η) /2 � η = � {X j } T X j /cσ 2 + {X j i;ζ i =j } T X i;ζ i =j/σ 2 � -1 ( {X j } T X j M j/cσ 2 + X j i;ζ i =j y i;ζ i =j/σ 2 )
which leads us to deduce the multivariate Gaussian distribution for the regression model coefficients as defined in the variable selection section 4.3.6. For σ 2 , we can write

π(σ 2 |β, α, y, X, ζ) ∝ (σ 2 ) -(n+s) /2 exp � - � γ j=1 � n i=1 (y i;ζ i =j -X j i;ζ i =j β M j ) T (y i;ζ i =j -X j i;ζ i =j β M j ) /2σ 2 � exp   - γ � j=1 (β M j -M j ) T ({X j } T X j )(β M j -M j ) /2cσ 2  
that results in inverse-gamma distribution with the parameters a, b as pointed out in 4.3.6 and the posterior density of α is the same as the Case 1.. The R program related to the Gibbs sampler algorithm that samples from the conditional

E π (β 0 |y, X) E π (β 1 |y, X) E π (β 2 |y, X) E π (σ 2 |y, X) E π (α 5 |y, X) Case 1.
5.18 -0.003 -0.039 0.54 0.96 Case 2.

5.21 -0.003 -0.051 0.56 0.77 MLE.

4.94 -0.002 -0.035 0.65 Table 5.9: Variable selection 5.7: Point estimate of the regression coefficients of the model M 5 , σ 2 and α 5 based on 10, 000 MCMC iterations when a 0 = 0.5.

posteriors of β, σ 2 and α is indicated in [Kamary 2016b]. The analyses spoken of in the Case 2. part of the variable selection section in Chapter 4 are based on this Gibbs sampler algorithm. However, running this program is time consuming even for a dataset with small sample sizes and with 10 4 MCMC iterations. For example, for the Caterpillar dataset, the convergence of the chains is achieved when we produce a large number of MCMC iterations. However, time system is much more than the one for the Case 1.. We can also use the Metropolis-within-Gibbs algorithm in this case where the acceptance probability of the proposal distributions of α, σ 2 and β M j ; j = 1, . . . , γ are based on the following posterior distributions

π(β M j |σ 2 , α, y, X) ∝ γ � j=1 α j exp � -(y-X j β M j ) T (y-X j β M j ) /2σ 2 � exp � -(β M j -M j ) T ({X j } T X j )(β M j -M j ) /2cσ 2 � π(σ 2 |β, α, y, X) ∝ (σ 2 ) -(n+s) /2 γ � j=1 α j exp � -(y-X j β M j ) T (y-X j β M j ) /2σ 2 � exp   - γ � j=1 (β M j -M j ) T ({X j } T X j )(β M j -M j ) /2cσ 2   π(α|β, σ 2 , y, X) ∝ γ � j=1 α j exp � -(y-X j β M j ) T (y-X j β M j ) /2σ 2 � γ � j=1 α a 0 -1 j .
The related code in R is pointed out in [Kamary 2016b]. After running the algorithm for the Caterpillar dataset, a graphical convergence check for the posterior draws of the regression coefficients is shown in Figures 5.15 and 5.16 that illustrate the Markov chains have stabilized and mixed very well in this case. The Bayes estimates of the regression components of the model M 5 and the maximum likelihood estimates are displayed in Table 5.9. The Bayesian outputs of both cases are very close to the maximum likelihood estimates while the posterior estimate of α 5 is also strongly in favor of M 5 in both cases. Note that in Case 1., the point estimate of α 5 is closer to 1 than in Case 2.. 5.8 when a 0 = .5.

Propriety of the posterior in the case study of Section 4

To prove the propriety of the posterior it is enough to prove the propriety of the subposterior distribution associated to each component since the parameter (θ, σ) is shared between the components. It is known that in the case of a Gaussian model N (θ, σ) the posterior associated to the prior π(θ, σ) = 1/σ is proper as soon as n ≥ 2 and at least 2 observations are distinct. We now show that this results extends to the case of a Gumbel(θ, σ 2 ) and of a Logistic(θ, σ). Let I denote the estimation model 5.8 when a 0 = .5.

marginal likelihood, in the Gumbel case.

I = � R×R + 1 σ n+1 exp � - n � i=1 (Y i -θ)/σ � exp � - n � i=1 e -(Y i -θ)/σ � dθdσ We set γ n = � n i=1 e -Y i /σ , then I(σ) = e -n Ȳn/σ � R exp(nθ/σ) exp(-γ n e θ/σ )dθ ∝ e -n Ȳn/σ σ � R u n-1 exp(-γ n u)du ∝ e -n Ȳn/σ σγ -n n ∝ exp � - n Ȳn σ � � n � i=1 e -Y i σ � -n . So I ∝ � R + σ -n exp � - n Ȳn σ � � n � i=1 e -Y i σ � -n dσ = � R + σ -n � n � i=1 e -1 σ (Yi-Ȳn) � -n dσ < +∞ if only if min i � Y i -Ȳn � < 0.
This is almost surely true when n ≥ 2. We now study the Logistic case, using similar computations, so that

I ∝ � R×R + e -n Ȳn σ σ (n+1) e θn/σ � i (1 + e -Y i /σ e θ/σ ) 2 dθdσ ∝ � R + e -n Ȳn σ σ n � R u n-1 � i (1 + e -Y i /σ u) 2 dudσ ≤ � R + e -n Ȳn σ σ n � R u n-1 (1 + u n-1 e -(n-1) Ȳn/σ max i e -(Y i -Ȳn)/σ ) 2 du ∝ � R + 1 σ n � R u n-1 (1 + u n-1 max i e -(Y i -Ȳn)/σ ) 2 du ∝ � R + 1 σ n e 2n min i (Y i -Ȳn)/σ dσ < +∞ if and only if min i � Y i -Ȳn � < 0.
Thus means the observations cannot be all equal.

Introduction

A mixture density is traditionally represented as a weighted average of densities from standard families, i.e.,

f (x|θ, p) = k � i=1 p i f (x|θ i ) k � i=1 p i = 1 . (6.1)
Each component of the mixture is characterized by a component-wise parameter θ i and the weights p i of those components translate the importance of each of those components in the model. This particular representation gives a separate meaning to each component through its parameter θ i , even though there is a well-known lack of identifiability in justifications, in particular when over-fitting mixtures [START_REF] Rousseau | [END_REF][START_REF] Kamary | [END_REF]], but most attempting to overcome the methodological hurdles in estimating mixture model [Chib 1995, Neal 1999, Berkhof 2003, Marin 2005, Frühwirth-Schnatter 2006, Lee 2009[START_REF] Mengersen | [END_REF].

In this paper, we introduce and study the global mean-variance reparameterisation (Section 6.2), which main consequence is to constrain all other parameters to a compact space. We study several possible parameterisations of that kind and demonstrate that the improper Jeffreys-like prior associated with them is proper. In Section 6.3, we propose some MCMC implementation to estimate the parameters of the mixture, discussing label switching (Section 6.3.2) and its resolution by tempering. Extensions to non-Gaussian mixtures are briefly discussed in Section 6.6.

Mixture representation

Mean-variance reparameterisation

Let us first recall how both mean and variance of a mixture distribution can be represented in terms of the mean and variance parameters of the component of the mixture:

Lemma 1 If µ i and σ 2 i denote the mean and variance of the distribution with density f (•|θ i ), respectively, the mean of the mixture distribution (6.1) is given by

E θ,p [X] = k � i=1 p i µ i
and its variance by

var θ,p (X) = k � i=1 p i σ 2 i + k � i=1 p i (µ 2 i -E θ,p [X] 2 )
Proof: The population mean given by

E θ,p [X] = k � i=1 p i E f (•|θ i ) [X] = k � i=1 p i µ i where E f (•|θ i ) [X]
is the expected value component i. Similarly, the population variance is given by

var θ,p (X) = k � i=1 p i E f (•|θ i ) [X 2 ] -E θ,p [X] 2 = k � i=1 p i (σ 2 i + µ 2 i ) -E θ,p [X] 2 ,
which concludes the proof � Chapter 6. Non-informative reparameterisations for location-scale mixtures

For any location-scale mixture, we then propose a reparameterisation of the mixture model that starts by scaling all parameters in terms of its global mean µ and global variance σ 2 . For instance, we can switch to the representation

µ i = µ + σα i and σ i = στ i (6.2)
of the component-wise parameters, where τ i > 0 and α i ∈ R. This is formally equivalent to the reparameterisation of [Mengersen 1996], except that they put no special meaning on the global mean and variance parameters. Once the global mean and variance are set, this imposes natural constraints on the other parameters of the model. For instance, setting the global variance to σ 2 implies that (µ 1 , . . . , µ k , σ 1 , . . . , σ k ) belongs to a specific ellipse conditional on the weights and σ 2 , by virtue of Lemma 1.

Considering the α i 's and the τ i 's in (6.2) as the new parameters of the components, the following result states that the global mean and variance parameters are the sole freely varying parameters. In other words, once both the global mean and variance are set, there exists a parameterisation such that all remaining parameters of a mixture distribution are restricted to a compact set, which is most helpful in selecting a non-informative prior distribution.

Lemma 2

The parameters α i and τ i in (6.2) are constrained by

k � i=1 p i α i = 0 and k � i=1 p i τ 2 i + k � i=1 p i α 2 i = 1 .
Proof: The result is a trivial consequence of Lemma 1. The population mean is

E θ,p [X] = k � i=1 p i µ i = k � i=1 p i (µ + σα i ) = µ + k � i=1 p i α i = µ
and the first constraint follows. The population variance is

var θ,p (X) = k � i=1 p i σ 2 i + k � i=1 p i (µ 2 i -E θ,p [X] 2 ) = k � i=1 p i σ 2 τ 2 i + k � i=1 p i p i (µ 2 + 2σµα i + σ 2 α 2 i -µ 2 ) = k � i=1 p i σ 2 τ 2 i + k � i=1 p i σ 2 α 2 i = σ 2
The last equation simplifies to the second constraint above. �

Reference priors

The constraints in Lemma 2 define a set of values of (p 1 , . . . , p k , α, . . . , α, τ, . . . , τ ) that is obviously compact. From a Bayesian perspective, this allows for the call to uniform and other non-informative proper priors, conditional on (µ, σ). Furthermore, since (µ, σ) is a location-scale parameter, we may invoke [Jeffreys 1939] to use the Jeffreys prior π(µ, σ) = 1/σ on this parameter, even though this is not the genuine Jeffreys prior for the mixture model [START_REF] Grazian | [END_REF]]. In the same spirit as [Robert 1998] who established properness of the posterior distribution derived by [Mengersen 1996], we now establish that this choice of prior produces a proper posterior distribution for a minimal sample size of two.

Theorem 3 The posterior distribution associated with the prior π(µ, σ) = 1/σ and with the likelihood derived from (6.1) is proper when the components f (•|µ, σ) are Gaussian densities, provided (a) proper distributions are used on the other parameters and (b) there are at least two observations in the sample.

Proof: When n = 1, it is easy to show that the Jeffreys posterior is not proper. The marginal likelihood is then

M k (x 1 ) = k � i=1 � p i f (x 1 |µ + σα i , σ 2 τ 2 i )π(µ, σ, p, α, τ ) d(µ, σ, p, α, τ ) = k � i=1 � �� p i √ 2πσ 2 τ i exp � -(x 1 -µ -σα i ) 2 2τ 2 i σ 2 � d(µ, σ) � π(p, α, τ ) d(p, α, τ ) = k � i=1 � �� ∞ 0 p i σ dσ � π(p, α, τ ) d(p, α, τ )
The integral against σ is then not defined.

For two data-points,

x 1 , x 2 ∼ � k i=1 p i f (µ + σα i , σ 2 τ 2 i ), the associated marginal likelihood is M k (x 1 , x 2 ) = � 2 � j=1 � k � i=1 p i f (x j |µ + σα i , σ 2 τ 2 i ) � π(µ, σ, p, α, τ ) d(µ, σ, p, α, τ ) = k � i=1 k � j=1 � p i p j f (x 1 |µ + σα i , σ 2 τ 2 i )f (x 2 |µ + σα j , σ 2 τ 2 j )π(µ, σ, p, α, τ ) d(µ, σ, p, α, τ ) .
If all those k 2 integrals are proper, the Jeffrey posterior distribution is proper. An arbitrary integral

(1 ≤ i, j ≤ k) in this sum leads to � p i p j f (x 1 |µ + σα i , σ 2 τ 2 i )f (x 2 |µ + σα j , σ 2 τ 2 j )π(µ, σ, p, α, τ ) d(µ, σ, p, α, τ ) = � � � p i p j 2πσ 3 τ i τ j exp � -(x 1 -µ -σα i ) 2 2τ 2 i σ 2 + -(x 2 -µ -σα j ) 2 2τ 2 j σ 2 � d(µ, σ) � π(p, α, τ ) d(p, α, τ ) = � � � ∞ 0 p i p j √ 2πσ 2 � τ 2 i + τ 2 j exp � -1 2(τ 2 i + τ 2 j ) � 1 σ 2 (x 1 -x 2 ) 2 + 2 σ (x 1 -x 2 )(α i -α j ) +(α i -α j ) 2 �� dσ � π(p, α, τ ) d(σ, p, α, τ ) .
Chapter 6. Non-informative reparameterisations for location-scale mixtures Substituting σ = 1/z, the above is integrated with respect to z, leading to

� � � ∞ 0 p i p j √ 2π � τ 2 i + τ 2 j exp � -1 2(τ 2 i + τ 2 j ) � z 2 (x 1 -x 2 ) 2 + 2z(x 1 -x 2 )(α i -α j ) +(α i -α j ) 2 �� dz � π(p, α, τ ) d(p, α, τ ) = � � � ∞ 0 p i p j √ 2π � τ 2 i + τ 2 j exp � -(x 1 -x 2 ) 2 2(τ 2 i + τ 2 j ) � z + α i -α j x 1 -x 2 � 2 � dz � π(p, α, τ ) d(p, α, τ ) = � p i p j |x 1 -x 2 | Φ   - α i -α j x 1 -x 2 |x 1 -x 2 | � τ 2 i + τ 2 j   π(p, α, τ ) d(p, α, τ ) ,
where Φ is the cumulative distribution function of the standardized Normal distribution. Given that the prior is proper on all remaining parameters of the mixture and that the integrand is bounded by 1/|x 1 -x 2 |, it integrates against the remaining components of θ.

Let us now consider the case n ≥ 3. Since the posterior π(θ|x 1 , x 2 ) is proper, it constitutes a proper prior when considering only the observations x 3 , . . . , x n . Therefore, the posterior is almost everywhere proper. �

Further reparameterisations

Before proposing relevant priors, let us note that the constraints in Lemma 2 suggest a new reparameterisation (among many possible ones): this reparameterisations uses the weights p i in the definition of the component parameters, as to achieve a more generic constraint. The component location and scale parameters in (6.2) can indeed be reparameterised as

α i = σγ i / √ p i and τ i = ση i / √ p i , leading to the mixture representation f (x|θ, p) = k � i=1 p i f (x|µ + σγ i / √ p i , ση i / √ p i ) , η i > 0 , (6.3) Given (p 1 , • • • , p k ), these new parameters are constrained by k � i=1 √ p i γ i = 0 and k � i=1 (η 2 i + γ 2 i ) = 1 ,
which means that (γ 1 , . . . , η k ) belongs to an hypersphere of R 2k intersected with an hyperplane of this space. Given these constraints, further simplifications via new reparameterisations can be contemplated, as for instance separating mean and variance parameters in (6.3) 6.2. Mixture representation 109 by introducing a radius ϕ such that

k � i=1 γ 2 i = ϕ 2 and k � i=1 η 2 i = 1 -ϕ 2 .
This choice naturally leads to a hierarchical prior where, e.g., ϕ 2 and (p 1 , . . . , p k ) are distributed from a Be(a 1 , a 2 ) and a Dir(α 0 , . . . , α 0 ) distributions, respectively, while the vectors (γ 1 , . . . , γ k ) and (η 1 , . . . , η k ) are uniformly distributed on the spheres of radius ϕ and � 1ϕ 2 , respectively, under the additional linear constraint

� k i=1 √ p i γ i = 0.
We now describe how this reparameterisation leads to a practical construction of the constrained parameter space, for an arbitrary number of components k. 6.2.3.1 Spherical coordinate representation of the γ's.

The vector (γ 1 , . . . , γ k ) belongs both to the hypersphere of radius ϕ and to the hyperplane orthogonal to ( √ p 1 , . . . , √ p k ). Therefore, (γ 1 , . . . , γ k ) can be expressed in terms of spherical coordinates within that hyperplane. Namely, if (� 1 , . . . , � k-1 ) denotes an orthonormal basis of the hyperplane, (γ 1 , . . . , γ k ) can be written as

(γ 1 , . . . , γ k ) = ϕ cos(� 1 )� 1 +ϕ sin(� 1 ) cos(� 2 )� 2 +. . .+ϕ sin(� 1 ) • • • sin(� k-2 )� k-1
with the angles � 1 , . . . , � k-3 in [0, π] and � k-2 in [0, 2π]. The s-th orthonormal base � s can be derived from the k-dimensional orthogonal vectors � � s where

� � 1,j = � - √ p 2 , j = 1 √ p 1 , j = 2 0, j > 2
and the s-th vector is given by

� � s,j =          -(p j p s+1 ) 1/2 � � � s l=1 p l � 1/2 , s > 1, j ≤ s � � s l=1 p l � 1/2 , s > 1, j = s + 1 0, s > 1, j > s + 1
Note the special case of k = 2 since the angle � 1 is then missing. In this special case, the mixture location parameter is defined by (γ 1 , γ 2 ) = ϕ� 1 and ϕ takes both positive and negative values. In the general setting, the parameter vector

(γ 1 • • • , γ k ) is a transform of (ϕ 2 , p 1 , • • • , p k , � 1 , • • • , � k-2 ). A natural reference prior for � is made of uniforms, � 1 , • • • , � k-3 ∼ U[0, π] and � k-2 ∼ U[0, 2π]
, although other choices are obviously possible and should be explored to test the sensitivity to the prior.

Chapter 6. Non-informative reparameterisations for location-scale mixtures 6.2.3.2 Dual representation of the η i 's.

For the component variance parameters, the vector (η 1 , • • • , η k ) belongs to the kdimension sphere of radius � 1ϕ 2 . A natural prior is then a Dirichlet distribution with common hyperparameter a,

π(η 2 1 , • • • , η 2 k , ϕ 2 ) = Dir(α, • • • , α) If k is small enough, (η 1 , • • • , η k )
can then be simulated from the corresponding posterior with no computational challenge. However, as k increases, sampling may become more delicate and benefits from a similar spherical reparameterisation. In this approach, the vector

(η 1 , • • • , η k ) is rewritten through spherical coordinates with angle components (ξ 1 , • • • , ξ k-1 )
,

η i =                  � 1 -ϕ 2 cos(ξ i ) , i = 1 � 1 -ϕ 2 i-1 � j=1 sin(ξ j ) cos(ξ i ) , 1 < i < k � 1 -ϕ 2 i-1 � j=1 sin(ξ j ) , i = k Unlike �, the support for all angles ξ 1 , • • • , ξ k-1 is limited to [0, π/2],
due to the positivity requirement on the η i 's. In this case, a reference prior on the angles is

(ξ 1 , • • • , ξ k-1 ) ∼ U([0, π/2] k-1 ) ,
while again other choices are possible.

MCMC implications

The Metropolis-within-Gibbs sampler

Given the reparameterisations introduced in Section 6.2, different MCMC implementations are possible and we investigate in this section some of these. To this effect, we distinguish between two cases: (i) only (µ 1 , • • • , µ k ) is expressed in spherical coordinates; and (ii) both the µ i 's and the σ i 's are associated with spherical coordinates.

Although the target density is similar to the target explored by early Gibbs samplers in [Diebolt 1990] and [START_REF] Gelman | [END_REF]], simulating directly the new parameters implies managing constrained parameter spaces. The hierarchical nature of the parameterisation also leads us to consider a block Gibbs sampler that coincides with this hierarchy. Since the corresponding full conditional posteriors are not in closed form, a Metropolis-within-Gibbs sampler is implemented here with random walk proposals. In this approach, the scales of the proposal distributions are automatically calibrated towards optimal acceptance rates [Roberts 1997, Roberts 2001, Roberts 2009, Rosenthal 2011]. Convergence of a simulated chain is assessed based on the rudimentary convergence monitoring technique of [Gelman 1992]. The description of the algorithm is provided by the pseudo-code version in Figure 6.1. Note that the Metropolis-within-Gibbs version does not rely on latent variables and complete likelihood as in [Tanner 1987] and [Diebolt 1990]. Following the adaptive MCMC method in Section 3 of [Roberts 2009], we derive the optimal scales associated with proposal densities, based on 10 batches with size 50. The scales ε are identified by a subscript with the corresponding parameter.

For the reparameterisation (i), all steps are the same except that steps 2.5 and 2.7 are combined together and that ((ϕ 2 ) (t) , (η 2 1 ) (t) , . . . , (η 2 k ) (t) ) is updated in the same manner. One potential proposal density is a Dirichlet distribution,

((ϕ 2 ) � , (η 2 1 ) � , . . . , (η 2 k ) � ) ∼ Dir((ϕ 2 ) (t-1) ε, (η 2 1 ) (t-1) ε, . . . , (η 2 k ) (t-1) ε) .
Alternative proposal densities will be discussed later along with simulation studies in Section 4.

Removing and detecting label switching

The standard parameterisation of mixture models contains weights {p i } k i=1 and component-wise parameters {θ i } k i=1 as shown in (6.1). The likelihood function is invariant under permutations of the component indices. If an exchangeable prior is chosen on weights and component-wise parameters, the posterior density reproduces the likelihood invariance and component labels are not identifiable. This phenomenon is called label switching and is well-studied in the literature [Celeux 2000, Stephens 2000, Frühwirth-Schnatter 2001, Frühwirth-Schnatter 2004[START_REF] Jasra | [END_REF]. This means that the posterior distribution consists of k! symmetric modes and a Markov chain with such target distribution is expected to explore all of them. However, a chain often fails and rather ends up exploring a particular mode.

In our reparameterisation of Gaussian mixture models, each component mean and variance are functions of angular and radius parameters with weights. The mapping between both parameterisations is a one-to-one map conditional on the weights. In other words, there are unique component-wise means and variances given particular values for angular and radius parameters and weights. Although the new parameterisation is not exchangeable, due to the choice of the orthogonal basis, adopting an exchangeable prior on the weights (e.g., a Dirichlet distribution with a common parameter) and uniform priors on all angular parameters leads to an exchangeable posterior on the natural parameters of the mixture. Therefore, label switching should also occur with this prior modeling.

When an MCMC chain manages to jump between modes, the inference on each of the mixture components becomes harder [Geweke 2007]. To get component-specific inference and to give a meaning to each component, various relabelling methods have been proposed in the literature (see, e.g., [Frühwirth-Schnatter 2004]). A first available alternative is to reorder labels so that the mixture weights are in increasing order [Frühwirth-Schnatter 2001]. A second alternative method proposed by, e.g., [Lee 2009] is that labels are reordered towards producing the shortest distance when the adaptive Metropolis-within-Gibbs sampler is used.

The scales of the various proposals are determined by aiming at [Roberts 1997] goal of an average acceptance rate of either 0.44 or 0.234 depending on the dimension of the simulated parameter. As shown in Table 6.1, an adaptive Metropolis-within-Gibbs strategy manages to recover acceptance rates close to optimal values. �

Having exposed how our sampler behaves we now discuss a second example, in which we briefly outline how this method may behave for a benchmark dataset with a slightly larger sample size.

Example 6.4.2 We now analyze the benchmark Old Faithful dataset, available from R, using the 272 observations of eruption times and a mixture with two components. The empirical mean and variance of the observations are (3.49, 1.30).

When using Proposal 1, the optimal scales ε µ , ε p , ε after 50, 000 burn-in iterations are 0.07, 501. 1, 802.19, respectively. The posterior distributions of the generated samples shown in Figure 6.4 demonstrate a strong concentration of (µ, σ 2 ) near the empirical mean and variance. Trace plots for the other parameters indicate a high dependence between successive iterations.There is a strong indication that the chain gets trapped into a single mode of the posterior density. In Section 6.5, we reanalyse this dataset when using parallel tempering. �

The general case

We now consider the general case of estimating a reparameterised mixture for any k when the variance vector (η 2 1 , . . . , η 2 k ) also has the spherical coordinate system as represented in Section 6.2.3. Example 6.4.3 We simulated 50 data points from the mixture 0.27N (-4.5, 1) + 0.4N (10, 1) + 0.33N (3, 1) .

Running our adaptive Metropolis-within-Gibbs algorithm shows that the simulated samples are quite close to the true values. However, the sampler has apparently visited only one of the posterior modes. This lack of label switching helps us in producing point estimates directly from this MCMC output [Geweke 2007] but this also shows an incomplete convergence of the MCMC sampler [Celeux 2000]. When considering the new parameters of this mixture, the single � plays a significant role in the lack of label switching since transforming � to π -� swaps first and second components. If we restrict the proposal on � to step 2.4 of the Metropolis-within-Gibbs algorithm, namely using only a uniform U(0, 2π) distribution, Figure 6.5 shows that the MCMC chains of the p i 's are both well-mixed and exhibiting strong exchangeability. However, the corresponding acceptance rate is quite low at 0.051. However, Bayesian inference for parameters related to individual components of the mixture using averaging over posterior draws is not possible in this case since the posterior means of the component specific parameters such as p, µ i , σ i ; i = 1, 2, 3 are the same for all components. We therefore revert to both methods of k-means clustering algorithm presented at the beginning of this section and removing label switching based on the distance between posterior sample and MAP estimate which are shown in left and right sides of Table 6.2, respectively. Bayesian estimations computed by both methods are almost identical and all parameters of the mixture distributions are accurately estimated in comparison with those of the true model with the acceptance rates of the proposal distributions of the Metropolis-within-Gibbs very close to the optimal ones. When computing point estimates of the natural parameters of the components, we obtain the maximum errors of 0.08 and 0.11 for µ and σ, respectively. The average absolute error over the 20 samples is quite low. Furthermore, when comparing the true and estimated mixtures, we can resort to the Kullback-Leibler divergence. For the 20 simulated samples, the maximum value is 0.02, which means an information loss of at most 2%. If we consider the upper bound introduced by [Sayyareh 2011] on Kullback-Leibler divergence, the obtained values indicates a good similarity between P true and P estimated and illustrates the consistency of the estimates resulting from our Metropolis-within-Gibbs algorithm. � Example 6.4.5 When an MCMC chain converges to a very small value for at least one component weight p i , this may lead to an extremely large mean or large variance in the corresponding component. This happens partly because there is hardly any information from the data for this component and partly because the new parameters are functions of 1/ √ p i . We may thus face extreme points in the simplex parameter spaces. This phenomenon is illustrated with the Galaxy dataset, a constant benchmark for mixture estimation [Roeder 1990[START_REF] Richardson | [END_REF] It is seen that the proposed prior is more dispersed for µ 1 and p 1 and is very skewed toward 0 for σ 1 with long tail. When k = 6, bayesm yields a more concentrated prior for p to accommodate all components and the proposed prior becomes dispersed to give flexible support on component-wise location and scale.

Parallel tempering

In Example 6.4.2 we have seen that for the Old Faithful dataset, the multimodality of the mixture model is not reproduced in the MCMC output, which means the adaptive Metropolis-within-Gibbs sampler cannot escape one of the modes. In this case, parallel tempering may be used [Marinari 1992, Neal 1996]. This method allows for better mixing in multimodal target distributions, when using straightforward Metropolis-Hastings algorithms fail [START_REF] Miasojedow | [END_REF]]. It is indeed designed For the prior by bayesm (dashed lines), hyperparameters are α 0 = 5 for k = 4 and α 0 = 25 for k = 6 while μ = 0 and ν = 3.

to overcome low probability regions between modal areas. Given the posterior density f (θ|x), we define tempered versions f β (θ|x) ∝ f (θ|x) β , where 0 ≤ β ≤ 1 is the inverse temperature and β = 1 corresponds to the original target distribution [Geyer 1991]. The tempered MCMC algorithm then runs a basic MCMC algorithm on a range of tempered distribution and, at each iteration, the current samples are considered for potential exchanges between adjacent temperatures, with a Metropolis-Hastings acceptance probability

α h = min � 1, f β h-1 (θ (t) h )f β h (θ (t) h-1 ) f β h-1 (θ (t) h-1 )f β h (θ (t) h ) � ,
as the chances of accepting a swap are higher for nearby temperatures. Proposal scales are calibrated by adaptive MCMC method and is used for all tempered versions of the target. Temperatures are chosen of the form 2 j (j = 1, . . .) and the sequence is determined according to the degree of symmetry in the distribution of the p i 's or when the minimum acceptance rate for swaps between adjacent temperatures is larger than a default threshold.

Example 6.5.1 Considering again the Old Faithful benchmark, we set this symmetry threshold to .1 and this acceptance threshold to 0.3. Using the same proposals as in Example 6.4.2 and N sim = 50, 000, the algorithm selects 4 temperatures, thus equal to 1, 2, 4, 8. Figure 6.11 demonstrates that the parallel tempering sampler visits all modes in the posterior distribution and that the mixing of the chains is greatly improved.

� mixtures Example 6.5.2 We now implement parallel tempering for a mixture of k = 3 components applied to a benchmark dataset from [START_REF][END_REF]]. This dataset is derived from an image of a car license plate, and made of 2625 observations. In [START_REF][END_REF]], a lack of label switching is observed when using a Gibbs sampler. Once again, this means each component can be estimated by its mean and standard deviation. The sample size is larger here and more likely to mixing problems. This is clearly exhibited in the six top plots of Figure 6.12 where the estimates provided for the three components are quite distinct. When implementing parallel tempering, the temperature increase stops when when all acceptance rates of swaps are above .4, meaning for this dataset 7 temperatures ranging from 1 to 64.

The six bottom plots of Figure 6.12 show that parallel tempering immensely improves the swaps between the posterior modes. The sample of �'s produced by parallel tempering visits a much larger region in (0, 2π), when compared with the highly peaked output of the original MCMC output.

The histograms in Figure 6.12 show that the posterior on p and η are now close to identical for each component. Two-dimensional plots also highlight this correct label switching behavior, which demonstrates better mixing and convergence of the produced chain. �

Conclusion

This paper has introduced a new parametrisation for mixtures of location-scale models. By constraining the parameters in terms of the global mean and global variance We have demonstrated that relatively standard simulation algorithms are able Chapter 6. Non-informative reparameterisations for location-scale mixtures to handle this new parametrisation and that they can manage the computing issues connected with label switching. In case of poor switching, we also established that parallel tempering can be easily implemented. As exhibited in the Ultimixt package, relabelling techniques are readily available.

While the extension to non-Gaussian cases with location-scale parameterisation (and beyond) is conceptually straightforward, considering this parameterisation in higher dimensions is delicate in terms of the covariance matrix. Indeed, even though we can easily set the global variance of the mixture as a parameter, reparameterising the component variances against this reference matrix remains an open question that we have not yet explored.

Chapter 7

Supplementary material:

Non-informative reparameterisations for location-scale mixtures Chapter 6 focuses on the reparametrisation of a mixture of Gaussian distributions with the purpose of using non-informative prior distributions on the parameters. We proceed here by some theoretical aspects, more data analyses and discuss the results.

The model of interest is a convex combination of the univariate Gaussian distributions defined by

f (x|µ i , σ i ) = k � i=1 p i N (x|µ i , σ i )
A feature of this combination of densities is that it allows to produce a probability density function because of preserving the properties of non negativity and integrating to 1. Multimodality of the produced density is another property which causes the "label switching" issue in the Bayesian analysis of the model as mentioned before while conditions for the number of modes have been explored by [Robertson 1969] and [Behboodian 1970]. As Bayesian methods enable the uncertainty in the model parameters to be directly quantified by examining the posterior distribution, they are useful for fitting these models to data. Despite that the mixture distributions have a range of applications, making an objective choice of prior for the component parameters is difficult in the case where no information is available to determine a subjective prior. Basically, assigning independent improper non-informative priors to the parameters of the mixture components results in improper posterior distribution as shown by [START_REF] Editeurs | [END_REF]] which unable these prior to be used for the mixtures. This difficulty motivated the idea of shifting the parameters of ith mixture component to two variability parameters α i and τ i by defining linear functions µ i = µ + σα i , σ i = στ i based on µ and σ acting as the intercept and slope of these linear equations. This change leads to a proper posterior derived from a Jeffreys prior for the global parameters of the mixtures, as demonstrated in Chapter 6. Since σ is positive, both µ i and σ i are increasing with respect to the values of α i and τ i . In a special case where (α i , τ i ); i = 1, . . . , k converges toward (0, 1), the mixture is Chapter 7. Supplementary material: Non-informative reparameterisations for location-scale mixtures transformed to a simple normal distribution with mean and standard deviation µ and σ, respectively.

Spherical coordinate concept

The modifications α i = σγ i/ √ p i and τ i = ση i/ √ p i bring about a hypersphere and hyperplane equations due to the constraints obtained from the mean and variance of the population that drive the resulting mixture model 6.3 more compact in terms of specifying the prior distribution for the resulting variability parameters. In addition, from the intersection of hypersphere � k i=1 γ 2 i = ϕ 2 centered at the origin and the hyperplane � k i=1 √ p i γ i = 0 which also passes through the origin, we deduce that γ i s belong to a circle of radius ϕ centered at the origin, as mentioned before. For example if we consider k = 3, ϕ = 0.5 and p = (0.35, 0.25, 0.4), we will have HP 3 : 0.59γ 1 + 0.5γ 2 + 0.63γ 3 = 0; HS 3 : γ 2 1 + γ 2 2 + γ 2 3 = 0.25, and a 3-dimensional graphical representation of HP 3 ∩ HS 3 is shown in Figure 7.1 which illustrates that the intersection is precisely a set of points in hyperplane at a distance of ϕ from the origin. For any k, here we represent the spherical coordinate of γ's obtained in Chapter 6 in more details. Suppose that

HP k : γ 1 √ p 1 + γ 2 √ p 2 + . . . + γ k √ p k = 0 HS k : γ 2 1 + γ 2 2 + . . . + γ 2 k = ϕ 2 Where � n i=1 p i = 1. From HP k , two vectors ( √ p 1 , . . . , √ p k ); Γ = (γ 1 , . . . , γ k ) 7.1. Spherical coordinate concept 129 are orthogonal. Let � k = ( √ p 1 , . . . , √ p k )/ � � � � k � i=1 p i = ( √ p 1 , . . . , √ p k ) (7.1)
be a unit-length vector. Since E denotes any Euclidean space of finite dimension k, the hyperplane HP k has dimension k -1 and we can find an orthonormal basis for HP k as (� 1 , . . . , � k-1 ) orthogonal to � k , [Gallier 2011]. An orthonormal basis including � k based on the orthogonal basis { �1 , . . . , �k } is given by

� 1 = (-√ p 2 , √ p 1 ,0,...,0,0) / √ p 1 +p 2 � 2 = (-√ p 1 p 3/ √ p 1 +p 2 ,- √ p 2 p 3/ √ p 1 +p 2 , √ p 1 +p 2 ,0,...,0) / √ p 1 +p 2 +p 3 � 3 = (-√ p 1 p 4/ √ p 1 +p 2 +p 3 ,- √ p 2 p 4/ √ p 1 +p 2 +p 3 ,- √ p 3 p 4/ √ p 1 +p 2 +p 3 , √ p 1 +p 2 +p 3 ,0,...,0) / √ p 1 +p 2 +p 3 +p 4 . . . � k-1 = (-√ p 1 p k/ √ � k-1 i=1 p i ,- √ p 2 p k/ √ � k-1 i=1 p i ,...,- √ p k-1 p k/ √ � k-1 i=1 p i , � � k-1 i=1 p i ) / � � k i=1 p i � k = ( √ p 1 , . . . , √ p k )
We can easily show that for all i ∈ 1, . . . , k; � i has unit-length and for all i � = j, dot product of � i and � j is zero. Using this orthonormal basis, any point on the hyperplane can therefore be expressed as

Γ = b 1 � 1 + b 2 � 2 + . . . + b k � k which gives γ 1 = -b 1 √ p 2/ √ p 1 +p 2 -b 2 √ p 1 p 3/ √ p 1 +p 2 √ � 3 i=1 p i -. . . -b k-1 √ p 1 p k/ � � k-1 i=1 p i � � k i=1 p i + b k √ p 1 γ 2 = b 1 √ p 1/ √ p 1 +p 2 -b 2 √ p 2 p 3/ √ p 1 +p 2 √ � 3 i=1 p i -. . . -b k-1 √ p 2 p k/ � � k-1 i=1 p i � � k i=1 p i + b k √ p 2 γ 3 = b 2 √ p 1 +p 2/ √ p 1 +p 2 +p 3 -b 3 √ p 3 p 4/ √ � 3 i=1 p i √ � 4 i=1 p i -. . . -b k-1 √ p 3 p k/ � � k-1 i=1 p i � � k i=1 p i + b k √ p 3 . . . γ k-1 = b k-2 � � k-2 i=1 p i/ � � k-1 i=1 p i -b k-1 √ p k-1 p k/ � � k-1 i=1 p i � � k i=1 p i + b k √ p k-1 γ k = b k-1 � � k-1 i=1 p i/ � � k i=1 p i + b k √ p k .
(γ 1 , . . . , γ k ) belongs to both HP k and HS k and thus replacing it in hyperplane results in Chapter 7. Supplementary material: Non-informative reparameterisations for location-scale mixtures

√ p 1 γ 1 + √ p 2 γ 2 + . . . + √ p k γ k = -b 1 √ p 1 p 2/ √ p 1 +p 2 + . . . + b k p 1 + b 1 √ p 1 p 2/ √ p 1 +p 2 + . . . + b k p 2 + . . . -b k-1 p k-1 √ p k/ � � k-1 i=1 p i � � k i=1 p i + b k p k-1 + b k-1 √ p k � � k-1 i=1 p i/ � � k i=1 p i + b k p k
and by canceling positive and negative similar terms we obtain � k i=1 b k p i = 0 that ends up with b k = 0. In this case, replacing γ i 's in HS k leads to a hypersphere of radius

ϕ in k -1-dimensional Euclidean space b 2 1 + b 2 2 + . . . + b 2 k-1 = ϕ 2
and thus any reparametrization of this object such as spherical coordinate system may be considered.

In the special case of k = 2, the orthonormal basis is defined by two following vectors

� 1 = (- √ p 2/ √ p 1 +p 2 , √ p 1/ √ p 1 +p 2 ); � 2 = ( √ p 1 , √ p 2 )
where p 1 = p; p 2 = 1p and we can therefore write

(γ 1 , γ 2 ) = b 1 � 1 + b 2 � 2 = (-b 1 � 1 -p + b 2 √ p, b 1 √ p + b 2 � 1 -p)
Replacing γ 1 and γ 2 above in both

HP 2 : √ pγ 1 + � 1 -pγ 2 = 0; HS 2 : γ 2 1 + γ 2 2 = ϕ 2 yields -b 1 � p(1 -p) + b 2 p + b 1 � p(1 -p) + b 2 (1 -p) = 0; b 2 = 0 b 2 1 (1 -p) + b 2 1 p = ϕ 2 ; b 2 1 = ϕ 2
b 2 1 = ϕ 2 can be considered as a sphere in 1-dimensional Euclidean space that represents a pair of points {-ϕ, +ϕ} which is the boundary of a line segment (a part of a line that is bounded by two distinct end points, and contains every point on the line between its endpoints). We can therefore rewrite γ i 's as

(γ 1 , γ 2 ) = (±ϕ � 1 -p, ±ϕ √ p).
Since η 2 1 + η 2 2 = 1ϕ 2 , spherical coordinate representation of η 1 and η 2 will be

η 1 = � 1 -ϕ 2 cos(ξ); η 2 = � 1 -ϕ 2 sin(ξ)
where ξ ∈ [0, π /2]. In example 6.4.1, the analyses related to 50 data points simulated from the mixture 0.65N (-8, 2) + 0.35N (-0.5, 1) are based on the case where a Dirichlet prior Dir(0.5, 0.5, 0.5) is assigned to (ϕ 2 , η 2 1 , η 2 2 ). It means that only the component means are expressed in spherical coordinates. Here, we reanalyze this model by considering the spherical coordinates above for η i 's. In this case, we place a beta prior on the parameter ϕ 2 with the same hyper parameters 0.5, 0.5. Figure 7.2 shows that the estimates of the marginal posterior distributions of means and standard deviations are symmetric and each µ i (σ i ) is very similar to one another due to the label switching phenomenon. As the component parameters are not identifiable marginally, estimating them on the basis of these MCMC output is not straightforward and we thus revert to the k-means clustering algorithm. The procedure is implemented by using the package Ultimixt in Chapter 8. The estimations of the component parameters displayed in Table 7.1 are quite similar to those of the true model while the calibration of the proposal scale results in the acceptance rates close to the optimal. On the other hand, since the estimate of the predictive density based on the MCMC output does not depend on the labelling of the components, Figure 7.3 shows that the estimated mixture is very smooth and unaffected by the label switching. Note that the black line is the estimate of the density that is obtained by averaging the simulated densities over the last 500 iterations. reparameterisations for location-scale mixtures The results of the example above illustrates that when k = 2, the Bayes estimates of the mixture parameters based on the spherical coordinate of η i 's are identical to the case where (ϕ 2 , η 2 1 , η 2 2 ) is supposed to be from Dirichlet distribution. In both cases, the parameters are accurately estimated as long as the convergence towards the stationary distribution is achieved as shown in Figure 7.2.

In the following, we apply the reparametrisation of the mixture distribution based on the spherical coordinate of the component parameters for some other datasets and summarize the resulting Bayesian analyses by implementing the functions of Ultimixt package, Chapter 8.

Data analyses

Two datasets Acidity dataset and Enzyme dataset were initially used by [START_REF] Richardson | [END_REF]] while Fishery dataset and Darwin's dataset are taken from [Frühwirth-Schnatter 2006]. For all datasets, we run the MCMC algorithm with 10 4 iterations and in Figures 7.4, 7.5, 7.6 and 7.7, the predictive mixture densities are computed once for the last 500 iterations another by averaging over 10 4 iterations (dark line in the figures).

Acidity data

The Acidity dataset is related to an acidity index measured in a sample of 155 lakes in the Northeastern United States. A histogram of the data points is shown in 7.4 and so a mixture of 2 components is well suited to model the data. In Figure 7.4 the histogram of the data is overlaid with the predictive density estimate obtained by applying Ultimixt package 8 that indicates the model represents the data well, with no need for any more components.

Enzyme data

Enzyme dataset concerns the distribution of enzymatic activity in the blood, for an enzyme involved in the metabolism of carcinogenic substances, among of group of 245 unrelated individuals with the purpose of identifying subgroups of slow or fast metabolizers as a marker of genetic polymorphism in the general population. This dataset has been reanalyzed by [van Havre 2014] who compared the probabilities that this data can be modeled by two or three components. [van Havre 2014] showed that fitting a mixture model with two components to the data is more likely than three components. We therefore analyze Enzyme dataset with our reparametrized mixture distribution by considering k = 2. The predictive density estimates shown in Figure 7.5 illustrate a good fitting of the mixture with two components to the data.

Darwin's data

Darwin's data consists of 15 observations of differences in heights between pairs of self-fertilized and cross-fertilized plants grown under the same condition. The histogram of the data overlaid with the predictive density estimates based on fitting a mixture of two normal components to the data is shown in Figure 7.7. Despite that the sample size is small, the mixture estimate fits very well the data. Table 7.2 shows the point estimates of the component parameters including mean and median for the datasets Acidity, Enzyme and Darwin. The table reveals very negligible difference between the estimates based on mean and median of the posterior draws and the errors between the point estimates of µ and x are zero for all datasets while for σ the error between the median of the posterior draws and the standard deviations of datasets are negligible (the error of posterior estimates based on the average of draws is also zero for all datasets). reparameterisations for location-scale mixtures

Fishery data

We fit a Gaussian mixture model on the Fishery data which consists of data on the lengths of 256 snappers. The heterogeneity in the data comes from the different age groups a fish might belong to depending if it comes from the current year's spawning or the previous, and so on.

Overfitting the Fishery dataset produces two possible alternate configurations with four or three components and we therefore analyze this data once by considering k = 4 another by modeling a mixture of k = 3 components. Table 7.3 displays the posterior parameters describing each configuration. When k = 4, we obtain two almost similar component means µ 3 = 3.6 and µ 4 = 3.2 and since one of them has an estimated weight close to zero this case can be considered as overfitted model. The Bayes estimates of the parameters of a mixture with three components exhibit a better performance as shown on the right of Table 7.3. We note here that the point estimates of global parameters µ and σ in the case where k = 4 are almost identical to the ones for k = 3. This indicates that the change of the component number of the mixture model does not impact the Bayesian inference of these two parameters. The radius ϕ decreases when k diminishes from 4 to 3 and since σ i 's are expressed by spherical coordinate with radius (1ϕ 2 ) 1 /2 , the component standard deviations naturally increase in comparison with the case of k = 4.

Figure 7.6 shows the posterior predictive densities based on MCMC output for the case where k = 4 that is almost identical to the one for k = 3. This illustrates that the estimated mixtures are unaffected by the overfitting phenomenon as mentioned before for Galaxy data in Chapter 6. Note that in both Tables 7.2 and 7.3, for the global parameters, the means and medians are computed based the obtained MCMC draws. These summary statistics are considered as the posterior estimations while for the angular and component-wise parameters the estimates are obtained by applying k-means clustering algorithm on the related posterior samples.

Parallel tempering algorithm

Last part of Chapter 6 deals with the analyses based on parallel tempering method. We showed that this method greatly improves mixing of MCMC chains. In some examples, we have seen that for a sample of size large enough, multimodality of the mixture model causes mixing problem or eventuates a good fit of one of the components to the extent that it is difficult for the Metropolis-within-Gibbs sampler to escape. In this case, using parallel tempering or replica exchange is suggested [START_REF] Gill | [END_REF]] which can be tracked down in a paper written by [START_REF] Swendsen | [END_REF]]. The applications of this method has used not only for problems in statistical physics but in chemistry, biology, engineering and materials science [Earl 2005]. The important reparameterisations for location-scale mixtures feature is preventing the Markov chain from sticking in minor modal areas for long periods of time. Parallel tempering MCMC as a method for generating candidate samples from all over a distribution is designed to overcome low probability regions between areas of importance. This means that a temperature parameter could be used to flatten out the target distribution such that the more temperature raised the more distribution flattens out and this makes the random chain more likely to mix quickly, for that temperature [START_REF] Lewandowski | [END_REF][START_REF] Gill | [END_REF], Earl 2005, Neal 1996, Wang 2011[START_REF] Li | [END_REF][START_REF] Swendsen | [END_REF]].

The method is that the target distribution is transformed to the Boltzmann distribution for a given temperature which is called "replicas". In order to simulate parameter θ from a posterior distribution non-standard f (θ|x), in temperature 1/β, the replicas is defined as f β (θ|x) = f (θ|x) β where 0 ≤ β ≤ 1. This means that the chains can be constructed from "tempered" versions of the target of interest by raising it to a power between 0 and 1, with 0 corresponding to a complete flattening of the distribution, and 1 corresponding to the desired target [START_REF] Altekar | [END_REF], Geyer 2011, Geyer 1991]. [START_REF] Miasojedow | [END_REF]] argue that for a target distribution f (), tempering of f () often provides better mixing within modes of the target distribution. However, this method is often more effective than the non-tempered approach [Hamze 2010].

Having this expression for the target distribution, we run the basic MCMC algorithm on each distribution and in each iteration, the current samples are considered probabilistically for exchanges between different temperature levels with probability α h as shown in parallel tempering algorithm below. Note that here only pairs between neighboring temperatures are considered for swapping because the chances of accepting a trade are more likely to be higher. Before performance of parallel tempering MCMC becomes optimal, we should tune the number of replicas h and their temperatures which are not actually evident and besides this, the simulation of multiple chains does increase the computation time. Several suggestions for the number of replicas and temperature of the replicas have been offered. An example is a geometric progression of temperatures [START_REF] Miasojedow | [END_REF], Earl 2005]. [START_REF] Miasojedow | [END_REF]] propose an adaptive algorithm in order to tune the temperature schedule and the parameters of the random-walk Metropolis kernel. Here, we consider that the proposal scale is automatically calibrated using adaptive Metropoliswithin-Gibbs algorithm for the target distribution according to the optimal acceptance rates as explained before and it is simultaneously used for the proposal distribution of all tempered versions of the target distribution. The number of temperatures is automatically chosen according to the degree of the symmetry of the generated samples of p to 0.5 in the case where the number of mixture components is k = 2. The number of temperatures is also determined according to the acceptance rate of the swaps between the neighboring temperatures, ar swap .

Let θ = (µ, σ, p, ϕ, ξ, �), ε θ be the scale of the proposal distribution in Metropoliswithin-Gibbs step. We suppose that δ 1 is the symmetry threshold and δ 2 acceptance rate of swaps threshold. We therefore define parallel tempering algorithm in 7.8. Note that if R 1 and R 2 indicate the ratios of the simulated samples for p, greater Chapter 8. Ultimixt package there does not exist a well-accepted reference Bayesian approach to such models. One reason for the difficulty is the general prohibition against using improper priors (Fruhwirth-Schnatter, 2006) due to the ill-posed nature of such statistical objects. [START_REF] Kamary | Non-informative reparameterisation 8.2. K.MixReparametrized function References[END_REF] took advantage of a mean-variance reparametrisation of a Gaussian mixture model to propose improper but valid reference priors in this setting. This R package implements the proposal and computes posterior estimates of the parameters of a Gaussian mixture distribution. The approach applies with an arbitrary number of components. The Ultimixt R package contains an MCMC algorithm function and further functions for summarizing and plotting posterior estimates of the model parameters for any number of components. Details Package: Ultimixt Type: Package Version: 1.0 Date: 2015-10-30 License: GPL (>=2.0) Beyond simulating MCMC samples from the posterior distribution of the Gaussian mixture model, this package also produces summaries of the MCMC outputs through numerical and graphical methods.

Note: The proposed parameterisation of the Gaussian mixture distribution is given by

f (x|µ, σ, p, ϕ, �, ξ) = k � i=1 p i f (x|µ + σγ i / √ p i , ση i / √ p i )
under the non-informative prior π(µ, σ) = 1/σ. Here, the vector of the

γ i = ϕΨ i � �, p � i
's belongs to an hypersphere of radius ϕ intersecting with an hyperplane. It is thus expressed in terms of spherical coordinates within that hyperplane that depend on k-2 angular coordinates � i . Similarly, the vector of

η i = � 1 -ϕ 2 Ψ i � ξ � i
's can be turned into a spherical coordinate in a k-dimensional Euclidean space, involving a radial coordinate � 1ϕ 2 and k -1 angular coordinates ξ i . A natural prior for � is made of uniforms, � 1 , . . . , � k-3 ∼ U [0, π] and � k-2 ∼ U [0, 2π], and for ϕ, we consider a beta prior Beta(α, α). A reference prior on the angles ξ is (ξ 1 , . . . , ξ k-1 ) ∼ U [0, π/2] k-1 and a Dirichlet prior Dir(α 0 , . . . , α 0 ) is assigned to the weights p 1 , . . . , p k of location-scale mixtures. arXiv.

See Also Ultimixt

Examples data(faithful) xobs=faithful[,1] estimate=K. MixReparametrized(xobs,k=2,alpha0=.5,alpha=.5,Nsim=1e4) 8.2 K.MixReparametrized function

K.MixReparametrized

Sample from a Gaussian mixture posterior associated with a noninformative prior and obtained by Metropolis-within-Gibbs sampling Description This function returns a sample simulated from the posterior distribution of the parameters of a Gaussian mixture under a non-informative prior. This prior is derived from a mean-variance reparameterisation of the mixture distribution, as proposed by Kamary et al. (2015). The algorithm is a Metropolis-within-Gibbs scheme with an adaptive calibration of the proposal distribution scales. Adaptation is driven by the formally optimal acceptance rates of 0.44 and 0.234 in one and larger dimensions, respectively (Roberts et al.,1997). This algorithm monitors the convergence of the MCMC sequences via Gelman's and Rubin's (1992) criterion.

Value

The output of this function is a list of the following variables, where the dimension of the vectors is the number of simulations: Details This function outputs posterior point estimates for all parameters of the mixture model. They mostly differ from the generaly useless posterior means. The output summarizes unimodal MCMC samples by computing measures of centrality, including mean and median, while multimodal outputs require a pre-processing, due to the label switching phenomenon (Jasra et al., 2005). The summary measures are then computed after performing a multi-dimensional k-means clustering [START_REF] Hartigan | A K-means clustering algorithm[END_REF] following the suggestion of Fruhwirth-Schnatter (2006). Author(s) Kaniav Kamary algorithm determines the scales such a way that the chain does not move neither too slowly (most of the proposals are accepted) nor very quickly (says most of the proposals will usually be rejected). This method asks the computer to automatically learn better parameter values while an algorithm runs. Under some conditions this method updates the scale of each proposal (different for each parameter) at each iteration with the intention of finding the best value. In other words, the scale of the proposal distribution is automatically calibrated according to the optimal acceptance rate (the fraction of the proposed moves which are accepted) and with a factor that decreases to 0 in such a way that the convergence conditions are held [Robert 2009a, Rosenthal 2011, Roberts 2001, Roberts 2009, Roberts 1997]. So when the acceptance rate on batches of 50 iterations (by default) is too high, the proposal variance is automatically increased whereas in the case of the acceptance rate too small, it will decrease by adding or subtracting min(0.01, 1/ √ t) after the t th iterations. [Roberts 1997] show that for a random walk Metropolis-Hastings with a d-dimensional target distribution which consists of i.i.d components, when the number of the parameters to simulate tends to ∞ the optimal acceptance rate is 0.234 and for a one dimensional problem, the optimal acceptance rate is approximately 0.44. With this description, the adaptive Metropolis-within-Gibbs algorithm can be summarized as follows: For simplicity's sake, let θ = (µ, σ, p, ϕ, ξ, �) be the mixture parameters and ε θ = (ε µ , ε σ , ε p , ε ϕ , ε ξ , ε � ) be the scale of the proposal distribution q θ () in Metropolis-within-Gibbs step from which the proposals of θ are generated. We also create an associated variable log(ε θ ) giving the logarithm of the standard deviation to be used when simulating a proposal increment to parameter θ. For a total number of iterations T , suppose that ar opt denotes the optimal acceptance rate of the proposals and after the j th batch of 50 and l th batch of 500 iterations, n j and N l are supposed to be the number of times that the proposals have been accepted while ar j = n j/50 and AR l = N l/500 are the related acceptance rate, respectively, for j = 1, . . . , T /50 and l = 1, . . . , T /500. If T t denotes the number of the batch of 500 iterations from which we start testing AR Tt towards optimal acceptance rate with threshold values δ, the adaptive Metropolis-within-Gibbs algorithm will be defined in Figure 9.1.

In the programs of Ultimixt, T t = 30 by default and so after 15, 000 iterations, the algorithm starts comparing the acceptance rate of the last 500 iterations with the optimal acceptance rate. The calibration step is terminated when the resulting acceptance rate is located in a small neighborhood of the optimal acceptance rate. Note that the total number of MCMC iterations at most is T = 30, 000 by default in adaptive Metropolis-within-Gibbs step of K.MixReparametrized function.

Convergence monitoring of the chain provided for all the parameters of the mixture distribution by applying [Gelman 1992] criterion;

Description of implementation 151

Adaptive Metropolis-within-Gibbs algorithm for reparameterised mixture:

1 Initialize θ (0) , ε

θ ; d 1 = 0; d 2 = 0; j = 1; l = 1. 2 For t = 1, . . . , T , the update of θ (t) and ε (t) θ follows 2.1 Generate a proposal θ � ∼ q(•|θ (t-1) ) and update θ (t) against π(•|x). This step is started after calibrating the proposal scales and in each 1000 additional iterations, [Gelman 1992] criterion is automatically computed for the posterior draws of the component parameters, each one based on 4 chains produced in parallel. The simulation is stopped when this criterion is close to 1 for all produced chains. [Gelman 1992] criterion is a convergence monitoring diagnostic and this method allows us not to continue simulating the parameters after achieving the convergence.

Thus the Ultimixt package calls function gelman.diag() from package coda for the convergence monitoring step.

SM.MixReparametrized and SM.MAP.MixReparametrized are able to analyze the output of the MCMC simulations by numerical methods. For the unimodal terms such as the mean µ, the variance σ and the radius coordinate ϕ, the draws are regarded as posterior draws and averaging and calculating median over these draws will be considered as the point estimator.

This method is not satisfactory in the case where the label switching problem occurs in the posterior draws because of the multimodality of the posterior distribution as discussed before. The label switching occurs for the parameters such as component weights, means and variances p i , µ i and σ i ; i = 1, . . . , k. In this case, the function SM.MixReparametrized eliminates the resulting unidentifiability using a method of post processing the MCMC draws by imposing a restriction on the ordering of the mixture component parameters, as, e.g., the constraint that the weights of the component distributions are ascending. After that, a clustering procedure is applied to permute the MCMC draws of component means and standard deviations by applying a two dimensional k-means type algorithm with k! clusters. This method is a combination of model identification methods such as identifiability constrain and unsupervised clustering suggested in [Frühwirth-Schnatter 2006]. The simulated draws of the angles are also summarized by applying two one-line k-means algorithm with k -2 and k -1 clusters to a sample of size T (k -2) for � i ; i = 1, . . . , k -2 and a sample of size T (k -1) for ξ i ; i = 1, . . . , k -1, respectively. Note that in order to avoid the problem caused by unidentifiability, the number of the mixture components should initially be properly chosen.

Using Plot.MixReparametrized function, the MCMC results can be visually analyzed by using trace plots and by plotting the estimated densities for the mixture distribution over the draws. With the functions of Ultimixt package the user can comfortably reproduce some of the results presented in simulation study section 6.4 of Chapter 6.

Application

In the following we illustrate the use of Ultimixt on an example for a simulated dataset with 50 data points for which the histogram is shown in Figure 9.2. The sample mean and standard deviation are 4.06, 3.38, respectively. We therefore fit a Gaussian mixture model on the data by choosing k = 5. In order to sample from the posterior distribution of the mixture model parameters, we apply the function K.MixReparametrized. To do so, we have to specify the hyper parameters α 0 , α and the total number of MCMC iterations. We therefore choose α 0 = 0.5, α = 0.5 and 10 4 iterations. The output is firstly summarized by Plot.MixReparametrized function. The distribution of the draws of µ and σ shown in Figures 9.3 are centered on the empirical mean and standard deviation of the dataset. The posterior distribution of each component parameter is evaluated by plotting boxplot in Figure 9.3 that helps us to quickly examine the output of K.MixReparametrized from a graphical way. Figure 9.2 also shows the related estimate of the predictive density that is efficiently fitted the data.

The evolution of the scales of the proposal distributions in the Metropolis-within-Gibbs algorithm is one of the outputs of Plot.MixReparametrized function that is displayed in Figure 9.4. This figure shows that the update of the proposal scales is terminated for 300 batch of 50 iterations of the adaptive Metropolis-within-Gibbs step. This means that the acceptance rates attain the stability over the optimal one after 15000 MCMC iterations.

Mots Clés

Distribution de mélange, Loi a priori noninformative, Analyse bayésienne, A priori impropre, Choix du modèle bayésien, Méthodes de MCMC.

Abstract

In the special case of the Bayesian approach, the solution of model comparison is the Bayes factor. The Bayes factor is however problematic in some cases. The first part of this thesis deals with a general review on non-informative priors and demonstrating the overall stability of posterior distributions. In the second part, we consider a novel paradigm for Bayesian testing of hypotheses and Bayesian model comparison. The idea is to define an alternative to the traditional construction of posterior probabilities that a given hypothesis is true and to replace the original testing problem with estimation that focus on the probability weight of a given model within a mixture model. In the last part, we construct a reference Bayesian analysis of mixtures of Gaussian distributions by creating a new parameterization centered on the mean and variance of those models itself. This enables us to develop a genuine noninformative prior for Gaussian mixtures with an arbitrary number of components. This part of the thesis is followed by the description of R package, Ultimixt, which implements a generic reference Bayesian analysis of unidimensional mixtures of Gaussian distributions obtained by a location-scale parameterization of model.

Figure 2 . 1 :

 21 Figure 2.1: Logistic cdfs across a few thousand simulations from the normal prior, when using the prior selected by [Seaman III 2012] (left) and the prior defined as the G-prior(right)

Figure 2 . 2 :

 22 Figure 2.2: Posterior distributions of the logistic parameter α when priors are N (0, σ) for σ = 10, 25, 100, 900, based on 10 4 MCMC simulations.

  Posterior estimates of the logistic parameters under a g-prior, a flat prior and Jeffreys' prior for the banknote benchmark. Posterior means and standard deviations remain quite similar under all priors. All point estimates are averages of MCMC samples of size 10 4 .

Figure 2 . 3 :

 23 Figure 2.3: Posterior distributions of the parameters of the logistic model when the prior is N (0, 900 2 ), g-prior, flat prior and Jeffreys' prior, respectively. The estimated posterior distributions are based on 10 4 MCMC iterations.

  Figure 2.4 compares the lognormal prior with mean -1 and standard deviation 0.5 and the gamma distribution G(4, 1).

Figure 2 . 4 :

 24 Figure 2.4: Comparison of lognormal and gamma priors for the standard deviation of the regression slope.

Figure 2 . 5 :

 25 Figure 2.5: Estimated posterior densities of the regression intercept (top left), slope (top right), standard deviation of the intercept (down left) and standard deviation of the slope (down right), respectively for 4 different normal regressions. All estimates based on 10 5 iterations simulated from Metropolis-withing-Gibbs algorithm.

Figure 2 . 6 :

 26 Figure 2.6: Priors induced on the evenness index: Four Dirichlet prior are assigned to θ with hyperparameters all equal to 0.1, 0.25, 0.5, 1, based on 10 4 simulations.

Figure 2 . 7 :

 27 Figure 2.7: Estimated posterior densities of H(θ) considering sample sizes of 50, 100, 250, 1000, 10, 000. They correspond to the priors on θ shown in Figure 2.6 and are based on 10 4 posterior simulations. The vertical line indicates the mode of all posteriors when sample size is large enough.

Figure 2 . 8 :

 28 Figure 2.8: Jeffreys' prior and estimated posterior densities of H(θ) considering sample sizes 50, 100, 250, 1000, 10, 000. The posterior distributions are based on 10 4 posterior draws. The vertical line indicates the mode of the posterior density when the sample size is 10 4 .

Figure 3 . 1 :

 31 Figure 3.1: Simulation of posterior distribution of α and β with a multivariate normal random walk when the proposal scale τ takes values 0.1, 1, 5 and a flat prior is assigned to α and β. From top to bottom: Sequence of 10 4 iterations; Empirical autocorrelation; Histogram of the last 9000 iterations compared with the target density.

Figure 3 . 2 :

 32 Figure 3.2: Simulation of posterior distribution of α and β with a multivariate normal random walk when the proposal scale τ takes values 0.1, 1, 5 and a g-prior is assigned to α and β. From top to bottom: Sequence of 10 4 iterations; Empirical autocorrelation; Histogram of the last 9000 iterations compared with the target density.

  π

Figure 3 . 3 :

 33 Figure 3.3: Simulation of posterior distribution of α and β with a multivariate normal random walk when the proposal scale τ takes values 0.1, 1, 5 and a Jeffreys prior is assigned to α and β. From top to bottom: Sequence of 10 4 iterations; Empirical autocorrelation; Histogram of the last 9000 iterations compared with the target density.

Figure 3 . 4 :

 34 Figure 3.4: Simulation of posterior distribution of the parameters of the normal regression model when standard deviations of the intercept and slope have log normal distribution. From top to bottom: Sequence of last 9000 iterations; Empirical autocorrelation; Histograms. True values of the parameters are indicated at the top of sequence plots.

Figures 3 .

 3 Figures 3.4 and 3.5 give an assessment of the convergence of the algorithm and show that for both priors, log normal and gamma distributions, the distributions of the chains visually cover the whole support of the target distribution with sufficient regularity for 10 4 MCMC iterations. The autocorrelation plots show high mixing behavior of the chains and from the histograms, the distributions of the generated

Figure 3 . 6 :

 36 Figure 3.6: Empirical density of simulated draws from conditional posterior density of β11 and β12 based on last 9000 iterations when (Top) log normal prior and (Bottom) gamma prior are assigned to the hyper parameter s j 's. True values of the parameters are indicated at the top of each graph.

Figure 4 . 1 :

 41 Figure 4.1: Gibbs sequences (α t ) on the first component weight for the mixture model αN (µ, 1) + (1α)N (0, 1) for a N (0, 1) sample of size N = 5, 10, 50, 100, 500, 10 3 (from top to bottom) based on 10 5 simulations. The y-range range for all series is (0, 1).

Figure 4

 4 Figure 4.2: Metropolis-Hastings sequences (α t ) on the first component weight for the mixture model αN (µ, 1) + (1α)N (0, 1) for a N (0, 1) sample of size N = 5, 10, 50, 100, 500, 10 3 (from top to bottom) based on 10 5 simulations.The y-range for all series is (0, 1).

  Figure 4.4 gives an assessment of the convergence of the Metropolis-Hastings for λ and the mixture model weight α even if the sample size is very small (n=5).

Figure 4

 4 Figure 4.3: Example 4.3.1: Boxplots of the posterior means (wheat) of λ and the posterior medians (dark wheat) of α for 100 Poisson P(4) datasets of size n = 1000 for a 0 = .0001, .001, .01, .1, .2, .3, .4, .5. Each posterior approximation is based on 10 4 Metropolis-Hastings iterations.

Figure 4 .

 4 Figure 4.5 highlights the convergence of the posterior means and posterior medians of α as the sample sizes n increase for the same Poisson P(4) simulated samples.The sensitivity of the posterior distribution of α on the hyperparameter a 0 is clearly expressed by that graph. While all posterior means and medians converge to 1 in this simulation, the impact of small values of a 0 on the estimates is such that we consider values a 0 ≤ .1 as having too strong and too lengthy an influence on the posterior distribution to be acceptable.We can also compare the outcome of a traditional (albeit invalid, since relying on improper priors) Bayesian analysis with our estimates of α. Figure 4.6 shows how the posterior probability of model M 1 and the posterior median of α relate as the sample size grows to 1000. The shaded areas indicate the range of all estimates of α, which varies between .2 and .8 for a 0 = .5 and between 0 and 1 for a 0 ≤ .1. This

  Figure 4.5 highlights the convergence of the posterior means and posterior medians of α as the sample sizes n increase for the same Poisson P(4) simulated samples.The sensitivity of the posterior distribution of α on the hyperparameter a 0 is clearly expressed by that graph. While all posterior means and medians converge to 1 in this simulation, the impact of small values of a 0 on the estimates is such that we consider values a 0 ≤ .1 as having too strong and too lengthy an influence on the posterior distribution to be acceptable.We can also compare the outcome of a traditional (albeit invalid, since relying on improper priors) Bayesian analysis with our estimates of α. Figure 4.6 shows how the posterior probability of model M 1 and the posterior median of α relate as the sample size grows to 1000. The shaded areas indicate the range of all estimates of α, which varies between .2 and .8 for a 0 = .5 and between 0 and 1 for a 0 ≤ .1. This

Figure 4 . 5 :

 45 Figure 4.5: Example 4.3.1: Posterior means (sky-blue) and medians (grey-dotted) of the posterior distributions on α, displayed over 100 Poisson P(4) datasets for sample sizes from 1 to 1000. The shaded and dotted areas indicate the range of the estimates. Each plot corresponds to a Beta prior on α with parameter a 0 = .1, .2, .3, .4, .5 and each posterior approximation is based on 10 4 iterations.

  Figure 4.6: Example 4.3.1: Comparison between the ranges of P(M 1 |x) (red dotted area) and of the posterior medians of α for 100 Poisson P(4) datasets with sample sizes n ranging from 1 to 1000 and for several values of the hyperparameter a 0 .

  Figure 4.8 displays the range of the posterior means and medians of α when either a 0 or n varies, showing the

Figure 4 . 7 :

 47 Figure 4.7: Example 4.3.1: Boxplots of the posterior medians of α for 100 geometric Geo(.1) datasets of size n = 5, 50, 500. Boxplots are plotted using four beta priors for α with a0 = .1, .2, .4, .5. Each posterior approximation is based on 10 4 iterations.

Figure 4 . 8 :

 48 Figure 4.8: Example 4.3.2: Boxplots of the posterior means (wheat) and medians of α (dark wheat), compared with a boxplot of the exact posterior probabilities of M 0 (gray) for a N (0, 1) sample, derived from 100 datasets for sample sizes equal to 15, 50, 100, 500. Each posterior approximation is based on 10 4 MCMC iterations.

Figure 4 . 9 :

 49 Figure 4.9: Example 4.3.2: Plots of ranges of log(n) log(1 -E[α|x]) (gray color) and log(1p(M 1 |x)) (red dotted) over 100 N (0, 1) samples as sample size n grows from 1 to 500. and α is the weight of N (0, 1) in the mixture model. The shaded areas indicate the range of the estimations and each plot is based on a Beta prior with a 0 = .1, .2, .3, .4, .5, 1 and each posterior approximation is based on 10 4 iterations.

  Figure 4.11 displays the

Figure 4

 4 Figure 4.10: Example 4.3.2: (left) Posterior distributions of the mixture weight α and (right) of their logarithmic transform log{α} under a Beta B(a 0 , a 0 ) prior when a 0 = .1, .2, .3, .4, .5, 1 and for a normal N (0, 2) sample of 10 3 observations. The MCMC outcome is based on 10 4 iterations.

Figure 4

 4 Figure 4.11: Example 4.3.3 Posterior distributions of the N (µ, 1) component weight α under a Beta B(a 0 , a 0 ) prior (left) for a 0 = 1 /n, .1, .2, .3, .4, .5, 1 with 10 2 N (1, 1) observations and (right) for a 0 = .1 with n = 5, 10, 50, 100, 500, 10 3 N (1, 1) observations. In both cases each posterior approximation is based on 10 5 MCMC iterations.

Figure 4 .

 4 Figure 4.12: Example 4.3.4: Ranges of posterior means (skyblue) and medians (dotted) of the weight α of model N (θ, 1) over 100 N (0, .7 2 ) datasets for sample sizes from 1 to 1000. Each estimate is based on a Beta prior with a 0 = .1, .3, .5 and 10 4 MCMC iterations.

Figure 4 .

 4 Figure 4.13: Example 4.3.4: (left) Boxplot of the posterior means (wheat) and medians (dark wheat) of α, and of the posterior probabilities of model N (µ, 1) over 100 N (0, .7 2 ) datasets for sample sizes n = 10, 40, 100, 500; (right) averages of the posterior means and posterior medians of α against the posterior probabilities P(M 1 |x) for sample sizes going from 1 to 1000. Each posterior approximation is based on 10 4 Metropolis-Hastings iterations.

  (a) Pima dataset (b) Data from logistic (c) Data from probit

Figure 4 .

 4 Figure 4.14: Example 4.3.5: Posterior distributions of α in favor of the logistic model based on 10 4 Metropolis-Hastings iterations where a 0 = .1, .2, .3, .4, .5.

Figure 4 . 15 :

 415 Figure 4.15: Example 4.3.6: (top) Posterior medians of the true model weight over 5 values of a 0 = .1, .2, .3, .4, .5 for sample sizes ranging from 1 to 10 4 and (bottom) from 1 to 200. Case 1 (left) and Case 2 (right) correspond to common and independent parameterisations of the mixture components. Each approximation is based on 10 4 Gibbs iterations.

  2 |y, X] 0.9836 0.9104 0.8043 0.7190 0.5190 Case 2: E[α 2 |y, X] 0.9611 0.9018 0.7743 0.6780 0.3905 Table 4.3: Example 4.3.6: Posterior means of α 2 , weight of model M 2 , when the sample size is n=30based on a single dataset, they are conclusive about the performances of the mixture approach.

Figure 4 .

 4 Figure 4.16: Example 4.3.6: log(1 -E(α 2 |y, X)) and log(1 -P(M 2 |y, X)) (red lines) over logarithm of the sample size for a 0 = .1, .2, .3, .4, .5. Each posterior approximation is based on 10 4 iterations.

Figure 4 .

 4 Figure 4.17: Case study: (left) Normal (solid), Gumbel (dashed) and Logistic (dotted) densities with (0, 1) parameter; (bottom) Boxplots of the posterior distributions of the weights under the 3 scenario: truth = Normal (left panel), truth = Gumbel (middle panel), truth = logistic (right panel).

Figure 4 .

 4 Figure 4.18: Case study: Boxplots of the posterior distributions of the Normal weight α 1 under the two scenarii: truth = Normal (left panel), truth = Gumbel (right panel), a 0 =0.01, 0.1, 1.0, 10.0 (from left to right in each panel) and n = 1, 000 observations.

Figure 4 .

 4 Figure 4.19: Case study: Boxplots of the posterior distributions of the Normal weight α 1 under the two scenarii: truth = Normal (left panel), truth = Gumbel (right panel), a 0 =0.01, 0.1, 1.0, 10.0 (from left to right in each panel) and n = 10, 000 simulated observations.

Figure 4 .

 4 Figure 4.20: Case study: (top panel) Histograms of the hazard (left) and log hazard

Figure 5 . 1 :

 51 Figure 5.1: Mixture weights distribution 5.1: (Left) Probability density function of beta prior for the weight of the mixture model. (right) Posterior probability of the mixture weights as a function of n 1 when a 0 = .5 and n 1 + n 2 = 100.

Figure 5 . 2 :

 52 Figure 5.2: Poisson versus geometric 5.2.1: Sequence of λ[t], p[t] and α[t] simulated by GSmix function with 10, 000 iterations for a dataset of size 30 from P(.5) when a 0 = .1.

Figure 5 . 3 :

 53 Figure 5.3: Poisson versus geometric 5.2.1: Sequence (Top), empirical autocorrelations using acf function in R (Middle) and histograms (Bottom) of λ[t], p[t] and α[t] with 10, 000 MCMC iterations for datasets of sizes 50, 1000 from P(.64) and P(10) when a 0 = .1.

Figure 5 . 4 :

 54 Figure 5.4: Poisson versus geometric 5.2.1: Sequence (Top), empirical autocorrelations using acf function in R (Middle) and histograms (Bottom) of λ[t], p[t] and α[t] with 10, 000 MCMC iterations for datasets of sizes 400, 600 from Geo(.4) and P(.59) when a 0 = .1.

Figure 5 . 5 :

 55 Figure 5.5: Poisson versus geometric 5.2.1: Evolution of α over the sample size when data is simulated from (left) Poisson distribution with parameter λ taking values from .5 to 10; (right) Geometric distribution with parameter p taking values from .06 to 1. Each estimation is based on 10, 000 MCMC iterations.

Figure 5 . 6 :

 56 Figure 5.6: Poisson versus geometric 5.2.2: Marginal posterior distributions of λ, p and α obtained by running MCMC algorithm for 50 times with a 0 = 0.1, 0.2, 0.3, 0.4, 0.5 for (Top) a sample of size 50 simulated from P(.64) when the true value of (β 1 , β 2 , δ 1 , δ 2 ) is (32, 50, 50, 32); (Bottom) a sample of size 500 simulated from Geo(.49) when the true value of (β 1 , β 2 , δ 1 , δ 2 ) is (505, 490, 490, 505). Each density is based on 10, 000 MCMC iterations.

  from N (θ, 2), we analyze the output of both algorithms. Comparison between the estimated marginal posterior distributions of θ derived from the outputs of both Metropolis-Hastings and Gibbs sampler is shown in Figures 5.7 and 5.8. Both methods accurately estimate the common location parameter of the normal distribution for any value of a 0 while the resulting posterior densities of α are identical and in favor of the true model in both cases even for samples of size 10.

Figure 5

 5 Figure 5.7: N (θ, 1) versus N (θ, 2) 5.3: Marginal posterior densities of θ (Top) and α (Bottom) obtained by running (dashed lines) Gibbs sampling algorithm and (solid lines) Metropolis-Hastings algorithm with different values of a 0 for samples of sizes 10, 510 and 1000 simulated from N (-10.15, 1), N (-0.03, 1) and N (9.97, 1), respectively. Each density is based on 10, 000 MCMC iterations and the vertical dotted line in θ plots corresponds to the true value.

Figure 5 . 8 :

 58 Figure 5.8: N (θ, 1) versus N (θ, 2) 5.3: Marginal posterior densities of θ (Top) and α (Bottom) obtained by running (dashed lines) the Gibbs sampling algorithm and (solid lines) the Metropolis-Hastings algorithm with different values of a 0 for samples of sizes 10, 510 and 1000 simulated from N (-9.72, 2), N (0.02, 2) and N (10.01, 2), respectively. Each density is based on 10, 000 MCMC iterations and the vertical dotted line in θ plots corresponds to the true value.

Figure 5

 5 Figure 5.9: N (0, 1) versus N (µ, 1) 5.4: Sequence, empirical autocorrelations using acf function in R and histograms of µ[t] and α[t] simulated by the Metropolis-within-Gibbs algorithm with 10, 000 iterations for: (Top) two datasets of sizes 5 (left) and 250 (right) from N (-9, 1) and N (1.95, 1), respectively; (Bottom) two datasets of sizes 850 (left) and 50, 000 (right) simulated from standard normal distribution. a 0 = .1.

  estimation model whatever the sample size is. By comparing the posterior distributions of α related to the datasets D.1, D.2, D.3 with N = 5, 15, 25 in the case where the datasets are from standard normal distribution, we can see that the variation of the posterior draws is high, especially for the sample of small size.

Figure 5 .

 5 Figure 5.10: N (0, 1) versus N (µ, 1) 5.4: Posterior distribution of α[t], the weight of the standard normal in the mixture model, under a beta prior with parameter a 0 = 0.1 for (left) 16 standard normal datasets and (right) 16 datasets simulated from N (µ, 1) when each posterior approximation is based on 10 4 MCMC iterations.Table 5.4 lists more details about the datasets.

Figure 5 .

 5 Figure 5.11: N (µ, 1) versus L(µ, √ 2) 5.5: (Top) Sequence of µ[t] and (Bottom) empirical densities of the posterior draws of α for 21 datasets from N (0, 1) (Left) and from L(0, √ 2) (Right), based on 10 4 iterations of MCMC algorithm when a 0 = 0.5.

Figure 5 . 12 :

 512 Figure 5.12: Logistic versus probit 5.6: (Top) Sequence of θ[t]; (Center) Autocorrelation over 10 4 iterations; Histogram over the last 9000 iterations. a 0 = .5.

Figure 5 . 13 :

 513 Figure 5.13: Logistic versus probit 5.6: Posterior distributions of α in favor of the logistic model based on 10 4 MCMC iterations where a 0 = .1, .2, .3, .4, .5. (Left) Two datasets of sizes 10 4 , 20, 000 simulated from logistic model; (Right) Two datasets of sizes 10 4 , 20, 000 simulated from probit model that are related to the analyses shown in Table 5.6.

Figure 5 .

 5 Figure 5.14: Case 1. Caterpillar dataset 5.7: (Top) Sequences of 1e5 Gibbs sampler iterations; Empirical autocorrelations using acf R function; Histograms of the last 90, 000 iterations. (Bottom) Histograms of the posterior distributions of α 1 , . . . , α 15 based on 10 5 MCMC iterations when a 0 = .1.

Figure 5 .

 5 Figure 5.15: Case 2. Caterpillar dataset 5.7: Sequences, autocorrelations and histograms of 10 4 Metropolis-within-Gibbs iterations for the regression coefficients β j of the model M j ; j = 1, . . . , 7 shown in Table5.8 when a 0 = .5.

Figure 5 .

 5 Figure 5.16: Case 2. Caterpillar dataset 5.7: Sequences, autocorrelations and histograms of 10 4 Metropolis-within-Gibbs sampler iterations for the regression coefficients β j of the model M j ; j = 8, . . . , 15 shown in Table5.8 when a 0 = .5.

Figure 6 . 2 :

 62 Figure 6.2: Example 6.4.1: Kernel estimates of the posterior densities of the parameters µ, σ, p, ϕ, η i , based on 10 parallel MCMC chains for Proposal 1 and 2 10 5 iterations, based on a single simulated sample of size 50. The true value of (ϕ, η 1 , η 2 ) is (0.91, 0.16, 0.38).

Figure 6 . 3 :

 63 Figure 6.3: Example 6.4.1:Comparison between MCMC samples from our Metropolis-within-Gibbs algorithm using Proposal 1 (solid line) or Proposal 2 (dashed line), with 90, 000 iterations and the same sample as in Figure6.2. The true value of (ϕ, η 1 , η 2 ) is (0.91, 0.16, 0.38).

Figure 6

 6 Figure 6.4: Old Faithful dataset (Example 6.4.2): Posterior distributions of the parameters of a two-component mixture distribution based on 50, 000 MCMC iterations.

Figure 6 . 5 :

 65 Figure 6.5: Example 6.4.3: (Left) Evolution of the sequence (� (t) ) and (Right) histograms of the simulated weights based on 10 5 iterations of an adaptive Metropoliswithin-Gibbs algorithm with independent proposal on �.

Figure 6 . 7 :

 67 Figure 6.7: Example 6.4.3: Estimated marginal posterior densities of component means and standard deviations, based on 10 5 MCMC iterations.

  simulated 20 samples of size 250. Calibration of the random walks is achieved after 10 4 for almost all samples.

Figure 6 . 8 :

 68 Figure 6.8: Galaxy dataset: (Left) Representation of 500 MCMC iterations as mixture distributions with the overlaid average curve for k = 6 components (dark line); (Right) mixture density estimate based on 15, 000 MCMC iterations for k = 6 components.

Figure 6 . 9 :

 69 Figure 6.9: Galaxy dataset: (left) Representation of 500 Metropolis-within-Gibbs iterations for the mixture estimation and the overlay curve (dark line) obtained by averaging over iterations; (right) The mixture density estimate to histogram of dataset computed by averaging over 50, 000 MCMC iterations.

  Figure 6.10: Galaxy dataset: Empirical prior densities based on 10 4 samples for µ 1 , σ 1 and p 1 when (Top) k = 6 and (Bottom) k = 4. For the proposed prior (solid lines), the priors induced are π(µ 1 ) ∝ π(σγ 1 / √ p 1 ) and π(σ 1 ) ∝ π(ση 1 / √ p 1 ).

Figure 6 .

 6 Figure 6.11: Faithful dataset 6.4.2: Posterior distribution of the mixture distribution parameters and comparison between the lowest and highest temperatures (target distribution and f (x|θ) 1/8 ) of parallel tempering outputs based on 50, 000 iterations.

Figure 6 . 12 :

 612 Figure 6.12: Licence dataset (Example 6.5.2): Comparison between Metropoliswithin-Gibbs and parallel tempering outputs: The distributions of the samples of 10 4 last points and corresponding 2 × 2 plots.

Figure 7 . 1 :

 71 Figure 7.1: Intersection between 3-dimensional hyperplane and hypersphere.

Figure 7 . 2 :

 72 Figure 7.2: Mixture of two normal distributions 7.1: Estimated marginal posterior densities of component means µ 1 , µ 2 and standard deviations σ 1 , σ 2 based on 10 5 MCMC iterations.

Figure 7 . 3 :

 73 Figure 7.3: Mixture of two normal distributions 7.1: Representation of 500 MCMC iterations as mixture distributions with the overlaid average curve (dark line).

Figure 7 . 4 :

 74 Figure 7.4: Acidity data 7.2.1.Figure7.5: Enzyme data 7.2.2.

Figure 7 . 5 :

 75 Figure 7.4: Acidity data 7.2.1.Figure7.5: Enzyme data 7.2.2.

Figure 7 . 6 :

 76 Figure 7.6: Fishery data 7.2.4.Figure 7.7: Darwin's data 7.2.3.

Figure 7 . 7 :

 77 Figure 7.6: Fishery data 7.2.4.Figure 7.7: Darwin's data 7.2.3.

  Usage K.MixReparametrized(xobs, k, alpha0, alpha, Nsim) Arguments xobs vector of the observations or dataset k number of components in the mixture model alpha0 hyperparameter of Dirichlet prior distribution of the mixture model weights which is .5 by default alpha hyperparameter of beta prior distribution of the radial coordinate which is .5 by default Nsim number of MCMC iterations after calibration step of proposal scales Details The output of this function contains a simulated sample for each parameter of the mixture distribution, the evolution of the proposal scales and acceptance rates over the number of iterations during the calibration stage, and their final values after calibration.

  mean global vector of simulated draws from the conditional posterior of the mixture model mean sigma global vector of simulated draws from the conditional posterior of the mixture model standard deviation weights matrix of simulated draws from the conditional posterior of the mixture model weights with a number of columns equal to the number of components k angles xi matrix of simulated draws from the conditional posterior of the angular coordinates of the component standard deviations with a number of columns equal to k -1 phi vector of simulated draws from the conditional posterior of the radian coordinate angles varpi matrix of simulated draws from the conditional posterior of the angular coordinates defined for component means with a number of columns equal to k -2 accept rat vector of resulting acceptance rates of the proposal distributions without calibration step of the proposal scales optimal para vector of resulting proposal scales after optimization obtained by adaptive MCMC adapt rat list of acceptance rates of batch of 50 iterations obtained when calibrating the proposal scales by adaptive MCMC. The number of columns depends on the number of proposal distributions. adapt scale list of proposal scales calibrated by adaptive MCMC for each batch of 50 iterations with respect to the optimal acceptance rate. The number of columns depends on the number of proposal distribution scales. component means matrix of MCMC samples of the component means of the mixture model with a number of columns equal to k component sigmas matrix of MCMC samples of the component standard deviations of the mixture model with a number of columns equal to k Note: The number of the iterations in this algorithm is automatically determined depending on the convergence of the generated samples for the means and standard deviations of the components. estimate=K.MixReparametrized(xobs, k=2, alpha0=.5, alpha=.5, Nsim=10000) 8.3 Plot.MixReparametrized function Plot.MixReparametrized plot of the MCMC output produced by K.MixReparametrized Description This is a generic function for a graphical rendering of theMCMCsamples produced by K.MixReparametrized function. The function draws boxplots for unimodal variables and for multimodal arguments after clustering them by applying a k-means algorithm. It also plots line charts for other variables.. Usage plot.MixReparametrized(xobs, estimate) Arguments xobs vector of the observations estimate output of the K. MixReparametrized function Details Boxplots are produced using the boxplot.default method. Value The output of this function consists of Description This is a generic function that summarizes the MCMC samples produced by K.MixReparametrized. The function invokes several estimation methods which choice depends on the unimodality or multimodality of the argument. Usage SM.MixReparametrized(xobs, estimate) Arguments xobs vector of the observations estimate output of K.MixReparametrized

  and median of simulated draws from the conditional posterior of the mixture model mean Sd vector of mean and median of simulated draws from the conditional posterior of the mixture model standard deviation Phi vector of mean and median of simulated draws from the conditional posterior of the radial coordinate Angles. 1. vector of means of the angular coordinates used for the component means in the mixture distribution Angles. 2. vector of means of the angular coordinates used for the component standard deviations in the mixture distribution weight.i vector of mean and median of simulated draws from the conditional posterior of the component weights of the mixture distribution; i = 1, . . . , k mean.i vector of mean and median of simulated draws from the conditional posterior of the component means of the mixture distribution; i = 1, . . . , k sd.i vector of mean and median of simulated draws from the conditional posterior of the component standard deviations of the mixture distribution; i = 1, . .. , kAcc rat vector of final acceptance rate of the proposal distributions of the algorithm with no calibration stage for the proposal scales Opt scale vector of optimal proposal scales obtained by calibration stage Note: For multimodal outputs such as the mixture model weights, component means, and component variances, for each MCMC draw, first the labels of the weights p i , i = 1, . . . , k and corresponding component means and standard deviations are permuted in such a way that p 1 ≤ . . . ≤ p k . Then the component means and standard deviations are jointly partitioned into k clusters by applying a standard k-means algorithm with k clusters to a sample of size T k (where T is the number of iterations), followingFruhwirth-Schnatter (2006) method. For each group, cluster centers are considered as parameter estimates.

2. 2

 2 If d 1 = 50 compute ar j � If ar j < ar opt do log(ε θ )min(0.01, 1 / √ t), d 1 = 0 and j = j + 1; � If ar j > ar opt do log(ε θ ) + min(0.01, 1 / √ t), d 1 = 0 and j = j + 1; � Otherwise d 1 = d 1 + 1.

2. 3

 3 If t ≥ T t and d 2 = 500 compute AR l � Compute d l (AR l , ar opt ) = |AR lar opt |; � Stop the algorithm if d l (AR l , ar opt ) < δ; � Otherwise d 2 = d 2 + 1 and l = l + 1. Go to line 2.1.

Figure 9 . 1 :

 91 Figure 9.1: Pseudo-code representation of the Adaptive Metropolis-within-Gibbs algorithm used in K.MixReparametrized function of Ultimixt. Note that for each t, line [2.1] is done according to a Metropolis-within-Gibbs algorithm step described in Figure 6.1 of Chapter 6.

Figure 9 . 2 :

 92 Figure 9.2: Mixture of 5 normal distributions 9.2: (Left) Histogram of a sample of size 50; (Right) Representation of 500 MCMC iterations as mixture distributions with the overlaid average curve (dark line) over 10 4 MCMC iterations obtained by applying Plot.MixReparametrized function.

  Finally, calling two functions SM.MixReparametrized and SM.MAP.MixReparametrized leads to the following results for the component specific parameters:########################### $'Angles. 1.' [1] 0.8438176 0.9496362 3.8469520 ############################## $'Angles. 2.' [1] 1.3154771 0.3632422 1.0788443 0.8003213 ############################## Component means,standard deviations and weights: weight weight.1 weight.2 weight.3 weight.4 Median 0.23842 0.3437 0.2859 0.04079 0.09119 mean 0.23283 0.3473 0.2860 0.04080 0.09307

Figure 9 . 3 :

 93 Figure 9.3: Mixture of 5 normal distributions 9.2: Posterior distributions of the global and component wise parameters of the mixture model by applying Plot.MixReparametrized function on the MCMC output.

Figure 9 . 4 :

 94 Figure 9.4: Mixture of 5 normal distributions 9.2: Evolution of the proposal scales and related acceptance rates over the number of the batch of size 50.

#

  -02 1.6e-02 2.2e+03 5.2e+02 3.1e-02 1.6e-01The output of the last two functions helps us to easily make inference for the parameter of the mixture model. If we compare the results of the functions SM.MixReparametrized and SM.MAP.MixReparametrized, we can see that both methods of k-means clustering algorithm and the one based on MAP estimate result in the same Bayesian inference about the mixture parameters. Résumé Dans le cas particulier de l'approche bayésienne, la solution à la comparaison de modèles est le facteur de Bayes. Le facteur de Bayes est très peu utilisé. La première partie de cette thèse traite d'un examen général sur les lois a priori non informatives, et montre la stabilité globale des distributions postérieures. Dans la deuxième partie de la thèse, nous considérons un nouveau paradigme pour les tests bayésiens d'hypothèses en définissant une alternative à la construction traditionnelle de probabilités a posteriori qu'une hypothèse est vraie. Cette méthode se fonde sur l'examen des modèles en compétition en tant que composants d'un modèle de mélange. Dans la dernière partie de la thèse, nous sommes intéressés à la construction d'une analyse bayésienne de référence pour mélanges de gaussiennes par la création d'une nouvelle paramétrisation centrée sur la moyenne et la variance de ces modèles, ce qui nous permet de développer une loi a priori non-informative pour les mélanges avec un nombre arbitraire de composants. Cette partie de la thèse est suivie par une package R nommée Ultimixt qui met en oeuvre une description de notre analyse bayésienne générique de mélanges de gaussiennes.

Table 2 .

 2 concludes with recommendations for prior checks. These recommendations are mostly sensible if mainly expressing the fact that some prior information is almost always 5: Posterior means of H(θ) for the priors shown in Figure2.6 and Jeffreys' prior on θ for sample sizes 50, 100, 250, 1000, 10, 000.

	Chapter 2. Reflecting about Selecting Noninformative Priors
	Sample size	50	100	250	1000 10,000
	Dirichlet prior when γ i = 0.1
	Posterior mean 0.308 0.336 0.403 0.383 0.395
	Dirichlet prior when γ i = 0.25
	Posterior mean 0.317 0.438 0.417 0.387 0.396
	Dirichlet prior when γ i = 0.5
	Posterior mean 0.378 0.368 0.423 0.387 0.397
		Dirichlet prior when γ i = 1
	Posterior mean 0.454 0.425 0.441 0.390 0.396
		Jeffreys' prior: γ i = 0.125
	Posterior mean 0.413 0.411 0.406 0.390 0.396
	Posterior s.d	0.058 0.057 0.037 0.018 0.006

Table 4 . 1 :

 41 Dataset Pima.tr: Posterior medians of the mixture model parameters.

	352 -4.06 .103	-2.51 .064
	.2 .427 -4.03 .103	-2.49 .064
	.3 .440 -4.02 .102	-2.49 .063
	.4 .456 -4.01 .102	-2.48 .063
	.5 .449 -4.05 .103	-2.51 .064

Table 4

 4

		998 4.940 1.480 2.460 .640 .003 7.617 1.777 3.547 .786
	.2	.972 4.935 1.490 2.459 .650 .039 7.606 1.778 3.542 .787
	.3	.918 4.942 1.484 2.463 .646 .088 7.624 1.781 3.550 .788
	.4	.872 4.945 1.485 2.464 .646 .141 7.616 1.791 3.547 .792
	.5	.836 4.947 1.489 2.465 .648 .186 7.596 1.782 3.537 .788

.2: Simulated dataset: Posterior medians of the mixture model parameters.

Table 4

 4

	True model		E(α j |y)			med(α j |y)		P(M j |y)
	a 0	.1	.3	.5	.1	.3	.5
	M 1	.952 .843	.791 1	1	.936 .465
	M 2	.983 .962	.786 .989 .994	.915 .411
	M 3	.976 .973	.821 1	1	.921 .494
	M 4	.991 .867	.902 1	.987	.934 .503
	M 5	.940 .952	.896 .978 .975	.909 .591
	M 6	.974 .939	.898 1	1	.940 .617
	M 7	.973 .899	.906 1	1	1	.888
	M 8	.991 .918	.924 1	1	1	.938
	M 9	.953 .940	.878 1	.993	.956 .505
	M 10	.951 .967	.849 .988 .988	.947 .663
	M 11	.958 .951	.820 1	.989	.971 .099
	M 12	.969 .964	.951 .995 .967	.943 .196
	M 13	.919 .951	.872 1	.962	.926 .547
	M 14	.952 .964	.890 .998 .981	.911 .126
	M 15	.991 .991	.955 1	.994	.908 .164

.4: Example 4.3.6: Comparison between posterior probabilities of the true linear models, posterior means and medians of the mixture model weights, averaged over 50 replicas of samples of size 500.

Table 5 . 1 :

 51 Poisson versus Geometric 5.2.1: Mean absolute error.

	estimation model

  θ, 1) versus N (θ, 2)

		Chapter 5. Supplementary material: Testing hypotheses as a mixture	81
																		estimation model
			D.1 D.1 D.2 D.2	D.3 D.4 D.3 D.4	D.5 D.6 D.5 D.6 D.7 D.8 D.7 D.8		D.9 D.10 D.11 D.9 D.10 D.11
	N N	10 10	30	30	40	40	50	50	70 70 90	90	110 310 110 310		510 810 510 810	1000 1000
	θ	θ	-10.6 -6.4 -9.3 -6.2	-3.9 -2.1 -3.9 -1.6	2	2.2 6.1	6	10 10.1 -3.9	-4.1		0	0	6	6	10	10
		sd. sd.	0.33 0.44 0.26 0.18	0.16 0.15 0.24 0.2	0.12 0.11 0.095 0.06 0.17 0.15 0.13 0.08 0.06 0.05 0.04 0.04 0.04	0.03
	0 = 0.1 θ md. -9.31 -6.18 -3.87 -1.62 1.97 6.12 10.12 -3.87 0.02 6.01 10.01 θ md. -10.65 -6.38 -3.92 -2.09 2.17 6.01 10.02 -4.07 0.03 5.97 9.97 gd. 1 1 1 1.02 1.01 1 1 1 1 1 1.01 sd. 0.26 0.09 0.07 0.14 0.09 0.07 0.04 0.02 0.02 0.01 0.02 gd. 1 1 1.01 1 1 1 1 1 1 1 1 0 = 0.1 sd. 0.13 0.11 0.08 0.12 0.06 0.16 0.09 0.09 0.02 0.02 0.03
	a a α md. 0.00 0.00 α md. 0.99 0.99	0.99 0.99 0.00 0.00 0.00 0.00 0.00 0.00 0.99 0.99 0.99 1	1 6e-5 4e-5 1	2e-5	0.99
		gd. gd.	1	1	1	1.01	1	1	1	1	1	1.001 1 1	1 1.01 1	1.01	1	1	1 1.01	1	1
		sd. sd.	0.33 0.44 0.24 0.19	0.16 0.15 0.22 0.21 0.17 0.14 0.13 0.08 0.12 0.11 0.09 0.06 0.06 0.05 0.05 0.04 0.04	0.03
	0 = 0.2 θ md. -9.33 -6.35 -3.88 -1.62 1.97 6.12 10.13 -3.88 0.01 6.01 10.01 θ md. -10.64 -6.38 -3.93 -2.1 2.18 6.02 10.01 -4.06 -0.02 5.96 10.01 gd. 1.002 1 1 1 1.02 1 1 1 1 1 1 sd. 0.27 0.11 0.11 0.17 0.11 0.08 0.06 0.03 0.02 0.014 0.02 gd. 1.01 1 1 1.02 1 1 1.01 1 1 1 1 0 = 0.2 sd. 0.16 0.14 0.11 0.14 0.08 0.18 0.12 0.10 0.03 0.02 0.03
	a a α md. 0.01 0.01 α md. 0.98 0.99	0.99 0.98 0.01 0.02 0.01 0.05 0.01 0.04 0.99 0.99 0.99 0.99 1e-3 1e-3 0.99 0.99 0.00	0.99
		gd. gd.	1.01 1	1	1	1	1	1	1	1	1	1	1	1	1.02 1 1		1	1	1.01 1.01	1	1
		sd. sd.	0.34 0.45 0.24 0.19	0.16 0.15 0.24 0.18 0.17 0.14 0.14 0.08 0.12 0.11 0.098 0.06 0.06 0.05 0.04 0.04 0.05	0.03
	0 = 0.3 θ md. -9.36 -6.18 -3.88 -1.62 1.98 6.11 10.12 -3.87 0.02 6.02 10.01 θ md. -10.54 -6.41 -3.91 -2.01 2.13 6.00 9.98 -4.09 0.01 6.01 9.97 gd. 1 1 1 1 1.01 1 1 1 1 1 1 sd. 0.28 0.13 0.12 0.17 0.12 0.10 0.07 0.03 0.03 0.02 0.02 gd. 1 1 1.01 1 1.01 1 1.01 1 1 1 1 0 = 0.3 sd. 0.18 0.15 0.11 0.16 0.1 0.19 0.13 0.10 0.04 0.03 0.04
	a a α md. 0.04 0.04 α md. 0.94 0.98	0.98 0.95 0.02 0.05 0.02 0.12 0.05 0.07 0.97 0.98 0.98 0.99 0.01 0.01 0.99 0.99 0.01	0.99
		gd. gd.	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1.02 1.01 1 1	1	1	1.01
		sd. sd.	0.34 0.45 0.24 0.19	0.16 0.15 0.22 0.2	0.12 0.11 0.10 0.06 0.16 0.14 0.14 0.08 0.06 0.05 0.05 0.03 0.04	0.03
	0 = 0.4 θ md. -9.37 -6.18 -3.88 -1.62 1.97 6.12 10.12 -3.87 0.02 6.00 10.01 θ md. -10.63 -6.37 -3.90 -2.13 2.17 6.01 10.01 -4.1 0.03 5.97 9.98 gd. 1 1 1 1 1 1.01 1 1 1 1 1 sd. 0.27 0.14 0.12 0.18 0.13 0.11 0.08 0.03 0.03 0.02 0.03 gd. 1.01 1 1 1 1 1 1.01 1 1.01 1 1.01 0 = 0.4 sd. 0.18 0.17 0.13 0.17 0.11 0.19 0.13 0.11 0.04 0.03 0.02
	a a α md. 0.06 0.06 α md. 0.92 0.96	0.97 0.92 0.05 0.08 0.04 0.17 0.07 0.09 0.95 0.96 0.98 0.99 0.01 0.01 0.99 0.99 0.02	0.99
		gd. gd.	1	1	1	1	1	1	1	1	1	1	1	1	1 1.01 1	1		1	1.02 1.01 1	1	1
		sd. sd.	0.34 0.45 0.24 0.19	0.16 0.15 0.21 0.19 0.16 0.14 0.13 0.08 0.12 0.11 0.10 0.06 0.06 0.05 0.04 0.04 0.04	0.03
	0 = 0.5 θ md. -9.35 -6.18 -3.87 -1.62 1.98 6.12 10.12 -3.79 -0.02 6.01 θ md. -10.65 -6.42 -3.89 -2.11 2.18 5.98 10 -4.07 -0.03 6 gd. 1 1.01 1 1 1 1 1 1 1 1 sd. 0.27 0.15 0.14 0.19 0.14 0.12 0.08 0.04 0.03 0.02 9.98 10.03 1 0.03 gd. 1.01 1 1 1 1 1 1.01 1 1 1 1 0 = 0.5 sd. 0.19 0.17 0.14 0.18 0.11 0.19 0.14 0.11 0.04 0.03 0.04
	a a α md. 0.09 0.08 α md. 0.88 0.95	0.96 0.89 0.07 0.11 0.05 0.18 0.09 0.12 0.93 0.95 0.97 0.99 0.02 0.01 0.99 0.99 0.02	0.98
		gd. gd.		1		1		1		1		1		1		1		1.02		1	1.01	1.01
		Table 5.2: N (θ, 1) versus N (θ, 2) 5.3: Posterior summaries; Datasets (D.1, ..., D.11)
		are simulated from N (θ, 1) and each point estimator is based on 10, 000 iterations of the
		Metropolis-Hastings algorithm.

Table (

 (

		a)
	Data: D.1	D.2 D.3	D.4 D.5	D.6 D.7	D.8
	N	5	15	25	35	55	65	75	95
	μ Mean -0.17 0.16 0.2 Sd. 0.83 0.82 0.9	0.09 -0.03 0.84 0.86	0.02 0.00 0.89 0.86	0.02 0.85
	Data: D.9 D.10 D.11 D.12 D.13 D.14 D.15 D.16
	N	250	450	650	750	850	950	1000	5e+4
	μ Mean -0.02 -0.02 -0.01 -0.02 -0.08 Sd. 0.86 0.81 0.83 0.85 0.88	0.05 -0.10 -0.03 0.89 0.80 0.90
	Table (b)
	Data: D.1	D.2 D.3	D.4 D.5	D.6 D.7	D.8
	N	5	15	25	35	45	55	65	250
	µ	-9	-8.3 -7.8	-6.8 -5.8	-5	-4.1	2
	μ Mean -8.93 -8.38 -7.75 -6.89 -5.75 -4.96 -4.06 Sd. 0.41 0.24 0.19 0.17 0.15 0.13 0.12	1.94 0.06
	Data: D.9 D.10 D.11 D.12 D.13 D.14 D.15 D.16
	N	350	450	550	650	750	850	950	1000
	µ	3	4	4.9	6	6.9	8	8.9	10
	μ Mean 3.01 Sd. 0.05	4.00 4.98 0.05 0.04	6.00 6.93 0.04 0.04	7.98 8.98 0.03 0.03	10.04 0.03
	Table 5.4:

N (0, 1) versus N (µ, 1) 5.4: Observation information and posterior summaries; Table (a) Datasets (D.1, ..., D.16) are simulated from N (0, 1); Table (b) Datasets (D.1, ..., D.16) are simulated from N (µ, 1). a 0 is supposed to be .1 and each point estimator is based on 10, 000 MCMC iterations.

  ). mixtures Proposal 1 ar µ ar σ ar p ar ϕ,η

		ε µ	ε p	ε
	0.40 0.47 0.45 0.24 0.56 77.06 99.94
	Proposal 2 ar µ ar σ ar p ar ϕ,η	ε µ	ε p	ε ϑ
	0.38 0.46 0.45 0.27 0.55 0.29	0.35
	Table 6.1: Example 6.4.1: Acceptance rate (ar) and corresponding proposal scale (ε)

  ], when we impose k = 6 components. The MCMC sample is again summarized by k-means clustering and MAP estimates, as presented in Section 6.3.2. The resulting means, medians and 95% credible intervals of the parameters of the mixture components

										mixtures
			6 components, k = 6		4 components, k = 4
				k-means clustering			k-means clustering
		p 1	p 2	p 3	p 4	p 5	p 6	p 1	p 2	p 3	p 4
	Median	0.01	0.08	0.13	0.43	0.05	0.24	0.56	0.27	0.06	0.10
	Mean	0.02	0.06	0.14	0.46	0.05	0.24	0.58	0.25	0.06	0.11
		µ 1	µ 2	µ 3	µ 4	µ 5	µ 6	µ 1	µ 2	µ 3	µ 4
	Median	25.95	9.72	22.06	19.83	32.71	22.87	20.19	21.52 32.79	9.72
	Mean	61.59	9.725	22.09	19.84	32.70	22.93	20.27	21.48 33.29	9.73
		σ 1	σ 2	σ 3	σ 4	σ 5	σ 6	σ 1	σ 2	σ 3	σ 4
	Median	4.53	4.91	1.91	0.52	2.86	0.65	0.52	1.62	3.00	1.05
	Mean	32.23	4.61	2.41	0.58	4.23	1.10	0.57	2.08	3.66	3.44
				MAP estimate				MAP estimate
		p 1	p 2	p 3	p 4	p 5	p 6	p 1	p 2	p 3	p 4
	Median	0.04	0.09	0.13	0.37	0.10	0.15	0.32	0.46	0.08	0.08
	Mean	0.05	0.09	0.10	0.39	0.14	0.22	0.34	0.43	0.13	0.09
	2.5%	< 10 -5	< 0.01	< 0.01	< 0.01	< 10 -3	< 0.01	0.04	0.02	0.01	0.04
	97.5%	0.2.1	0.13	0.69	0.39	0.56	0.68	0.87	0.82	0.51	0.15
		µ 1	µ 2	µ 3	µ 4	µ 5	µ 6	µ 1	µ 2	µ 3	µ 4
	Median	30.96	9.70	21.75	19.73	20.61	23.12	19.84	22.17 28.23	9.71
	Mean	67.26	8.18	21.58	18.73	20.84	24.33	19.83	22.34 29.03	9.50
	2.5%	22.87	-9.28	19.60	9.68	12.83	21.29	17.59	20.14 22.27	9.17
	97.5%	606.16	10.21	23.44	20.47	25.69	33.07	21.47	26.87 36.20	10.21
		σ 1	σ 2	σ 3	σ 4	σ 5	σ 6	σ 1	σ 2	σ 3	σ 4
	Median	4.82	0.54	1.76	0.60	3.41	1.73	0.69	2.22	3.22	0.53
	Mean	20.53	2.05	2.06	0.73	15.59	2.34	0.96	3.23	4.15	0.91
	2.5%	0.79	0.30	0.31	0.19	0.41	0.17	0.29	0.87	0.68	0.29
	97.5%	198.23	17.28	7.63	2.13	35.95	7.62	2.44	9.62	10.57	1.34
			Gibbs sampler (bayesm)		Gibbs sampler (bayesm)
		p 1	p 2	p 3	p 4	p 5	p 6	p 1	p 2	p 3	p 4
		0.17	0.09	0.14	0.23	0.19	0.19	0.33	0.31	0.18	0.18
		µ 1	µ 2	µ 3	µ 4	µ 5	µ 6	µ 1	µ 2	µ 3	µ 4
		19.59	21.97	20.83	22.21	19.93	22.00	20.53	21.05 21.75	20.90
		σ 1	σ 2	σ 3	σ 4	σ 5	σ 6	σ 1	σ 2	σ 3	σ 4
		0.35	0.23	0.22	0.24	0.26	0.31	0.22	0.19	0.21	0.27
			EM estimate (mixtools)		EM estimate (mixtools)
		p 1	p 2	p 3	p 4	p 5	p 6	p 1	p 2	p 3	p 4
		0.04	0.08	0.17	0.41	0.09	0.20	0.52	0.33	0.04	0.11
		µ 1	µ 2	µ 3	µ 4	µ 5	µ 6	µ 1	µ 2	µ 3	µ 4
		24.27	9.71	22.33	19.88	33.04	24.26	19.72	22.72 33.04	10.14
		σ 1	σ 2	σ 3	σ 4	σ 5	σ 6	σ 1	σ 2	σ 3	σ 4
		0.08	0.42	0.44	.70	0.19	8.33	0.62	1.77	0.92	2.73
	Table 6.3: Galaxy dataset: Estimates of the parameters of a mixture of 6 and 4
	components.

Table 7 . 1 :

 71 Mixture of two normal distributions 7.1: Point estimates, proposal scales and acceptance rates.

  Table 7.2: Data analyses 7.2: Point estimators of the parameters of a mixture of two components. Each estimate is obtained based on 10 4 MCMC iterations. Med indicates the estimate based on the median of draws and two values behind Obs. are the mean and standard deviation of the dataset x =

	Med Mean Med Mean Med Mean Med Mean	Angular & component-wise parameters Acidity data Enzyme data Darwin's data ξ ξ ξ 0.71 1.39 0.13 0.75 1.4 0.12 p 1 p 2 p 3 p 1 p 2 p 3 0.40 0.60 0.41 0.59 0.85 0.15 0.41 0.59 0.40 0.60 0.84 0.16 µ 1 µ 2 µ 1 µ 2 µ 1 µ 2 6.24 4.32 1.25 0.19 21.68 12.0 6.23 4.33 1.24 0.19 21.63 12.2 σ 1 σ 2 σ 1 σ 2 σ 1 σ 2 0.53 0.38 0.52 0.08 0.33 1.5 0.54 0.37 0.53 0.08 0.47 1.4	Acidity Enzyme Darwin	Obs. Med Mean Obs. Med Mean Obs. Med Mean	Global parameters 5.1 1.04 µ σ ϕ 5.1 1.04 0.90 5.1 1.04 0.89 0.62 0.62 µ σ ϕ 0.62 0.63 0.84 0.62 0.62 0.83 20.2 3.6 µ σ ϕ 20.2 3.7 0.92 20.2 3.6 0.92
		� n j=1	xj /n; s = (	� n j=1	(xj -x) /(n-1)) 1 /2 , respectively.

Table 7 .

 7 3: Fishery data 7.2: Point estimators of the parameters of a mixture of two components. Each estimate is obtained based on 10 4 MCMC iterations. Med indicates the estimate based on the median of draws and two value behind Obs. are the mean and standard deviation of the dataset x =

	� n j=1	xj /n; s = (	� n

j=1

(xj -x) /(n-1)) 1 /2 , respectively.

Using a Bayes factor to test for the number of components in the mixture (6.1) as in[START_REF] Richardson | [END_REF]] would be possible. However, the outcome would fail to answer the original question of selecting between both (or more) models.

For this example, the conditional evidence π(x|ζ) can easily be derived in closed form, which means that a random walk on the allocation space {1, 2} n could be implemented. We did not follow that direction, as it seemed unlikely such a random walk would have been more efficient than a Metropolis-Hastings algorithm on the parameter space only.

are simulated from N (θ, 2) and each point estimator is based on 10, 000 iterations of the Metropolis-Hastings algorithm.

Acknowledgments

Chapter 6. Non-informative reparameterisations for location-scale mixtures such models, due to the invariance of the sum by permutation of the indices. This issue relates to the equally well-known "label switching" phenomenon in the Bayesian approach to the model, which pertains both to inference and to simulation of the corresponding posterior [Celeux 2000, Stephens 2000, Frühwirth-Schnatter 2001, Frühwirth-Schnatter 2004[START_REF] Jasra | [END_REF]. From this Bayesian viewpoint, the choice of the prior distribution on the component parameters is quite open, the only constraint being that the corresponding posterior is proper [Diebolt 1994, Frühwirth-Schnatter 2004]. [Diebolt 1994] and [Wasserman 1999] discussed the alternative approach of imposing proper posteriors on improper priors by banning almost empty components from the likelihood function. While consistent, this approach induces dependence between the observations, higher computational costs and is not handling overfitting very well. It has therefore seen little following.

The prior distribution on the weights p i is equally open for choice, but a standard version is a Dirichlet distribution with common hyperparameter a, Dir(a, . . . , a). Recently, [START_REF] Rousseau | [END_REF]] demonstrated that the choice of this hyperparameter a relates to the inference on the total number of components, namely that a small enough value of a manages to handle over-fitted mixtures in a convergent manner. In a Bayesian non-parametric modeling, [Griffin 2010] showed that the prior on the weights may have a higher impact when inferring about the number of components, relative to the prior on the component-specific parameters. As indicated above, the prior distribution on the θ i 's has received less attention and conjugate choices are most standard, since they facilitate simulation via Gibbs samplers [Diebolt 1990, Escobar 1995[START_REF] Richardson | [END_REF] if not estimation, since posterior moments remain unavailable in closed form. In addition, [START_REF] Richardson | [END_REF]] among others proposed data-based priors that derive some hyperparameters as functions of the data, towards an automatic scaling of such priors. An R package, bayesm [START_REF] Rossi | [END_REF]] incorporates some of those ideas. In the case when θ i = (µ i , σ i ) is a location-scale parameter, [Mengersen 1996] proposed a reparameterisation of (6.1) that express each component as a local perturbation of the previous one, namely (i > 1)

with µ 1 and σ 1 being the reference values. Based on this reparameterisation, [Robert 1998] established that a particular improper prior on (µ 1 , σ 1 ) still leads to a proper prior. We propose here to modify further this reparameterisation towards using the global mean and global variance of the mixture distribution as reference location and scale, respectively. This modification has foundational consequences in terms of using improper and non-informative priors over mixtures, in sharp contrast with the existing literature (see, e.g. [Diebolt 1993, Diebolt 1994, O'Hagan 1994, Wasserman 1999]). Bayesian computing for mixtures covers a wide variety of proposals, starting with the introduction of the Gibbs sampler [Diebolt 1990[START_REF] Gelman | [END_REF], Escobar 1995], some concerned with approximations [Roeder 1990, Wasserman 1999] and MCMC features [START_REF] Richardson | [END_REF], Celeux 2000, Casella 2002], and others with asymptotic mixtures Metropolis-within-Gibbs algorithm for reparameterised mixture model

2 For t = 1, . . . , T , the update of (µ (t) , σ (t) , p (t) , ϕ (t) , ξ 1) , p (t-1) , ϕ (t-1) , ξ (t-1) , � (t-1) ).

2.2 Generate a proposal log(σ) � ∼ N (log(σ (t-1) ), ε σ ) and update σ (t) against π(•|x, µ (t) , p (t-1) , ϕ (t-1) , ξ (t-1) , � (t-1) ). 1) , ξ (t) , � (t) ). ε p + 1), and update p (t) against π(•|x, µ (t) , σ (t) , ϕ (t) , ξ (t) , � (t) ).

Generate proposals ξ

2.6

Generate proposals ξ

Figure 6.1: Pseudo-code representation of the Metropolis-within-Gibbs algorithm used in this paper for the reparameterisation (ii) based on two sets of spherical coordinates. For simplicity's sake, we denote p (t) = (p

MCMC implications 113

between the current posterior sample and the (or a) maximum posterior probability (MAP) estimate.

Let us denote by h the map from our reparameterisation to the standard parameterisation of (6.1), i.e.,

with its inverse h -1 available as well. We also denote by S k the set of permutations of {1, . . . , k}. Then, given an MCMC sample {p (t) , θ (t) } T t=1 , the above relabelling technique procedure follows;

1. Reparameterise the MCMC sample {p (t) , θ (t) } T t=1 into component-wise means and standard deviations via the function h, resulting in {µ

2. Find the MAP estimate by computing the posterior values of the sample; denote the solution as

3. Reorder (µ

k , p (t) ) as

where δ j = arg min δ∈S k �δ(µ

The resulting permutation is then denoted λ (t) ∈ S k . Label switching occurrences in an MCMC sequence can be monitored via the changes in the sequence λ (1) , . . . , λ (T ) . If the chain fails to switch modes, the sequence is likely to remain at the same permutation. On the opposite, if a chain moves between some of the k! symmetric posterior modes, the λ (t) 's are expected to vary.

We proceed here by a simulation studies section and all algorithms used in this section are publicly available within the R package Ultimixt [Kamary 2015]. The package Ultimixt contains functions that implement adaptive determination of optimal scales and convergence monitoring based on [Gelman 1992] criterion. In addition, Ultimixt includes functions that summarize the simulations and compute point estimates of each parameter, such as posterior mean and median. It also produces an estimated mixture density in numerical and graphical formats. The output further includes graphical representations of the generated parameter samples. For the potentially unimodal parameters µ, σ and ϕ, averaging and calculating the median over the generated chains directly returns valid point estimators, as those parameters are not subjected to label switching. For the other parameters (component weights, means and variances), since label switching is a possible issue, we need to postprocess the MCMC draws as discussed earlier, by first relabelling these simulations. We then derive point estimates by clustering over the parameter space, using k-mean clustering [Hastie 2001].

Chapter 6. Non-informative reparameterisations for location-scale mixtures

Simulation studies

In this section, we examine the performances of the above Metropolis-within-Gibbs algorithm, when applied to both reparameterisations defined above. We also consider the special case k = 2 in Section 6.4.1. All simulations were conducted using the package Ultimixt [Kamary 2015].

The case k = 2

In this specific case, we do not have to simulate any angle. Two straightforward proposals are compared over simulation experiments. One is based on Beta and Dirichlet proposals:

(this will be called Proposal 1) and another one is based on Gaussian random walks:

and θ * = 1 + exp(ϑ * 1 ) + exp(ϑ * 2 ) (which will be called Proposal 2). The global parameters are proposed using Normal and Inverse-Gamma proposals

where x and σ2 are sample mean and variance respectively. We present below some analyses and also explain how MCMC methods can be used to fit the reparameterised mixture distribution.

Example 6.4.1 In this experiment, a dataset of size 50 is simulated from the mixture 0.65N (-8, 2) + 0.35N (-0.5, 1), which implies that while the true value of (ϕ, η 1 , η 2 ) is (0.91, 0.16, 0.38). We also run the Metropolis-within-Gibbs algorithm based on Proposal 2 using the same simulated dataset for comparison purposes. As shown in Figure 6.3, the outputs for both proposals are quite similar but Proposal 1 produces more symmetric chains on p, ϕ, η 1 , η 2 , thus suggesting higher mixing abilities.

Chapter 6. Non-informative reparameterisations for location-scale mixtures are displayed in Table 6.3. Unsurprisingly, global mean and standard deviation are quite similar to the empirical estimates. Table 6.3 also displays estimates based on the Gibbs sampler of bayesm [START_REF] Rossi | [END_REF]] and on the EM algorithms of mixtools [Benaglia 2009], with our approach being produced by Ultimixt [Kamary 2015].

Obtaining very close estimations for two component means µ i , as µ 1 = 19.59 and µ 5 = 19.93, and µ 2 = 21.97 and µ 6 = 22 and µ 4 = 22.21 for bayesm, and µ 1 = 24.27 and µ 6 = 24.26 for mixtools, signals that overfitting occurs: there are more components than supported by the data. With our analysis, overfitting is handled in a different way: the mean of one or more component weights is close to zero. For instance, we obtained estimates of p 1 very close to zero, inducing estimates for µ 1 and σ 1 of 61.59 and 32.23 for µ 1 and σ 1 (obtained by k-means clustering) and of 67.26 and 20.53 (using MAP estimates), as shown in Table 6.3. If we examine the MCMC sequences in detail, the minimum simulated value for the first component weight and the corresponding first component mean and standard deviation are 1.045 10 -6 , 449.25 and 284.34, respectively. Such extreme values are produced because of the extremely small weight. However, such large values have no impact on the resulting estimate of the mixture itself. This is clearly exhibited in Figure 6.8 for the Galaxy dataset, which shows that extreme values have no effect on the predictive density plots due to the small weights. Using our modeling, the resulting density estimate is remarkably smooth when considering that the number of observations is 82 and a number of components equal to 6.

If we repeat running the algorithm on the Galaxy dataset for 50, 000 iterations and a smaller number of components, for instance k = 4, summary and model fit statistics are provided in Table 6.3. In this case, extreme values do not occur and the predictive density plots show that a four component model fits the data equally well as displayed in Figure 6.9. The posterior estimates of the component parameters computed by three methods (k-means clustering, MAP, and EM estimates) are almost similar, while the Gibbs sampler results from bayesm yield two very close estimates of component means, µ 2 = 21.05 and µ 4 = 20.90 in this case.

The common priors for the standard parameters are

where IW(ν, 3) is the Inverse-Wishart distribution with the scale parameter of 3 and the degrees of freedom of ν. Unknown hyperparameters μ, σ R , α 0 and ν are given by bayesm from the empirical estimation of data and, the comparison of the proposed priors and the prior obtained from bayesm are graphically presented in Figure 6.10.

Parallel tempering algorithm for reparameterised mixture when k = 2:

1 Generate initial values θ (0) ; 2 Compute ε θ by implementing the adaptive Metropolis-within-Gibbs algorithm 9.1.

4 For t = 1, . . . , T , the update of θ (t) follows

h with probability

6 Otherwise compute β = (β, 2 -N β ), N β = N β + 1 and go to line 4. than and less than 0.5 to the total number of iterations, respectively, the degree of symmetry of draws of p to 0.5 is computed by

The output of the algorithm 7.8 consists of the chain of samples for θ, the number of replicas and the temperature set, 1 /β. Both examples 6.4.2 and 6.5.2 illustrate the good performance of this algorithm in terms of improving the mixing of the MCMC chains.

Chapter 8

Ultimixt package

Ultimixt

Type Package Title Bayesian Analysis of a Non-Informative Parametrization for Gaussian Mixture Distributions Version 1.0 Date 2015-12-10 Author Kaniav Kamary, Kate Lee Maintainer Kaniav Kamary <kamary@ceremade.dauphine.fr> Depends coda, gtools, graphics, grDevices, stats Description A generic reference Bayesian analysis of unidimensional mixtures of Gaussian distributions obtained by a location-scale parameterisation of the model is implemented. Included functions can be applied to produce a Bayesian analysis of Gaussian mixtures with an arbitrary number of components, with no need to define the prior distribution. License GPL (>=2.0) R topics documented: The label of components of θ i (before the above transform) is defined by

The number of label switching occurrences is defined by the number of changes in τ * .

Usage

Supplementary material: Ultimixt package

Mixture models as a popular method for modeling unobserved heterogeneity, find application in a very wide number of applied fields. As mentioned before, they can be used to model a statistical population with subpopulations with the densities on the subpopulations as the mixture components and the proportions of each subpopulation in the overall population as the weights. For maximum likelihood estimation of the mixture model parameters, the Expectation-Maximization (EM) algorithm is most frequently used by the classical community which is provided in R by packages such as mclust [Fraley 2002], flexmix [Gruen 2008] and mixtools [Benaglia 2009].

Bayesian estimation has become feasible with the advent of Markov chain Monte Carlo (MCMC) simulation and R packages such as BayesMix [Gruen 2015] and bayesm [START_REF] Rossi | [END_REF]] have been made for estimating univariate Gaussian finite mixtures with MCMC methods. However, the model class that is implemented in these packages only allows informative structure of the prior distributions. This chapter involves the algorithms of Ultimixt package described in Chapter 8.

Description of implementation

The first function of the package, K.MixReparametrized, provides the functionality for estimating univariate Gaussian mixture models regarding the non-informative parameterization expressed in Chapter 6 with MCMC methods. Within a given model class users can modify the prior specification of the mixture distribution weights and the number of the components for developing a suitable model for the dataset. K.MixReparametrized function shares the following features:

Calibration step which consists of determining the scales of the proposal distributions by applying adaptive Metropolis-within-Gibbs algorithm;

The motivation of this step is to avoid using the method of trial and error in order to obtain proposal scales. In some special cases, the method of trial and error can be useful to make some ideas about which value could work well in the sense that the proposal distribution results in good mixing MCMC chains. However, trial and error method is time consuming and sometimes fails to attain a satisfactory value to get the convergence with the movements in the support of the target distribution. Adaptive MCMC