Nonparametric and Variable-Dimension Bayesian Mixture Models: Analysis, Comparison, and New Methods by Jeffrey W. Miller B.Sc., Georgia Institute of Technology; Atlanta, GA, 2001 M.Sc., Stanford University; Stanford, CA, 2002 M.Sc., Brown University; Providence, RI, 2010 A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in The Division of Applied Mathematics at Brown University PROVIDENCE, RHODE ISLAND May 2014 c Copyright 2014 by Jeffrey W. Miller This dissertation by Jeffrey W. Miller is accepted in its present form by The Division of Applied Mathematics as satisfying the dissertation requirement for the degree of Doctor of Philosophy. Date Matthew T. Harrison, Ph.D., Advisor Recommended to the Graduate Council Date Stuart Geman, Ph.D., Reader Date Steven N. MacEachern, Ph.D., Reader Approved by the Graduate Council Date Peter M. Weber, Dean of the Graduate School iii Vita Jeffrey W. Miller was born on June 20, 1979, in Wilmington, Delaware, and grew up in West Chester, Pennsylvania. He received a Bachelor of Science in Mechanical Engineering from the Georgia Institute of Technology in 2001, a Master of Science in Mechanical Engineering from Stanford University in 2002, and a Master of Science in Mathematics from Brown University in 2010. He served as an officer in the U.S. Air Force during 2002–2005, at the Air Force Research Lab at Tyndall Air Force Base, Florida. During 2005–2008, he worked at Draper Laboratory in Cambridge, Massachusetts. As a PhD candidate in the Division of Applied Mathematics at Brown University, Jeffrey has been advised and mentored by Matthew T. Harrison and Stuart Geman. In addition to several teaching assistantships, he has been the primary instructor for two courses, Introduction to Machine Learning (summer 2011), and Information Theory (fall 2011), and in 2012 he received the Presidential Award for Excellence in Teaching, awarded to one Brown University graduate student each year. He has also received awards for his research, including the Sigma Xi Outstanding Graduate Student Award in 2013, 1st Prize for Poster by a Young Participant at the 9th Conference on Bayesian Nonparametrics in 2013, and the IBM Thomas J. Watson Research Center Student Research Award at the New England Statistics Symposium in 2011. iv Peer-reviewed publications J. W. Miller and M. T. Harrison, Inconsistency of Pitman–Yor process mixtures for the number of components, 2014+ (Accepted pending minor revisions), arXiv:1309.0024. J. W. Miller and M. T. Harrison, A simple example of Dirichlet process mixture inconsistency for the number of components, Advances in Neural Information Processing Systems (NIPS), Vol. 26, 2013. J. W. Miller and M. T. Harrison, Exact sampling and counting for fixed-margin matrices, The Annals of Statistics, Vol. 41, No. 3, 2013, pp. 1569-1592. J. W. Miller, Reduced criteria for degree sequences, Discrete Mathematics, Vol. 313, Issue 4, 2013, pp. 550-562. v Preface and Acknowledgments In the fall of 2011, I decided to try my hand at a problem that Stu Geman had mentioned several times — in his words, a “stone in his shoe” that had been bothering him for some time. He had noticed that Dirichlet process mixtures (DPMs) were often being used to infer the number of components in a finite mixture, a practice which seemed highly questionable, due to the nature of the prior. Still, it was not obvious how the posterior would behave — would it, nonetheless, be consistent for the number of components? After finding some special cases in which inconsistency did provably occur, I worked for several months with Matt Harrison, trying to determine the limiting behavior of the posterior. Eventually, we found that, although precisely determining the limiting behavior eluded us, we could still prove inconsistency under very general conditions. Simultaneously, we started exploring the natural alternative of putting a prior on the number of components in a finite mixture, i.e., using a variable-dimension mix- ture. Naively unaware of the previous work with this model, but being familiar with the DPM, the natural thing to do was to develop an inference algorithm resembling the DPM algorithms. This turned out to be a novel development, quite different from previous inference algorithms for the model, such as reversible jump. In developing this DPM-style inference algorithm, it became clear that a variable- vi dimension mixture could possess many of the same properties as a nonparametric mixture, in particular, representation in terms of a simple restaurant process. This led naturally to the idea of pursuing such a representation for variable-dimension alternatives to other popular nonparametric models, and to my surprise, alterna- tives to the hierarchical Dirichlet process and the Indian buffet process worked out beautifully in much the same way. There are many people, without whom, this work would hardly have been possi- ble. I would like to thank my parents for their continual support in all my endeavors, and their interest in understanding and keeping up with my work. My advisor, Matt Harrison, has been as close to ideal as I could imagine — always available and respon- sive, pushing me to get to the bottom of things and comprehensively address each question, keenly understanding tricky problems and creatively suggesting solutions, and wisely helping to choose projects worth pursuing. It has been an honor and a joy to work with Stu Geman, who, in addition to providing the spark which developed into this thesis, has been an encouraging mentor in many ways; his boundless curios- ity is an inspiration. Special thanks also to Steve MacEachern for many insightful comments and suggestions which turned out to be most helpful, even though I often didn’t realize it until later! Thanks to my “lab”-mates, Dan Klein and Dahlia Nadkarni, for many enjoyable group meetings, sharing ideas and working through problems as they arise. Thanks to the regulars of the machine learning reading group, Mike Hughes, Jason Pacheco, Dae Il Kim, Thomas Wiecki, Imri Sofer, and Mark Homer, for many discussions about Bayesian nonparametrics and machine learning. To my best friends, Lauren and Arthur Sugden, thanks for countless evenings of conversations about math, science, and life in general. vii Thanks to Tamara Broderick, Dave Duvenaud, and Konstantina Palla for sug- gesting that a variable-dimension alternative to the IBP might also be interesting to pursue. The nonparametric Bayesian community has been most welcoming and enjoyable to interact with, and numerous discussions with its members have made this work immeasurably better; in particular, I would like to thank Erik Sudderth, David Dunson, Surya Tokdar, Mike West, Vinayak Rao, Ryan Adams, Annalisa Cerquetti, Peter M¨ uller, Long Nguyen, Levi Boyles, Dan Roy, Zoubin Ghahramani, Yee Whye Teh, Amar Shah, Nick Foti, and Debdeep Pati. viii Contents Vita iv Preface and Acknowledgments vi 1 Introduction 1 2 A simple example of Dirichlet process mixture inconsistency for the number of components 4 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3 Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3.1 Dirichlet process mixture model . . . . . . . . . . . . . . . . . 10 2.3.2 Specialization to the standard normal case . . . . . . . . . . . 12 2.4 Simple example of inconsistency . . . . . . . . . . . . . . . . . . . . . 13 2.5 Severe inconsistency . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3 Inconsistency of Pitman–Yor process mixtures for the number of components 19 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.1.1 A motivating example . . . . . . . . . . . . . . . . . . . . . . 20 3.1.2 Overview of results . . . . . . . . . . . . . . . . . . . . . . . . 24 3.1.3 Discussion / related work . . . . . . . . . . . . . . . . . . . . 26 3.1.4 Intuition for the result . . . . . . . . . . . . . . . . . . . . . . 28 3.1.5 Organization of this chapter . . . . . . . . . . . . . . . . . . . 29 3.2 Model distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.2.1 Gibbs partitions . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.2.2 Gibbs partition mixtures . . . . . . . . . . . . . . . . . . . . . 33 3.3 Inconsistency theorem . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.4 Application to discrete or bounded cases . . . . . . . . . . . . . . . . 39 3.5 Exponential families and conjugate priors . . . . . . . . . . . . . . . . 42 3.5.1 Exponential families . . . . . . . . . . . . . . . . . . . . . . . 42 3.5.2 Conjugate priors . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.6 Application to exponential families . . . . . . . . . . . . . . . . . . . 44 3.7 Marginal inequalities . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 ix 3.8 Marginal inequality for subsets of the data . . . . . . . . . . . . . . . 52 3.9 Miscellaneous proofs . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.10 Exponential family properties . . . . . . . . . . . . . . . . . . . . . . 58 3.11 Bounds on the Laplace approximation . . . . . . . . . . . . . . . . . 60 3.12 Capture lemma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4 The mixture of finite mixtures model 70 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.1.1 Mixture of finite mixtures (MFM) . . . . . . . . . . . . . . . . 73 4.2 Properties of MFMs . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.2.1 Partition distribution . . . . . . . . . . . . . . . . . . . . . . . 76 4.2.2 Equivalent models . . . . . . . . . . . . . . . . . . . . . . . . 78 4.2.3 Various basic properties . . . . . . . . . . . . . . . . . . . . . 80 4.2.4 The coefficients Vn (t) . . . . . . . . . . . . . . . . . . . . . . . 81 4.2.5 Self-consistent marginals . . . . . . . . . . . . . . . . . . . . . 84 4.2.6 Restaurant process / P´olya urn scheme . . . . . . . . . . . . . 85 4.2.7 Random discrete measure formulation . . . . . . . . . . . . . . 87 4.2.8 Density estimates . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.2.9 Stick-breaking representation . . . . . . . . . . . . . . . . . . 91 4.2.10 Proofs and details . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.3 Asymptotics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.3.1 Posterior asymptotics . . . . . . . . . . . . . . . . . . . . . . . 96 4.3.2 Relationship between the number of clusters and number of components . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 4.3.3 Distribution of the cluster sizes under the prior . . . . . . . . 105 4.3.4 Asymptotics of Vn (t) . . . . . . . . . . . . . . . . . . . . . . . 107 4.4 Inference algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 4.4.1 Inference with conjugate priors . . . . . . . . . . . . . . . . . 110 4.4.2 Inference with non-conjugate priors . . . . . . . . . . . . . . . 112 4.4.3 Justification of the non-conjugate sampler . . . . . . . . . . . 116 5 Experiments with the MFM 123 5.1 Univariate normal mixtures . . . . . . . . . . . . . . . . . . . . . . . 124 5.1.1 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 5.1.2 Model description . . . . . . . . . . . . . . . . . . . . . . . . . 127 5.1.3 Approximate inference . . . . . . . . . . . . . . . . . . . . . . 128 5.1.4 Density estimation . . . . . . . . . . . . . . . . . . . . . . . . 129 5.1.5 Clustering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 5.1.6 Number of components and clusters . . . . . . . . . . . . . . . 147 5.1.7 Mixing issues . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 5.2 Multivariate skew-normal mixtures . . . . . . . . . . . . . . . . . . . 158 5.2.1 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159 5.2.2 Model description . . . . . . . . . . . . . . . . . . . . . . . . . 162 5.2.3 Approximate inference . . . . . . . . . . . . . . . . . . . . . . 162 5.2.4 Density estimation . . . . . . . . . . . . . . . . . . . . . . . . 163 x 5.2.5 Clustering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 5.2.6 Number of components and clusters . . . . . . . . . . . . . . . 166 6 Combinatorial stochastic processes for other variable-dimension models 169 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170 6.2 Hierarchical mixture of finite mixtures . . . . . . . . . . . . . . . . . 170 6.2.1 A hierarchical variable-dimension mixture model . . . . . . . . 172 6.2.2 Franchise process . . . . . . . . . . . . . . . . . . . . . . . . . 174 6.2.3 Equivalence of the models . . . . . . . . . . . . . . . . . . . . 178 6.2.4 Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179 6.3 Mixture of finite feature models . . . . . . . . . . . . . . . . . . . . . 179 6.3.1 A distribution on binary matrices . . . . . . . . . . . . . . . . 181 6.3.2 A simple urn process . . . . . . . . . . . . . . . . . . . . . . . 181 6.3.3 Equivalence classes of matrices . . . . . . . . . . . . . . . . . 183 6.3.4 Buffet process . . . . . . . . . . . . . . . . . . . . . . . . . . . 186 6.3.5 Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 190 7 Conclusion 191 xi List of Tables 4.1 Summary of known posterior consistency results . . . . . . . . . . . . 99 5.1 RMS of the differences between pairwise probability matrices . . . . . 146 5.2 MFM posterior on number of components for the classic galaxy data 153 5.3 Parameters of the skew-normal mixture used for data simulation . . . 160 xii List of Figures 2.1 Posterior on the number of clusters for a univariate normal DPM . . 7 3.1 DPM posteriors on the number of clusters . . . . . . . . . . . . . . . 22 3.2 Partition sampled from the posterior of a bivariate Gaussian DPM . . 23 3.3 CDF of the size of the first cluster, given two clusters, in a DPM . . . 29 4.1 Graphical model for the MFM . . . . . . . . . . . . . . . . . . . . . . 74 4.2 Alternative graphical models for the MFM . . . . . . . . . . . . . . . 79 4.3 Graphical models for random discrete measures for the MFM . . . . . 87 5.1 Histograms of datasets used . . . . . . . . . . . . . . . . . . . . . . . 126 5.2 Estimated densities for the four component distribution . . . . . . . . 131 5.3 Estimated densities for the Shapley galaxy dataset . . . . . . . . . . . 132 5.4 Estimated densities for classic galaxy, Shapley galaxy, and SLC . . . 133 5.5 Estimated Hellinger distances to the true density . . . . . . . . . . . 134 5.6 Test set log-likelihood for MFM and DPM . . . . . . . . . . . . . . . 136 5.7 Sample clusterings from the posterior, on standard normal data . . . 139 5.8 Sample clusterings from the posterior, on four component data . . . . 140 5.9 Sample clusterings from the posterior, on the classic galaxy data . . . 141 5.10 Sample clusterings from the posterior, on the Shapley galaxy data . . 142 5.11 Sample clusterings from the posterior, on the SLC data . . . . . . . . 143 5.12 Means of the sorted cluster sizes, given the number of clusters . . . . 145 5.13 Pairwise probability matrices . . . . . . . . . . . . . . . . . . . . . . 148 5.14 Posteriors on number of clusters for standard normal data . . . . . . 150 5.15 Posteriors on number of clusters for four component data . . . . . . . 151 5.16 Posteriors on number of clusters for the classic galaxy dataset . . . . 151 5.17 Posteriors on number of clusters for the Shapley galaxy dataset . . . 152 5.18 Posteriors on number of clusters for the SLC dataset . . . . . . . . . 152 5.19 Posteriors on the number of components for the MFM . . . . . . . . . 154 5.20 Traceplots of the number of clusters . . . . . . . . . . . . . . . . . . . 157 5.21 Contour plot and scatterplot of a bivariate skew-normal . . . . . . . . 160 5.22 Skew-normal mixture components used for simulations . . . . . . . . 161 5.23 Contour plots of density estimates for skew-normal mixtures . . . . . 165 5.24 Hellinger distances to the true density for skew-normal mixtures . . . 166 5.25 Sample clusterings from the posterior for skew-normal mixtures . . . 167 5.26 Posteriors of # of components & clusters for skew-normal mixtures . 168 6.1 Graphical models for two representations of the HMFM . . . . . . . . 173 6.2 Graphical model for combinatorial representation of the HMFM . . . 178 6.3 Sample matrices from the MFFM buffet process . . . . . . . . . . . . 188 xiii Chapter One Introduction 2 Nonparametric Bayesian models have found an extraordinarily wide range of appli- cations, including astronomy, meteorology, epidemiology, gene expression profiling, haplotype inference, medical image analysis, survival analysis, econometrics, phylo- genetics, species delimitation, computer vision, classification, document modeling, cognitive science, natural language processing, and perhaps more. Many nonparametric Bayesian models can be viewed as an infinite-dimensional limit of a family of finite-dimensional models. However, another way to construct a flexible Bayesian model is to put a prior on the dimension — that is, to use a variable- dimension model — for example, putting a prior on the number of components in a finite mixture. This thesis (a) analyzes some of the differences and similarities between the nonparametric and variable-dimension approaches, using theory and empirical studies, (b) develops new inference algorithms for these models, and (c) proposes new variable-dimension models and explores their properties. Primarily, we focus on the Dirichlet process mixture (DPM) and a variable- dimensional alternative that we refer to as the mixture of finite mixtures (MFM) model. One of the main differences between DPMs and MFMs is the behavior of the posterior on the number of clusters. In Chapter 2, we use a simple example to show that, on data from a finite mixture, the DPM posterior on the number of clusters may fail to concentrate at the true number of components. Further, in Chapter 3, we generalize this inconsistency result to a large class of nonparametric mixtures, including DPMs and Pitman–Yor process mixtures over a wide range of families of component distributions. On the other hand, it is known that the MFM posterior on the number of components concentrates at the true number, assuming the data comes from a finite mixture with the same family of component distributions as that used in the model (Nobile, 1994). (Note: The material in Chapters 2 and 3 appears in Miller and Harrison (2013a) and Miller and Harrison (2013b), respectively, so 3 there is some overlap in their introductions and definitions; on the other hand, this has the benefit of making these two chapters self-contained.) In Chapter 4, we study the properties of the MFM model, finding that it has many of the same attractive features as the DPM: a simple partition distribution, restau- rant process, random discrete measure representation, and in certain special cases, a simple stick-breaking representation, as well as similar exchangeability properties. As a result, many of the same approximate inference algorithms used for nonpara- metric mixtures can be easily adapted to the MFM. We discuss some of the existing literature on asymptotic results regarding consistency and rates of convergence for the MFM and DPM, for estimating the density, the mixing distribution, and the number of components. We also derive some asymptotic properties of the MFM. Chapter 5 is an empirical study of the posterior behavior of the MFM model, compared with the DPM. We demonstrate similarities and differences with regard to density estimation, clustering, and the posteriors on the number of clusters and components. These experiments use univariate normal mixtures and bivariate skew- normal mixtures on a variety of real and simulated data sets. In Chapter 6, we propose two new variable-dimension models, and derive some of their basic properties. For modeling data in separate but related groups, we propose the hierarchical mixture of finite mixtures (HMFM) as an alternative to the hierar- chical Dirichlet process (HDP) mixture model. For modeling latent binary features, we propose the mixture of finite feature models (MFFM) as an alternative to the Indian buffet process (IBP). As with the MFM, these variable-dimension models ex- hibit some of the same appealing characteristics as their nonparametric counterparts, in particular, simple distributions on discrete structures, exchangeability properties, and representation via combinatorial stochastic processes. Chapter Two A simple example of Dirichlet process mixture inconsistency for the number of components 5 For data assumed to come from a finite mixture with an unknown number of com- ponents, it has become common to use Dirichlet process mixtures (DPMs) not only for density estimation, but also for inferences about the number of components. The typical approach is to use the posterior distribution on the number of clusters — that is, the posterior on the number of components represented in the observed data. However, it turns out that this posterior is not consistent — it does not con- centrate at the true number of components. In this chapter1 , we give an elementary proof of this inconsistency in what is perhaps the simplest possible setting: a DPM with normal components of unit variance, applied to data from a “mixture” with one standard normal component. Further, we show that this example exhibits se- vere inconsistency: instead of going to 1, the posterior probability that there is one cluster converges (in probability) to 0. 2.1 Introduction It is well-known that Dirichlet process mixtures (DPMs) of normals are consistent for the density — that is, given data from a sufficiently regular density p0 the posterior converges to the point mass at p0 (see Ghosal (2010) for details and references). However, it is easy to see that this does not necessarily imply consistency for the number of components, since for example, a good estimate of the density might include superfluous components having vanishingly small weight. Despite the fact that a DPM has infinitely many components with probability 1, it has become common to apply DPMs to data assumed to come from finitely many components or “populations”, and to apply the posterior on the number of clus- 1 The material in this chapter appeared in the proceedings of the NIPS 2013 conference; see Miller and Harrison (2013a). 6 ters (in other words, the number of components used in the process of generating the observed data) for inferences about the true number of components; see Huelsen- beck and Andolfatto (2007), Medvedovic and Sivaganesan (2002), Otranto and Gallo (2002), Xing et al. (2006), Fearnhead (2004) for a few prominent examples. Of course, if the data-generating process very closely resembles the DPM model, then it makes sense to use this posterior for inferences about the number of clusters (but beware of misspecification; see Section 2.2). However, in the examples cited, the authors evaluated the performance of their methods on data simulated from a fixed finite number of components or populations, suggesting that they found this to be more realistic than a DPM for their applications. Therefore, it is important to understand the behavior of this posterior when the data comes from a finite mixture — in particular, does it concentrate at the true number of components? In this chapter, we give a simple example in which a DPM is applied to data from a finite mixture and the posterior distribution on the number of clusters does not concentrate at the true number of components. In fact, DPMs exhibit this type of inconsistency under very general conditions (see Chapter 3) — however, the aim of this chapter is brevity and clarity. To that end, we focus our attention on a special case that is as simple as possible: a “standard normal DPM”, that is, a DPM using univariate normal components of unit variance, with a standard normal base measure (prior on component means). The rest of the chapter is organized as follows. In Section 2.2, we address several pertinent questions and consider some suggestive experimental evidence. In Section 2.3, we formally define the DPM model under consideration. In Section 2.4, we give an elementary proof of inconsistency in the case of a standard normal DPM on data from one component, and in Section 2.5, we show that on standard normal data, a standard normal DPM is in fact severely inconsistent. 7 Figure 2.1: Prior (red x) and estimated posterior (blue o) of the number of clusters in P data, for a univariate normal DPM on n i.i.d. samples from (a) N (0, 1), the observed and (b) 2k=−2 51 N (4k, 12 ). The DPM had concentration parameter α = 1 and a Normal–Gamma base measure on the mean and precision: N (µ | 0, 1/cλ)Gamma(λ | a, b) with a = 1, b = 0.1, and c = 0.001. Estimates were made using a collapsed Gibbs sampler (MacEachern, 1994), with 104 burn-in sweeps and 105 sample sweeps; traceplots and running averages were used as convergence diagnostics. Each plot shown is an average over 5 independent runs. 2.2 Discussion It should be emphasized that these results do not diminish, in any way, the utility of Dirichlet process mixtures as a flexible prior on densities, i.e., for Bayesian density estimation. In addition to their widespread success in empirical studies, DPMs are backed by theoretical guarantees showing that in many cases the posterior on the density concentrates at the true density at the minimax-optimal rate, up to a logarithmic factor (see Ghosal (2010) and references therein). Many researchers (e.g. West et al. (1994), Onogi et al. (2011), among others) have empirically observed that the DPM posterior on the number of clusters tends to 8 overestimate the number of components, in the sense that it tends to put its mass on a range of values greater or equal to the true number. Figure 2.1 illustrates this effect for univariate normals, and similar experiments with different families of component distributions yield similar results. Thus, while our theoretical results in Sections 2.4 and 2.5 (and in Chapter 3) are asymptotic in nature, experimental evidence suggests that the issue is present even in small samples. It is natural to think that this overestimation is due to the fact that the prior on the number of clusters diverges as n → ∞, at a log n rate. However, this does not seem to be the main issue — rather, the problem is that DPMs strongly prefer having some tiny clusters and will introduce extra clusters even when they are not needed (see Chapter 3 for an intuitive explanation of why this is the case). In fact, many researchers have observed the presence of tiny extra clusters (e.g. West et al. (1994), Onogi et al. (2011)), but the reason for this has not previously been well understood, often being incorrectly attributed to the difficulty of detecting components with small weight. These tiny extra clusters are rather inconvenient, especially in clustering applications, and are often dealt with in an ad hoc way by simply removing them. It might be possible to consistently estimate the number of components in this way, but this remains an open question. A more natural solution is the following: if the number of components is un- known, put a prior on the number of components. For example, draw the number of components k from a probability mass function p(k) on {1, 2, . . . } with p(k) > 0 for all k, draw mixing weights π = (π1 , . . . , πk ) (given k), draw component parameters θ1 , . . . , θk i.i.d. (given k and π) from an appropriate prior, and draw X1 , X2 , . . . i.i.d. (given k, π, and θ1:k ) from the resulting mixture. This approach has been widely used (Nobile, 1994, Richardson and Green, 1997, Green and Richardson, 2001, Nobile and 9 Fearnside, 2007). Under certain conditions, the posterior on the density has been shown to concentrate at the true density at the minimax-optimal rate, up to a loga- rithmic factor, for any sufficiently regular true density (Kruijer et al., 2010). Strictly speaking, as defined, such a model is not identifiable, but it is fairly straightforward to modify it to be identifiable by choosing one representative from each equivalence class. Subject to a modification of this sort, it can be shown (see Nobile (1994)) that under very general conditions, when the data is from a finite mixture of the chosen family, such models are (a.e.) consistent for the number of components, the mixing weights, the component parameters, and the density. Also see McCullagh and Yang (2008) for an interesting discussion about estimating the number of components. However, as a practical matter, when dealing with real-world data, one would not expect to find data coming exactly from a finite mixture of a known family (except, perhaps, in rare circumstances). Unfortunately, even for a model as in the preceding paragraph, the posterior on the number of components will typically be highly sensi- tive to misspecification, and it seems likely that in order to obtain robust estimators, the problem itself may need to be reformulated. We urge researchers interested in the number of components to be wary of this robustness issue, and to think carefully about whether they really need to estimate the number of components, or whether some other measure of heterogeneity will suffice. 2.3 Setup In this section, we define the Dirichlet process mixture model under consideration. 10 2.3.1 Dirichlet process mixture model The DPM model was introduced by Ferguson (1983) and Lo (1984) for the pur- pose of Bayesian density estimation, and was made practical through the efforts of several authors (see Escobar and West (1998) and references therein); for a more complete list of references, see Chapter 3. In this chapter, we will use p(·) to denote probabilities under the DPM model (as opposed to other probability distributions that will be considered). The core of the DPM is the so-called Chinese restaurant process (CRP), which defines a certain probability distribution on partitions. Given n ∈ {1, 2, . . . } and t ∈ {1, . . . , n}, let At (n) denote the set of all ordered partitions (A1 , . . . , At ) of {1, . . . , n} into t nonempty sets. In other words, n t [ o At (n) = (A1 , . . . , At ) : A1 , . . . , At are disjoint, Ai = {1, . . . , n}, |Ai | ≥ 1 ∀i . i=1 The CRP with concentration parameter α > 0 defines a probability mass function on A(n) = nt=1 At (n) by setting S t αt Y p(A) = (|Ai | − 1)! α(n) t! i=1 for A ∈ At (n), where α(n) = α(α + 1) · · · (α + n − 1). Note that since t is a function of A, we have p(A) = p(A, t). (It is more common to see this distribution defined in terms of unordered partitions {A1 , . . . , At }, in which case the t! does not appear in the denominator — however, for our purposes it is more convenient to use the distribution on ordered partitions (A1 , . . . , At ) obtained by uniformly permuting the parts. This does not affect the prior or posterior on t.) 11 Consider the hierarchical model t αt Y p(A, t) = p(A) = (|Ai | − 1)!, (2.3.1) α(n) t! i=1 t Y p(θ1:t | A, t) = π(θi ), and i=1 t Y Y p(x1:n | θ1:t , A, t) = pθi (xj ), i=1 j∈Ai where π(θ) is a prior on component parameters θ ∈ Θ, and {pθ : θ ∈ Θ} is a parametrized family of distributions on x ∈ X for the components. Typically, X ⊂ Rd and Θ ⊂ Rk for some d and k. Here, x1:n = (x1 , . . . , xn ) with xi ∈ X , and θ1:t = (θ1 , . . . , θt ) with θi ∈ Θ. This hierarchical model is referred to as a Dirichlet process mixture (DPM) model. P The prior on the number of clusters t under this model is pn (t) = A∈At (n) p(A, t). We use Tn (rather than T ) to denote the random variable representing the number of clusters, as a reminder that its distribution depends on n. Note that we distinguish between the terms “component” and “cluster”: a component is part of a mixture distribution (e.g. a mixture ∞ P i=1 πi pθi has components pθ1 , pθ2 , . . . ), while a cluster is the set of indices of data points coming from a given component (e.g. in the DPM model above, A1 , . . . , At are the clusters). Since we are concerned with the posterior distribution p(Tn = t | x1:n ) on the number of clusters, we will be especially interested in the marginal distribution on 12 (x1:n , t), given by X Z p(x1:n , Tn = t) = p(x1:n , θ1:t , A, t) dθ1:t A∈At (n) X t Z  Y Y  = p(A) pθi (xj ) π(θi ) dθi A∈At (n) i=1 j∈Ai X t Y = p(A) m(xAi ) (2.3.2) A∈At (n) i=1 where for any subset of indices S ⊂ {1, . . . , n}, we denote xS = (xj : j ∈ S) and let m(xS ) denote the single-cluster marginal of xS , Z Y  m(xS ) = pθ (xj ) π(θ) dθ. (2.3.3) j∈S 2.3.2 Specialization to the standard normal case In this chapter, for brevity and clarity, we focus on the univariate normal case with unit variance, with a standard normal prior on means — that is, for x ∈ R and θ ∈ R, 1 pθ (x) = N (x | θ, 1) = √ exp(− 12 (x − θ)2 ), and 2π 1 π(θ) = N (θ | 0, 1) = √ exp(− 12 θ2 ). 2π It is a straightforward calculation to show that the single-cluster marginal is then n 1  1 1  X 2  m(x1:n ) = √ p0 (x1:n ) exp xj , (2.3.4) n+1 2 n + 1 j=1 13 where p0 (x1:n ) = p0 (x1 ) · · · p0 (xn ) (and p0 is the N (0, 1) density). When pθ (x) and π(θ) are as above, we refer to the resulting DPM as a standard normal DPM. 2.4 Simple example of inconsistency In this section, we prove the following result, exhibiting a simple example in which a DPM is inconsistent for the number of components: even when the true number of components is 1 (e.g. N (µ, 1) data), the posterior probability of Tn = 1 does not converge to 1. Interestingly, the result applies even when X1 , X2 , . . . are identically equal to a constant c ∈ R. To keep it simple, we set α = 1; for more general results, see Chapter 3. Theorem 2.4.1. If X1 , X2 , . . . ∈ R are i.i.d. from any distribution with E|Xi | < ∞, then with probability 1, under the standard normal DPM with α = 1 as defined above, p(Tn = 1 | X1:n ) does not converge to 1 as n → ∞. Proof. Let n ∈ {2, 3, . . . }. Let x1 , . . . , xn ∈ R, A ∈ A2 (n), and ai = |Ai | for i = 1, 2. Define sn = nj=1 xj and sAi = j∈Ai xj for i = 1, 2. Using Equation 2.3.4 P P and noting that 1/(n + 1) ≤ 1/(n + 2) + 1/n2 , we have √ m(x1:n )  1 s2  n  1 s2  n  1 s2  n n+1 = exp ≤ exp exp . p0 (x1:n ) 2n+1 2n+2 2 n2 Pn The second factor equals exp( 12 x2n ), where xn = 1 n j=1 xj . By the convexity of x 7→ x2 ,  s 2 a + 1  s 2 a + 1  s 2 n 1 A1 2 A2 ≤ + , n+2 n + 2 a1 + 1 n + 2 a2 + 1 14 and thus, the first factor is less or equal to  1 s2 1 s2A2  √ √ m(xA1 ) m(xA2 ) A1 exp + = a1 + 1 a2 + 1 . 2 a1 + 1 2 a2 + 1 p0 (x1:n ) Hence, √ √ m(x1:n ) a1 + 1 a2 + 1 ≤ √ exp( 12 x2n ). (2.4.1) m(xA1 ) m(xA2 ) n+1 Consequently, we have p(x1:n , Tn = 2) (a) X m(xA1 ) m(xA2 ) = n p(A) p(x1:n , Tn = 1) m(x1:n ) A∈A2 (n) √ (b) X n+1 ≥ n p(A) p p exp(− 21 x2n ) A∈A2 (n) |A1 | + 1 |A 2 | + 1 √ (c) X (n − 2)! n + 1 ≥ n √ √ exp(− 12 x2n ) n! 2! 2 n A∈A (n): 2 |A1 |=1 (d) 1 ≥ √ exp(− 12 x2n ), 2 2 where step (a) follows from applying Equation 2.3.2 to both numerator and de- nominator, plus using Equation 2.3.1 (with α = 1) to see that p(A) = 1/n when A = ({1, . . . , n}), step (b) follows from Equation 2.4.1 above, step (c) follows since all the terms in the sum are nonnegative and p(A) = (n − 2)!/n! 2! when |A1 | = 1 (by Equation 2.3.1, with α = 1), and step (d) follows since there are n partitions A ∈ A2 (n) such that |A1 | = 1. If X1 , X2 , . . . ∈ R are i.i.d. with µ = EXj finite, then by the law of large numbers, 15 1 Pn Xn = n j=1 Xj → µ almost surely as n → ∞. Therefore, p(X1:n , Tn = 1) p(X1:n , Tn = 1) p(Tn = 1 | X1:n ) = P∞ ≤ t=1 p(X1:n , Tn = t) p(X1:n , Tn = 1) + p(X1:n , Tn = 2) 1 a.s. 1 ≤ 2 −−→ < 1. 1 + 2√1 2 exp(− 12 X n ) 1 1 + 2√2 exp(− 12 µ2 ) Hence, almost surely, p(Tn = 1 | X1:n ) does not converge to 1. 2.5 Severe inconsistency In the previous section, we showed that p(Tn = 1 | X1:n ) does not converge to 1 for a standard normal DPM on any data with finite mean. In this section, we prove that in fact, it converges to 0 on standard normal data. This vividly illustrates that improperly using DPMs in this way can lead to entirely misleading results. The key step in the proof is an application of Hoeffding’s strong law of large numbers for U-statistics. Theorem 2.5.1. If X1 , X2 , . . . ∼ N (0, 1) i.i.d. then Pr p(Tn = 1 | X1:n ) −→ 0 as n → ∞ under the standard normal DPM with concentration parameter α = 1. Proof. For t = 1 and t = 2 define p(X1:n , Tn = t) Rt (X1:n ) = n3/2 . p0 (X1:n ) 16 Our method of proof is as follows. We will show that Pr R2 (X1:n ) −−−→ ∞ n→∞ (or in other words, for any B > 0 we have P(R2 (X1:n ) > B) → 1 as n → ∞), and we will show that R1 (X1:n ) is bounded in probability: R1 (X1:n ) = OP (1) (or in other words, for any ε > 0 there exists Bε > 0 such that P(R1 (X1:n ) > Bε ) ≤ ε for all n ∈ {1, 2, . . . }). Putting these two together, we will have p(X1:n , Tn = 1) p(X1:n , Tn = 1) R1 (X1:n ) Pr p(Tn = 1 | X1:n ) = P∞ ≤ = −−−→ 0. t=1 p(X1:n , Tn = t) p(X1:n , Tn = 2) R2 (X1:n ) n→∞ First, let’s show that R2 (X1:n ) → ∞ in probability. For S ⊂ {1, . . . , n} with |S| ≥ 1, define h(xS ) by m(xS ) 1  1 1  X 2  h(xS ) = =p exp xj , p0 (xS ) |S| + 1 2 |S| + 1 j∈S where m is the single-cluster marginal as in Equations 2.3.3 and 2.3.4. Note that √ when 1 ≤ |S| ≤ n − 1, we have n h(xS ) ≥ 1. Note also that Eh(XS ) = 1 since Z Z Eh(XS ) = h(xS ) p0 (xS ) dxS = m(xS ) dxS = 1, using the fact that m(xS ) is a density with respect to Lebesgue measure. For k ∈ {1, . . . , n}, define the U-statistics 1 X Uk (X1:n ) = n  h(XS ) k |S|=k 17 where the sum is over all S ⊂ {1, . . . , n} such that |S| = k. By Hoeffding’s strong law of large numbers for U-statistics (Hoeffding, 1961), a.s. Uk (X1:n ) −−−→ Eh(X1:k ) = 1 n→∞ for any k ∈ {1, 2, . . . }. Therefore, using Equations 2.3.1 and 2.3.2 we have that for any K ∈ {1, 2, . . . } and any n > K, X m(XA1 ) m(XA2 ) R2 (X1:n ) = n3/2 p(A) p0 (X1:n ) A∈A2 (n) X √ =n p(A) n h(XA1 ) h(XA2 ) A∈A2 (n) X ≥n p(A) h(XA1 ) A∈A2 (n) n−1 X X (k − 1)! (n − k − 1)! =n h(XS ) k=1 n! 2! |S|=k n−1 X n 1 X = h(XS ) 2k(n − k) nk  k=1 |S|=k n−1 X n = Uk (X1:n ) k=1 2k(n − k) K X n ≥ Uk (X1:n ) k=1 2k(n − k) K a.s. X 1 HK log K −−−→ = > n→∞ k=1 2k 2 2 where HK is the K th harmonic number, and the last inequality follows from the standard bounds (Graham et al., 1989) on harmonic numbers: log K < HK ≤ log K + 1. Hence, for any K, log K lim inf R2 (X1:n ) > almost surely, n→∞ 2 18 and it follows easily that a.s. R2 (X1:n ) −−−→ ∞. n→∞ Convergence in probability is implied by almost sure convergence. Now, let’s show that R1 (X1:n ) = OP (1). By Equations 2.3.1, 2.3.2, and 2.3.4, we have p(X1:n , Tn = 1) √ m(X1:n ) R1 (X1:n ) = n3/2 = n p0 (X1:n ) p0 (X1:n ) √ n n 1 n  1 X 2  =√ exp √ Xi ≤ exp(Zn2 /2) n+1 2n+1 n i=1 √ P where Zn = (1/ n) ni=1 Xi ∼ N (0, 1) for each n ∈ {1, 2, . . . }. Since Zn = OP (1) then we conclude that R1 (X1:n ) = OP (1). This completes the proof. Chapter Three Inconsistency of Pitman–Yor process mixtures for the number of components 20 3.1 Introduction In many applications, a finite mixture is a natural model, but it can be difficult to choose an appropriate number of components. To circumvent this choice, investiga- tors are increasingly turning to Dirichlet process mixtures (DPMs), and Pitman–Yor process mixtures (PYMs), more generally. While these models may be well-suited for Bayesian density estimation, some investigators are using them for inferences about the number of components, by considering the posterior on the number of compo- nents represented in the observed data. We show that this posterior is not consistent — that is, on data from a finite mixture, it does not concentrate at the true number of components. This result applies to a large class of nonparametric mixtures, in- cluding DPMs and PYMs, over a wide variety of families of component distributions, including essentially all discrete families, as well as continuous exponential families satisfying mild regularity conditions (such as multivariate Gaussians).1 3.1.1 A motivating example In population genetics, determining the “population structure” is an important step in the analysis of sampled data. As an illustrative example, consider the impala, a species of antelope in southern Africa. Impalas are divided into two subspecies: the common impala occupying much of the eastern half of the region, and the black-faced impala inhabiting a small area in the west. While common impalas are abundant, the number of black-faced impalas has been decimated by drought, poaching, and declining resources due to human and livestock expansion. To assist conservation efforts, Lorenzen et al. (2006) collected samples from 216 impalas, and analyzed the 1 The material in this chapter appears in a manuscript that has been submitted for publication; see Miller and Harrison (2013b). 21 genetic variation between/within the two subspecies. A key part of their analysis consisted of inferring the population structure — that is, partitioning the data into distinct populations, and in particular, determining how many such populations there are. To infer the impala population structure, Lorenzen et al. employed a widely-used tool called Structure (Pritchard et al., 2000) which, in the simplest version, models the data as a finite mixture, with each component in the mixture corresponding to a distinct population. Structure uses an ad hoc method to choose the number of components, but this comes with no guarantees. Seeking a more principled approach, Pella and Masuda (2006) proposed using a Dirichlet process mixture (DPM). Now, in a DPM, the number of components is infinite with probability 1, and thus the posterior on the number of components is always, trivially, a point mass at infinity. Consequently, Pella and Masuda instead employed the posterior on the number of clusters (that is, the number of compo- nents used in generating the data observed so far) for inferences about the number of components. (The terms “component” and “cluster” are often used interchange- ably, but we make the following crucial distinction: a component is part of a mixture distribution, while a cluster is the set of indices of datapoints coming from a given component.) This DPM approach was implemented in a software tool called Struc- turama (Huelsenbeck and Andolfatto, 2007), and demonstrated on the impala data of Lorenzen et al.; see Figure 3.1(a). Structurama has gained acceptance within the population genetics commu- nity, and has been used in studies of a variety of organisms, from apples and av- ocados, to sardines and geckos (Richards et al., 2009, Chen et al., 2009, Gonzalez and Zardoya, 2007, Leach´e and Fujita, 2010). Studies such as these can carry signif- icant weight, since they may be used by officials to make informed policy decisions 22 (a) Posterior for impala data (b) Posterior for Gaussian data Figure 3.1: Estimated DPM posterior distribution of the number of clusters: (a) For the impala data of Lorenzen et al. (n = 216 datapoints). Our empirical results, shown here, agree with those of Huelsenbeck and Andolfatto. (b) For bivariate Gaussian data from a four-component mixture; see Figure 3.2. Each plot is the average over 10 independently-drawn datasets. (Lines drawn for visualization purposes only.) (For (a) and (b), estimates were made via Gibbs sampling, with 104 burn-in sweeps and 105 sample sweeps.) regarding agriculture, conservation, and public health. More generally, in a number of applications the same scenario has played out: a finite mixture seems to be a natural model, but requires the user to choose the number of components, while a Dirichlet process mixture offers a convenient way to avoid this choice. For nonparametric Bayesian density estimation, DPMs are indeed attractive, since the posterior on the density exhibits nice convergence properties; see Section 3.1.3. However, in several applications, investigators have drawn inferences from the posterior on the number of clusters — not just the density — on the assumption that this is informative about the number of components. Further examples include gene expression profiling (Medvedovic and Sivaganesan, 2002), haplotype inference (Xing et al., 2006), econometrics (Otranto and Gallo, 2002), and evaluation of inference algorithms (Fearnhead, 2004). Of course, if the data-generating process is well- modeled by a DPM (and in particular, there are infinitely many components), then it is sensible to use this posterior for inference about the number of components represented so far in the data — but that does not seem to be the perspective of these investigators, since they measure performance on simulated data coming from 23 Figure 3.2: A typical partition sampled from the posterior of a Dirichlet process mixture of bivariate Gaussians, on simulated data from a four-component mixture. Different clusters have different marker shapes (+, ×, O, M,◦, 2) and different colors. Note the tiny “extra” clusters (◦ and 2), in addition to the four dominant clusters. finitely many components or populations. (To be clear, we should note that although several investigators have used DPMs in this way, most of those who are familiar with nonparametric models are well aware that this usage is questionable.) Therefore, it is important to understand the properties of this procedure. Sim- ulation results give some cause for concern; for instance, Figures 3.1(b) and 3.2 display results for data from a mixture of two-dimensional Gaussians with four com- ponents. Partitions sampled from the posterior often have tiny “extra” clusters, and the posterior on the number of clusters does not appear to be concentrating as the number of datapoints n increases. This raises a fundamental question that has not been addressed in the literature: With enough data, will this posterior eventually concentrate at the true number of components? In other words, is it consistent? 24 3.1.2 Overview of results In this chapter, we prove that under fairly general conditions, when using a Dirichlet process mixture, the posterior on the number of clusters will not concentrate at any finite value, and therefore will not be consistent for the number of components in a finite mixture. In fact, our results apply to a large class of nonparametric mixtures including DPMs, and Pitman–Yor process mixtures (PYMs) more generally, over a wide variety of families of component distributions. Before treating our general results and their prerequisite technicalities, we would like to highlight a few interesting special cases that can be succinctly stated. The terminology and notation used below will be made precise in later sections. To reiterate, our results are considerably more general than the following corollary, which is simply presented for the reader’s convenience. Corollary 3.1.1. Consider a Pitman–Yor process mixture with component distribu- tions from one of the following families: (a) Normal(µ, Σ) (multivariate Gaussian), (b) Exponential(θ), (c) Gamma(a, b), (d) Log-Normal(µ, σ 2 ), or (e) Weibull(a, b) with fixed shape a > 0, along with a base measure that is a conjugate prior of the form in Section 3.5.2, or 25 T (f ) any discrete family {Pθ } such that θ {x 6 ∅ (e.g., Poisson, : Pθ (x) > 0} = Geometric, Negative Binomial, Binomial, Multinomial, etc.), along with any continuous base measure. Consider any t ∈ {1, 2, . . . }, except for t = N in the case of a Pitman–Yor process with parameters σ < 0 and ϑ = N |σ|. If X1 , X2 , . . . are i.i.d. from a mixture with t components from the family used in the model, then the posterior on the number of clusters Tn is not consistent for t, and in fact, lim sup p(Tn = t | X1 , . . . , Xn ) < 1 n→∞ with probability 1. This is implied by Theorems 3.3.4, 3.4.1, and 3.6.2. These more general theorems apply to a broad class of partition distributions, handling Pitman–Yor processes as a special case, and they apply to many other families of component distributions: Theorem 3.6.2 covers a large class of exponential families, and Theorem 3.4.1 covers families satisfying a certain boundedness condition on the densities (including any case in which the model and data distributions have one or more point masses in common, as well as many location–scale families with scale bounded away from zero). Dirichlet processes are subsumed as a further special case, being Pitman–Yor processes with parameters σ = 0 and ϑ > 0. Also, the assumption of i.i.d. data from a finite mixture is much stronger than what is required by these results. Regarding the exception of t = N when σ < 0 in Corollary 3.1.1: posterior con- sistency at t = N is possible, however, this could only occur if the chosen parameter N just happens to be equal to the actual number of components, t. On the other hand, consistency at any t can (in principle) be obtained by putting a prior on N ; see Section 3.1.3 below. In a similar vein, some investigators place a prior on the 26 concentration parameter ϑ in a DPM, or allow ϑ to depend on n; we conjecture that inconsistency can still occur in these cases, but in this work, we examine only the case of fixed σ and ϑ. 3.1.3 Discussion / related work We would like to emphasize that this inconsistency should not be viewed as a defi- ciency of Dirichlet process mixtures, but is simply due to a misapplication of them. As flexible priors on densities, DPMs are superb, and there are strong results showing that in many cases the posterior on the density converges in L1 to the true density at the minimax-optimal rate, up to a logarithmic factor (see Scricciolo (2012), Ghosal (2010) and references therein). Further, Nguyen (2013) has recently shown that the posterior on the mixing distribution converges in the Wasserstein metric to the true mixing distribution. However, these results do not necessarily imply consistency for the number of components, since any mixture can be approximated arbitrarily well in these metrics by another mixture with a larger number of components (for in- stance, by making the weights of the extra components infinitesimally small). There seems to be no prior work on consistency of DPMs (or PYMs) for the number of components in a finite mixture (aside from Miller and Harrison (2013a) — appearing here as Chapter 2 — in which we discuss the very special case of a DPM on data from a univariate Gaussian “mixture” with one component of known variance). In the context of “species sampling”, several authors have studied the Pitman– Yor process posterior (see James (2008), Jang et al. (2010), Lijoi et al. (2007) and references therein), but this is very different from our situation — in a species sam- pling model, the observed data is drawn directly from a measure with a Pitman–Yor process prior, while in a PYM model, the observed data is drawn from a mixture 27 with such a measure as the mixing distribution. Rousseau and Mengersen (2011) proved an interesting result on “overfitted” mix- tures, in which data from a finite mixture is modeled by a finite mixture with too many components. In cases where this approximates a DPM, their result implies that the posterior weight of the extra components goes to zero. In a rough sense, this is complementary to our results, which involve showing that there are always some nonempty (but perhaps small) extra clusters. Empirically, many investigators have noticed that the DPM posterior tends to overestimate the number of components (e.g. West et al. (1994), Lartillot and Philippe (2004), Onogi et al. (2011)), and such observations are consistent with our theoretical results. This overestimation seems to occur because there are typ- ically a few tiny “extra” clusters. Among researchers using DPMs for clustering, this is an annoyance that is often dealt with by pruning such clusters — that is, by simply ignoring them when calculating statistics such as the number of clusters. It may be possible to obtain consistent estimators in this way, but this remains an open question; Rousseau and Mengersen’s (2011) results may be applicable here. Under the (strong) assumption that the family of component distributions is correctly specified, one can obtain posterior consistency by simply putting a prior on the number of components in a finite mixture (Nobile, 1994). (It turns out that putting a prior on N in a PYM with σ < 0, ϑ = N |σ| is a special case of this (Gnedin and Pitman, 2006).) That said, it seems likely that such estimates will be severely affected by misspecification of the model, which is inevitable in most applications. Robustness to model misspecification seems essential for reliable estimation of the number of components, for real-world data. Other approaches have also been proposed for estimating the number of components (Henna, 1985, Keribin, 28 2000, Leroux, 1992, Ishwaran et al., 2001, James et al., 2001, Henna, 2005), and steps toward addressing the issue of robustness have been taken by Woo and Sriram (2006, 2007). 3.1.4 Intuition for the result To illustrate the intuition behind this inconsistency, consider a Dirichlet process with concentration parameter ϑ = 1. (Similar reasoning applies for any Pitman–Yor process with σ ≥ 0, but the σ < 0 case is somewhat different.) It is tempting to think that the prior on the number of clusters is the culprit, since (as is well-known) it diverges as n → ∞. Surprisingly, this does not seem to be the main reason why inconsistency occurs. Instead, the right intuition comes from examining the prior on partitions, given the number of clusters. The prior on ordered partitions A = (A1 , . . . , At ) is p(A) = (n! t!)−1 ti=1 (ai −1)!, where t is the number of parts (i.e. clusters) and ai = |Ai | is the Q size of the ith part. (The t! comes from uniformly permuting the parts; see Section 3.2.1.) Since there are n!/(a1 ! · · · at !) such partitions with part sizes (a1 , . . . , at ), the conditional distribution of the sizes (a1 , . . . , at ) given t is proportional to a−1 −1 1 · · · at P (subject to the constraint that ai = n). See Figure 3.3 for the case of t = 2. The key observation is that, for large n, this conditional distribution is heavily concentrated in the “corners”, where one or more of the ai ’s is small. By pursuing this line of thought, one can show that the probability of drawing a partition with t + 1 parts and one or more of the ai ’s equal to 1 is, at least, the same order of magnitude (with respect to n) as the probability of drawing a partition with t parts. This leads to the basic idea of the proof — if the likelihood of the data is 29 Figure 3.3: The cumulative distribution function of the conditional distribution of a1 given that t = 2, for a Dirichlet process with ϑ = 1. As n increases, the distribution becomes concentrated at the extremes. also the same order of magnitude, then the posterior probability of t + 1 will not be too much smaller than that of t. Roughly speaking, the posterior will always find it reasonably attractive to split off one element as a singleton. Handling the likelihood is the difficult part, and occupies the bulk of the proof. 3.1.5 Organization of this chapter In Section 3.2, we define Gibbs partition mixture models, which includes Pitman–Yor and Dirichlet process mixtures as special cases. In Section 3.3, we prove a general inconsistency theorem for Gibbs partition mixtures satisfying certain conditions. In Section 3.4, we apply the theorem to cases satisfying a certain boundedness condition on the densities, including discrete families as a special case. In Section 3.5, we introduce notation for exponential families and conjugate priors, and in Section 3.6, we apply the theorem to cases in which the mixture is over an exponential family satisfying some regularity conditions. The remainder of the chapter proves the key lemma used in this application. In Section 3.7, we obtain certain inequalities involving the marginal density under an exponential family with conjugate prior. In Section 3.8, we prove the key lemma of Section 3.6: an inequality involving the marginal density of any sufficiently large subset of the data. Sections 3.9, 3.10, 3.11, and 3.12 contain various supporting results. 30 3.2 Model distribution Our analysis involves two probability distributions: one which is defined by the model, and another which gives rise to the data. In this section, we describe the model distribution. Building upon the Dirichlet process (Ferguson, 1973, Blackwell and MacQueen, 1973, Antoniak, 1974), Dirichlet process mixtures were first studied by Antoniak (1974) (who also considered the mixtures of Dirichlet processes arising in the poste- rior), Berry and Christensen (1979), Ferguson (1983), and Lo (1984), and were later made practical through the efforts of a number of authors (Escobar, 1988, MacEach- ern, 1994, Escobar and West, 1995, West, 1992, West et al., 1994, Neal, 1992, Liu, 1994, Bush and MacEachern, 1996, MacEachern and M¨ uller, 1998, MacEachern, 1998, Escobar and West, 1998, MacEachern et al., 1999, Neal, 2000). Pitman–Yor process mixtures (Ishwaran and James, 2003) are a generalization of DPMs based on the Pitman–Yor process (Perman et al., 1992, Pitman and Yor, 1997), also known as the two-parameter Poisson–Dirichlet process. The model we consider is, in turn, a generalization of PYMs based on the family of Gibbs partitions (Pitman, 2006, Gnedin and Pitman, 2006). 3.2.1 Gibbs partitions We will use p(·) to denote probabilities and probability densities under the model. Our model specification begins with a distribution on partitions, or more precisely, on ordered partitions. Given n ∈ {1, 2, . . . } and t ∈ {1, . . . , n}, let At (n) denote the set of all ordered partitions (A1 , . . . , At ) of {1, . . . , n} into t nonempty sets (or 31 “parts”). In other words, n t [ o At (n) = (A1 , . . . , At ) : A1 , . . . , At are disjoint, Ai = {1, . . . , n}, |Ai | ≥ 1 ∀i . i=1 Sn For each n ∈ {1, 2, . . . }, consider a probability mass function (p.m.f.) on t=1 At (n) of the form t Y p(A) = vn (t) wn (|Ai |) (3.2.1) i=1 for A ∈ At (n), where vn : {1, . . . , n} → [0, ∞) and wn : {1, . . . , n} → [0, ∞). This induces a distribution on t in the natural way, via p(t | A) = I(A ∈ At (n)). (Throughout, we use I to denote the indicator function: I(E) is 1 if E is true, and 0 otherwise.) It follows that p(A) = p(A, t) when A ∈ At (n). Although it is more common to use a distribution on unordered partitions {A1 , . . . , At }, for our purposes it is more convenient to work with the corresponding distribution on ordered partitions (A1 , . . . , At ) obtained by uniformly permuting the parts. This does not affect the distribution of t. Under this correspondence, any p.m.f. as in Equation 3.2.1 corresponds to a member of the class of “exchangeable partition probability functions”, or EPPFs (Pitman, 2006). In particular, for any given n it yields an EPPF in “Gibbs form”, and a random partition from such an EPPF is called a Gibbs partition (Pitman, 2006). (Note: We do not assume that, as n varies, the sequence of p.m.f.s in Equation 3.2.1 necessarily satisfies the marginal- ization property referred to as “consistency in distribution”.) For example, to obtain the partition distribution for a Dirichlet process, we can 32 choose ϑt vn (t) = and wn (a) = (a − 1)! (3.2.2) ϑn↑1 t! where ϑ > 0 and xn↑δ = x(x+δ)(x+2δ) · · · (x+(n−1)δ), with x0↑δ = 1 by convention. (The t! in the denominator appears since we are working with ordered partitions.) More generally, to obtain the partition distribution for a Pitman–Yor process, we can choose (ϑ + σ)t−1↑σ vn (t) = and wn (a) = (1 − σ)a−1↑1 (3.2.3) (ϑ + 1)n−1↑1 t! where either σ ∈ [0, 1) and ϑ ∈ (−σ, ∞), or σ ∈ (−∞, 0) and ϑ = N |σ| for some N ∈ {1, 2, . . . } (Ishwaran and James, 2003). When σ = 0, this reduces to the partition distribution of a Dirichlet process. When σ < 0 and ϑ = N |σ|, it is the partition distribution obtained by drawing q = (q1 , . . . , qN ) from a symmetric N - dimensional Dirichlet with parameters |σ|, . . . , |σ|, sampling assignments Z1 , . . . , Zn i.i.d. from q, and removing any empty parts (Gnedin and Pitman, 2006). Thus, in this latter case, t is always in {1, . . . , N }. 33 3.2.2 Gibbs partition mixtures Consider the hierarchical model t Y p(A, t) = p(A) = vn (t) wn (|Ai |), (3.2.4) i=1 t Y p(θ1:t | A, t) = π(θi ), i=1 t Y Y p(x1:n | θ1:t , A, t) = pθi (xj ), i=1 j∈Ai where π is a prior density on component parameters θ ∈ Θ ⊂ Rk for some k, and {pθ : θ ∈ Θ} is a parametrized family of densities on x ∈ X ⊂ Rd for some d. Here, x1:n = (x1 , . . . , xn ) with xi ∈ X , θ1:t = (θ1 , . . . , θt ) with θi ∈ Θ, and A ∈ At (n). Assume that π is a density with respect to Lebesgue measure, and that {pθ : θ ∈ Θ} are densities with respect to some sigma-finite Borel measure λ on X , such that (θ, x) 7→ pθ (x) is measurable. (Of course, the distribution of x under pθ (x) may be discrete, continuous, or neither, depending on the nature of λ.) For x1 , . . . , xn ∈ X and J ⊂ {1, . . . , n}, define the single-cluster marginal, Z Y  m(xJ ) = pθ (xj ) π(θ) dθ, (3.2.5) Θ j∈J where xJ = (xj : j ∈ J), and assume m(xJ ) < ∞. By convention, m(xJ ) = 1 when J = ∅. Note that m(xJ ) is a density with respect to the product measure λ` on X ` , where ` = |J|, and that m(xJ ) can (and often will) be positive outside the support of λ` . Definition 3.2.1. We refer to such a hierarchical model as a Gibbs partition mixture model. 34 (Note: This is, perhaps, a slight abuse of the term “Gibbs partition”, since we allow vn and wn to vary arbitrarily with n.) In particular, it is a Dirichlet process mixture model when vn and wn are as in Equation 3.2.2, or more generally, a Pitman– Yor process mixture model when vn and wn are as in Equation 3.2.3. We distinguish between the terms “component” and “cluster”: a component of a mixture is one of the distributions used in it (e.g. pθi ), while a cluster is the set of indices of datapoints coming from a given component (e.g. Ai ). The prior on P the number of clusters under such a model is pn (t) = A∈At (n) p(A). We use Tn , rather than T , to denote the random variable representing the number of clusters, as a reminder that its distribution depends on n. Since we are concerned with the posterior p(Tn = t | x1:n ) on the number of clusters, we will be especially interested in the marginal density of (x1:n , t), given by X Z p(x1:n , Tn = t) = p(x1:n , θ1:t , A, t) dθ1:t A∈At (n) X t Z  Y Y  = p(A) pθi (xj ) π(θi ) dθi A∈At (n) i=1 j∈Ai X t Y = p(A) m(xAi ). (3.2.6) A∈At (n) i=1 As usual, the posterior p(Tn = t | x1:n ) is not uniquely defined, since it can be modified arbitrarily on any subset of X n having probability zero under the model distribution. For definiteness, we will employ the usual version of this posterior, p(x1:n , Tn = t) p(x1:n , Tn = t) p(Tn = t | x1:n ) = = P∞ 0 p(x1:n ) t0 =1 p(x1:n , Tn = t ) whenever the denominator is nonzero, and p(Tn = t | x1:n ) = 0 otherwise (for 35 notational convenience). 3.3 Inconsistency theorem The essential ingredients in the main theorem are Conditions 3.3.1 and 3.3.2 below. For each n ∈ {1, 2, . . . }, consider a partition distribution as in Equation 3.2.1. For n > t ≥ 1, define wn (a) vn (t) cwn = max and cvn (t) = , a∈{2,...,n} a wn (a − 1) wn (1) vn (t + 1) with the convention that 0/0 = 0 and y/0 = ∞ for y > 0. Condition 3.3.1. Assume lim supn→∞ cwn < ∞ and lim supn→∞ cvn (t) < ∞, given some particular t ∈ {1, 2, . . . }. For Pitman–Yor processes, Condition 3.3.1 holds for all relevant values of t, by Proposition 3.3.3 below. Now, consider a collection of single-cluster marginals m(·) as in Equation 3.2.5. Given n ≥ t ≥ 1, x1 , . . . , xn ∈ X , and c ∈ [0, ∞), define 1 ϕt (x1:n , c) = min |SA (x1:n , c)| A∈At (n) n where SA (x1:n , c) is the set of indices j ∈ {1, . . . , n} such that the part A` containing j satisfies m(xA` ) ≤ c m(xA` rj )m(xj ). Condition 3.3.2. Given a sequence of random variables X1 , X2 , . . . ∈ X , a collection of single-cluster marginals m(·), and t ∈ {1, 2, . . . }, assume sup lim inf ϕt (X1:n , c) > 0 with probability 1. c∈[0,∞) n→∞ 36 Note that Condition 3.3.1 involves only the partition distributions, while Condi- tion 3.3.2 involves only the data distribution and the single-cluster marginals. Proposition 3.3.3. Consider a Pitman–Yor process. If σ ∈ [0, 1) and ϑ ∈ (−σ, ∞) then Condition 3.3.1 holds for any t ∈ {1, 2, . . . }. If σ ∈ (−∞, 0) and ϑ = N |σ|, then it holds for any t ∈ {1, 2, . . . } except N . Proof. This is a simple calculation. See Section 3.9. Theorem 3.3.4. Let X1 , X2 , . . . ∈ X be a sequence of random variables (not neces- sarily i.i.d.). Consider a Gibbs partition mixture model. For any t ∈ {1, 2, . . . }, if Conditions 3.3.1 and 3.3.2 hold, then lim sup p(Tn = t | X1:n ) < 1 with probability 1. n→∞ If, further, the sequence X1 , X2 , . . . is i.i.d. from a mixture with t components, then with probability 1 the posterior of Tn (under the model) is not consistent for t. Proof. This follows easily from Lemma 3.3.5 below. See Section 3.9. Lemma 3.3.5. Consider a Gibbs partition mixture model. Let n > t ≥ 1, x1 , . . . , xn ∈ X , and c ∈ [0, ∞). If ϕt (x1:n , c) > t/n, cwn < ∞, and cvn (t) < ∞, then Ct (x1:n , c) p(Tn = t | x1:n ) ≤ , 1 + Ct (x1:n , c) where Ct (x1:n , c) = t c cwn cvn (t)/(ϕt (x1:n , c) − t/n). Proof. To simplify notation, let us denote ϕ = ϕt (x1:n , c), C = Ct (x1:n , c), and SA = SA (x1:n , c) for A ∈ At (n). Given J ⊂ {1, . . . , n} such that |J| ≥ 1, define hJ = wn (|J|) m(xJ ). 37 S  For A ∈ At (n), let RA = SA r i:|Ai |=1 Ai , that is, RA consists of those j ∈ SA such that the size of the part A` containing j is greater than 1. Note that |RA | ≥ |SA | − t ≥ nϕ − t > 0. (3.3.1) For any j ∈ RA , the part A` containing j satisfies hA` = wn (|A` |) m(xA` ) (3.3.2) ≤ cwn |A` | wn (|A` | − 1) wn (1) c m(xA` rj ) m(xj ) ≤ n c cwn hA` rj hj . Given j ∈ RA , define B(A, j) to be the element B of At+1 (n) such that Bi = Ai r j for i = 1, . . . , t, and Bt+1 = {j} (that is, remove j from whatever part it belongs to, and make {j} the (t + 1)th part). Define  YA = B(A, j) : j ∈ RA . Now, using Equations 3.3.1 and 3.3.2, for any A ∈ At (n) we have t t Y 1 X X Y hAi = hA` hAi (3.3.3) i=1 |RA | `=1 j∈R ∩A i6=` A ` t 1 X X Y ≤ n c cwn hA` rj hj hAi nϕ − t `=1 j∈RA ∩A` i6=` t+1 c cwn X Y = hBi (A,j) ϕ − t/n j∈R i=1 A t+1 c cwn X hY i = hBi I(B ∈ YA ). ϕ − t/n i=1 B∈At+1 (n) 38 For any B ∈ At+1 (n),  # A ∈ At (n) : B ∈ YA ≤ t, (3.3.4) since there are only t parts that Bt+1 could have come from. Therefore, t (a) X Y p(x1:n , Tn = t) = p(A) m(xAi ) A∈At (n) i=1 t (b) X Y = vn (t) hAi A∈At (n) i=1 t+1 (c) c cwn X X hY i ≤ vn (t) hBi I(B ∈ YA ) ϕ − t/n i=1 A∈At (n) B∈At+1 (n) t+1 c cwn X hY i  = vn (t) hBi # A ∈ At (n) : B ∈ YA ϕ − t/n i=1 B∈At+1 (n) t+1 (d) c cwn cvn (t) X hY i ≤ vn (t + 1) hBi t ϕ − t/n i=1 B∈At+1 (n) t+1 t c cwn cvn (t) X Y = p(B) m(xBi ) ϕ − t/n i=1 B∈At+1 (n) = C p(x1:n , Tn = t + 1), where (a) is from Equation 3.2.6, (b) is from Equation 3.2.4 and the definition of hJ above, (c) follows from Equation 3.3.3, and (d) follows from Equation 3.3.4. If p(Tn = t | x1:n ) = 0, then trivially p(Tn = t | x1:n ) ≤ C/(C + 1). On the other hand, if p(Tn = t | x1:n ) > 0, then p(x1:n , Tn = t) > 0, and therefore p(x1:n , Tn = t) p(Tn = t | x1:n ) = P∞ 0 t0 =1 p(x1:n , Tn = t ) p(x1:n , Tn = t) C ≤ ≤ . p(x1:n , Tn = t) + p(x1:n , Tn = t + 1) C +1 39 3.4 Application to discrete or bounded cases By Theorem 3.3.4, the following result implies inconsistency in a large class of PYM models, including essentially all discrete cases (or more generally anything with at least one point mass) and a number of continuous cases as well. Theorem 3.4.1. Consider a family of densities {pθ : θ ∈ Θ} on X along with a prior π on Θ and the resulting collection of single-cluster marginals m(·) as in Equation 3.2.5. Let X1 , X2 , . . . ∈ X be a sequence of random variables (not necessarily i.i.d.). If there exists U ⊂ X such that n 1X (1) lim inf I(Xj ∈ U ) > 0 with probability 1, and n→∞ n j=1 n p (x) o θ (2) sup : x ∈ U, θ ∈ Θ < ∞ (where 0/0 = 0, y/0 = ∞ for y > 0), m(x) then Condition 3.3.2 holds for all t ∈ {1, 2, . . . }. Proof. Suppose U ⊂ X satisfies (1) and (2), and let t ∈ {1, 2, . . . }. Define c = sup pm(x)  θ (x) : x ∈ U, θ ∈ Θ . Let n > t and x1 , . . . , xn ∈ X . Now, for any x ∈ U and θ ∈ Θ, we have pθ (x) ≤ c m(x). Hence, for any J ⊂ {1, . . . , n}, if j ∈ J and xj ∈ U then Z h Y i m(xJ ) = pθ (xj ) pθ (xi ) π(θ) dθ ≤ c m(xj )m(xJrj ). (3.4.1) Θ i∈Jrj  Thus, letting R(x1:n ) = j ∈ {1, . . . , n} : xj ∈ U , we have R(x1:n ) ⊂ SA (x1:n , c) for any A ∈ At (n), and hence, ϕt (x1:n , c) ≥ n1 |R(x1:n )|. 40 Therefore, by (1), with probability 1, 1 lim inf ϕt (X1:n , c) ≥ lim inf |R(X1:n )| > 0. n→∞ n→∞ n The preceding theorem covers a fairly wide range of cases; here are some exam- ples. Consider a model with {pθ }, π, λ, and m(·), as in Section 3.2. (i) Finite sample space. Suppose X is a finite set, λ is counting measure, and m(x) > 0 for all x ∈ X . Then choosing U = X , conditions (1) and (2) of The- orem 3.4.1 are trivially satisfied, regardless of the distribution of X1 , X2 , . . . . (Note that when λ is counting measure, pθ (x) and m(x) are p.m.f.s on X .) It is often easy to check that m(x) > 0 by using the fact that this is true when- ever {θ ∈ Θ : pθ (x) > 0} has nonzero probability under π. This case covers, for instance, Multinomials (including Binomials), and the population genetics model from Section 3.1.1. We should mention a subtle point here: when X is finite, mixture identifia- bility might only hold up to a certain maximum number of components (e.g., Teicher (1963), Proposition 4, showed this for Binomials), making consistency impossible in general — however, consistency might still be possible within that identifiable range. Regardless, our result shows that PYMs are not consistent anyway. iid Now, suppose P is a probability measure on X , and X1 , X2 , . . . ∼ P . Let us abuse notation and write P (x) = P ({x}) and λ(x) = λ({x}) for x ∈ X . (ii) One or more point masses in common. If there exists x0 ∈ X such that P (x0 ) > 0, λ(x0 ) > 0, and m(x0 ) > 0, then it is easy to verify that 41 conditions (1) and (2) are satisfied with U = {x0 }. (Note that λ(x0 ) > 0 implies pθ (x0 ) ≤ 1/λ(x0 ) for any θ ∈ Θ.) (iii) Discrete families. Case (ii) essentially covers all discrete families — e.g., Poisson, Geometric, Negative Binomial, or any power-series distribution (see Sapatinas (1995) for mixture identifiability of these) — provided that the data is i.i.d.. For, suppose X is a countable set and λ is counting measure. By case (ii), the theorem applies if there is any x0 ∈ X such that m(x0 ) > 0 and P (x0 ) > 0. If this is not so, the model is extremely misspecified, since then the model distribution and the data distribution are mutually singular. (iv) Continuous densities bounded on some non-null compact set. Suppose there exists c ∈ (0, ∞) and U ⊂ X compact such that (a) P (U ) > 0, (b) x 7→ pθ (x) is continuous on U for all θ ∈ Θ, and (c) pθ (x) ∈ (0, c] for all x ∈ U , θ ∈ Θ. Then condition (1) is satisfied due to item (a), and condition (2) follows eas- ily from (b) and (c) since m(x) is continuous (by the dominated convergence theorem) and positive on the compact set U , so inf x∈U m(x) > 0. This case covers, for example, the following families (with any P ): (a) Exponential(θ), X = (0, ∞), (b) Gamma(a, b), X = (0, ∞), with variance a/b2 bounded away from zero, (c) Normal(µ, Σ), X = Rd , (multivariate Gaussian) with det(Σ) bounded away from zero, and (d) many location–scale families with scale bounded away from zero (for in- stance, Laplace(µ, σ) or Cauchy(µ, σ), with σ ≥ ε > 0). 42 The examples listed in item (iv) are indicative of a deficiency in Theorem 3.4.1: condition (2) is not satisfied in some important cases, such as multivariate Gaussians with unrestricted covariance. Showing that Condition 3.3.2 still holds, for many exponential families at least, is the objective of the remainder of the chapter. 3.5 Exponential families and conjugate priors 3.5.1 Exponential families In this section, we make the usual definitions for exponential families and state the regularity conditions to be assumed. Consider an exponential family of the following form. Fix a sigma-finite Borel measure λ on X ⊂ Rd such that λ(X ) 6= 0, let s : X → Rk be Borel measurable, and for θ ∈ Θ ⊂ Rk , define a density pθ with respect to λ by setting pθ (x) = exp(θT s(x) − κ(θ)) where Z κ(θ) = log exp(θT s(x)) dλ(x). X Let Pθ be the probability measure on X corresponding to pθ , that is, Pθ (E) = R p (x) dλ(x) for E ⊂ X measurable. Any exponential family on Rd can be written E θ in the form above by reparametrizing if necessary, and choosing λ appropriately. We will assume the following (very mild) regularity conditions. Conditions 3.5.1. Assume the family {Pθ : θ ∈ Θ} is: (1) full, that is, Θ = {θ ∈ Rk : κ(θ) < ∞}, 43 (2) nonempty, that is, Θ 6= ∅, (3) regular, that is, Θ is an open subset of Rk , and (4) identifiable, that is, if θ 6= θ0 then Pθ 6= Pθ0 . Most commonly-used exponential families satisfy Conditions 3.5.1, including mul- tivariate Gaussian, Gamma, Poisson, Exponential, Geometric, and others. (A no- table exception is the Inverse Gaussian, for which Θ is not open.) Let M denote the moment space, that is, M = {Eθ s(X) : θ ∈ Θ} where Eθ denotes expectation under Pθ . Finiteness of these expectations is guaran- teed, thus M ⊂ Rk ; see Section 3.10 for this and other well-known properties that we will use. 3.5.2 Conjugate priors Given an exponential family {Pθ } as above, let n o Ξ = (ξ, ν) : ξ ∈ Rk , ν > 0 s.t. ξ/ν ∈ M , and consider the family {πξ,ν : (ξ, ν) ∈ Ξ} where  πξ,ν (θ) = exp ξ T θ − νκ(θ) − ψ(ξ, ν) I(θ ∈ Θ) is a density with respect to Lebesgue measure on Rk . Here, Z  ψ(ξ, ν) = log exp ξ T θ − νκ(θ) dθ. Θ 44 In Section 3.10, we note a few basic properties of this family — in particular, it is a conjugate prior for {Pθ }. Definition 3.5.2. We will say that an exponential family with conjugate prior is well-behaved if it takes the form above, satisfies Conditions 3.5.1, and has (ξ, ν) ∈ Ξ. 3.6 Application to exponential families In this section, we apply Theorem 3.3.4 to prove that in many cases, a PYM model using a well-behaved exponential family with conjugate prior will exhibit inconsis- tency for the number of components. Conditions 3.6.1. Consider an exponential family with sufficient statistics function s : X → Rk and moment space M. Given a probability measure P on X , let X ∼ P and assume: (1) E|s(X)| < ∞, (2) P(s(X) ∈ M) = 1, and (3) P(s(X) ∈ L) = 0 for any hyperplane L that does not intersect M. Throughout, we use | · | to denote the Euclidean norm. Here, a hyperplane refers to a set L = {x ∈ Rk : xT y = b} for some y ∈ Rk r {0}, b ∈ R. In Theorem 3.6.2 below, it is assumed that the data comes from a distribution P satisfying Conditions 3.6.1. In Proposition 3.6.3, we give some simple conditions under which, if P is a finite mixture from the exponential family under consideration, then Conditions 3.6.1 hold. 45 The following theorem follows almost immediately from Lemma 3.8.4, the proof of which will occupy most of the remainder of the chapter. Theorem 3.6.2. Consider a well-behaved exponential family with conjugate prior (as in Definition 3.5.2), along with the resulting collection of single-cluster marginals m(·). Let P be a probability measure on X satisfying Conditions 3.6.1 (for the s and iid M from the exponential family under consideration), and let X1 , X2 , . . . ∼ P . Then Condition 3.3.2 holds for any t ∈ {1, 2, . . . }. Proof. Let t ∈ {1, 2, . . . } and choose c according to Lemma 3.8.4 with β = 1/t. We will show that for any n > t, if the event of Lemma 3.8.4 holds, then ϕt (X1:n , c) ≥ 1/(2t). Since with probability 1, this event holds for all n sufficiently large, it will follow that with probability 1, lim inf n ϕt (X1:n , c) ≥ 1/(2t) > 0. So, let n > t and x1 , . . . , xn ∈ X , and assume the event of Lemma 3.8.4 holds. Let A ∈ At (n). There is at least one part A` such that |A` | ≥ n/t = βn. Then, by assumption there exists RA ⊂ A` such that |RA | ≥ 21 |A` | and for any j ∈ RA , m(xA` ) ≤ c m(xA` rj ) m(xj ). Thus, RA ⊂ SA (x1:n , c), hence |SA (x1:n , c)| ≥ |RA | ≥ 1 2 |A` | ≥ n/(2t). Since A ∈ At (n) was arbitrary, ϕt (x1:n , c) ≥ 1/(2t). This theorem implies inconsistency in several important cases. In particular, it can be verified that each of the following is well-behaved (when put in canonical form and given the conjugate prior in Section 3.5.2) and, using Proposition 3.6.3 below, that if P is a finite mixture from the same family then P satisfies Conditions 3.6.1: (a) Normal(µ, Σ) (multivariate Gaussian), (b) Exponential(θ), 46 (c) Gamma(a, b), (d) Log-Normal(µ, σ 2 ), and (e) Weibull(a, b) with fixed shape a > 0. Combined with the cases covered by Theorem 3.4.1, these results are fairly compre- hensive. Proposition 3.6.3. Consider an exponential family {Pθ : θ ∈ Θ} satisfying Con- Pt ditions 3.5.1. If X ∼ P = i=1 πi Pθ(i) for some θ(1), . . . , θ(t) ∈ Θ and some π1 , . . . , πt ≥ 0 such that ti=1 πi = 1, then P (1) E|s(X)| < ∞, and (2) P(s(X) ∈ M) = 1. If, further, the exponential family is continuous (that is, the underlying measure λ is absolutely continuous with respect to Lebesgue measure on X ), X ⊂ Rd is open and connected, and the sufficient statistics function s : X → Rk is real analytic (that is, each coordinate function s1 , . . . , sk is real analytic), then (3) P(s(X) ∈ L) = 0 for any hyperplane L ⊂ Rk . Proof. See Section 3.9. Sometimes, Condition 3.6.1(3) will be satisfied even when Proposition 3.6.3 is not applicable. In any particular case, it may be a simple matter to check this condition by using the characterization of M as the interior of the closed convex hull of support(λs−1 ) (see Proposition 3.10.1(8) in Section 3.10). 47 3.7 Marginal inequalities Consider a well-behaved exponential family with conjugate prior (as in Definition 3.5.2). In this section, we use some simple bounds on the Laplace approximation (see Section 3.11) to prove certain inequalities involving the marginal density (from Equation 3.2.5), n Z Y  m(x1:n ) = pθ (xj ) πξ,ν (θ) dθ Θ j=1 of x1:n = (x1 , . . . , xn ), where xj ∈ X . Of course, it is commonplace to apply the Laplace approximation to m(X1:n ) when X1 , . . . , Xn are i.i.d. random variables. In contrast, our application of it is considerably more subtle. For our purposes, it is necessary to show that the approximation is good not only in the i.i.d. case, but in fact whenever the sufficient statistics are not too extreme. We make extensive use of the exponential family properties in Section 3.10, often without mention. We use f 0 to denote the gradient and f 00 to denote the Hessian of a (sufficiently smooth) function f : Rk → R. For µ ∈ M, define fµ (θ) = θT µ − κ(θ),  L(µ) = sup θT µ − κ(θ) , θ∈Θ  θµ = argmax θT µ − κ(θ) , θ∈Θ and note that θµ = κ0−1 (µ) (Proposition 3.10.1). L is known as the Legendre transform of κ. Note that L(µ) = fµ (θµ ), and L is C ∞ smooth on M (since L(µ) = θµT µ − κ(θµ ), θµ = κ0−1 (µ), and both κ and κ0−1 are C ∞ smooth). De- 48 fine Pn ξ+ s(xj ) j=1 µx1:n = (3.7.1) ν+n (cf. Equation 3.10.1), and given x1:n such that µx1:n ∈ M, define e 1:n ) = (ν + n)−k/2 exp (ν + n) L(µx1:n ) ,  m(x where k is the dimension of the sufficient statistics function s : X → Rk . The first of the two results of this section provides uniform bounds on m(x1:n )/m(x e 1:n ). Here, m(x e 1:n ) is only intended to approximate m(x1:n ) up to a multiplicative constant — a better approximation could always be obtained via the usual asymptotic form of the Laplace approximation. Proposition 3.7.1. Consider a well-behaved exponential family with conjugate prior. For any U ⊂ M compact, there exist C1 , C2 ∈ (0, ∞) such that for any n ∈ {1, 2, . . . } and any x1 , . . . , xn ∈ X satisfying µx1:n ∈ U , we have m(x1:n ) C1 ≤ ≤ C2 . m(x e 1:n ) Proof. Assume U 6= ∅, since otherwise the result is trivial. Let V = κ0−1 (U ) = {θµ : µ ∈ U }. It is straightforward to show that there exists ε ∈ (0, 1) such that Vε ⊂ Θ where Vε = {θ ∈ Rk : d(θ, V ) ≤ ε}. (Here, d(θ, V ) = inf θ0 ∈V |θ − θ0 |.) Note that Vε is compact, since κ0−1 is continuous. 49 Given a symmetric matrix A, define λ∗ (A) and λ∗ (A) to be the minimal and maximal eigenvalues, respectively, and recall that λ∗ , λ∗ are continuous functions of the entries of A. Letting α = min λ∗ (κ00 (θ)) and β = max λ∗ (κ00 (θ)), θ∈Vε θ∈Vε we have 0 < α ≤ β < ∞ since Vε is compact and λ∗ (κ00 (·)), λ∗ (κ00 (·)) are continuous and positive on Θ. Letting Z −fµ (θµ ) γ = sup e exp(fµ (θ))dθ = sup e−L(µ) eψ(µ,1) µ∈U Θ µ∈U we have 0 < γ < ∞ since U is compact, and both L (as noted above) and ψ(µ, 1) (by Proposition 3.10.2) are continuous on M. Define h(µ, θ) = fµ (θµ ) − fµ (θ) = L(µ) − θT µ + κ(θ) for µ ∈ M, θ ∈ Θ. For any µ ∈ M, we have that h(µ, θ) > 0 whenever θ ∈ Θ r {θµ }, and that h(µ, θ) is strictly convex in θ. Letting Bε (θµ ) = {θ ∈ Rk : |θ − θµ | ≤ ε}, it follows that δ := inf inf h(µ, θ) = inf inf h(µ, θµ + εu) µ∈U θ∈ΘrBε (θµ ) µ∈U u∈Rk :|u|=1 is positive, as the minimum of a positive continuous function on a compact set. Now, applying the Laplace approximation bounds in Corollary 3.11.2 with α, β, γ, δ, ε as just defined, we obtain c1 , c2 ∈ (0, ∞) such that for any µ ∈ U we have (taking E = Θ, f = −fµ , x0 = θµ , A = αIk×k , B = βIk×k ) R Θ exp(tfµ (θ))dθ c1 ≤ −k/2 ≤ c2 t exp(tfµ (θµ )) 50 for any t ≥ 1. We prove the result with Ci = ci e−ψ(ξ,ν) for i = 1, 2. Let n ∈ {1, 2, . . . } and x1 , . . . , xn ∈ X such that µx1:n ∈ U . Choose t = ν + n. By integrating Equation 3.10.1, we have Z −ψ(ξ,ν)  m(x1:n ) = e exp tfµx1:n (θ) dθ, Θ and meanwhile, e 1:n ) = t−k/2 exp tfµx1:n (θµx1:n ) .  m(x Thus, combining the preceding three displayed equations, m(x1:n ) 0 < C1 = c1 e−ψ(ξ,ν) ≤ ≤ c2 e−ψ(ξ,ν) = C2 < ∞. m(x e 1:n ) The second result of this section is an inequality involving a product of marginals. Proposition 3.7.2 (Splitting inequality). Consider a well-behaved exponential fam- ily with conjugate prior. For any U ⊂ M compact there exists C ∈ (0, ∞) such that we have the following: For any n ∈ {1, 2, . . . }, if A ⊂ {1, . . . , n} and B = {1, . . . , n} r A are nonempty, 1 P and x1 , . . . , xn ∈ X satisfy |A| j∈A s(xj ) ∈ U and µxB ∈ U , then m(x1:n )  ab k/2 ≤C m(xA )m(xB ) ν+n where a = ν + |A| and b = ν + |B|. Proof. Let U 0 be the convex hull of U ∪ {ξ/ν}. Then U 0 is compact (as the convex hull of a compact set in Rk ) and U 0 ⊂ M (since U ∪ {ξ/ν} ⊂ M and M is convex). 51 We show that the result holds with C = C2 exp(C0 )/C12 , where C1 , C2 ∈ (0, ∞) are obtained by applying Proposition 3.7.1 to U 0 , and C0 = ν sup |(ξ/ν − y)T L0 (y)| + ν sup |L(y)| < ∞. (3.7.2) y∈U 0 y∈U 0 Since L is convex (being a Legendre transform) and smooth, then for any y, z ∈ M we have 1 L(y + ρ(z − y)) − L(y) = (z − y)T L0 (y)  inf ρ∈(0,1) ρ (by e.g. Rockafellar (1970) 23.1) and therefore for any ρ ∈ (0, 1), L(y) ≤ L((1 − ρ)y + ρz) − ρ(z − y)T L0 (y). (3.7.3) Choosing y = µx1:n , z = ξ/ν, and ρ = ν/(n + 2ν), we have Pn 2ξ + j=1 s(xj ) aµxA + bµxB (1 − ρ)y + ρz = = . (3.7.4) 2ν + n a+b Note that µxA , µxB , µx1:n ∈ U 0 , by taking various convex combinations of ξ/ν, 1 0 P |A| j∈A s(xj ), µxB ∈ U . Thus, (ν + n)L(µx1:n ) = (a + b)L(y) − νL(y) (a) ≤ (a + b)L((1 − ρ)y + ρz) − (a + b)ρ(z − y)T L0 (y) − νL(y) (b)  aµ + bµ  xA xB ≤ (a + b)L + C0 a+b (c) ≤ aL(µxA ) + bL(µxB ) + C0 , where (a) is by Equation 3.7.3, (b) is by Equations 3.7.2 and 3.7.4, and (c) is by the 52 convexity of L. Hence, (ν + n)k/2 m(x e 1:n ) ≤ (ab)k/2 m(x e A )m(x e B ) exp(C0 ), so by our choice of C1 and C2 , m(x1:n ) C2 m(x e 1:n ) C2 exp(C0 )  ab k/2 ≤ 2 ≤ . m(xA )m(xB ) C1 m(x e A )m(x e B) C12 n+ν 3.8 Marginal inequality for subsets of the data In this section, we prove Lemma 3.8.4, the key lemma used in the proof of Theorem 3.6.2. First, we need a few supporting results. Given y1 , . . . , yn ∈ R` (for some ` > 0), β ∈ (0, 1], and U ⊂ R` , define Y  1 X  Iβ (y1:n , U ) = I yj ∈ U , |A| j∈A A⊂{1,...,n}: |A|≥βn where as usual, I(E) is 1 if E is true, and 0 otherwise. Lemma 3.8.1 (Capture lemma). Let V ⊂ Rk be open and convex. Let Q be a probability measure on Rk such that: (1) E|Y | < ∞ when Y ∼ Q, (2) Q(V ) = 1, and (3) Q(L) = 0 for any hyperplane L that does not intersect V . iid If Y1 , Y2 , . . . ∼ Q, then for any β ∈ (0, 1] there exists U ⊂ V compact such that a.s. Iβ (Y1:n , U ) −−→ 1 as n → ∞. 53 Proof. The proof is rather long, but not terribly difficult; see Section 3.12. Proposition 3.8.2. Let Z1 , Z2 , . . . ∈ Rk be i.i.d.. If β ∈ (0, 1] and U ⊂ Rk such a.s. that P(Zj 6∈ U ) < β/2, then Iβ (Y1:n , [ 21 , 1]) −−→ 1 as n → ∞, where Yj = I(Zj ∈ U ). 1 Pn a.s. Proof. By the law of large numbers, n I(Zj 6∈ U ) −−→ P(Zj 6∈ U ) < β/2. j=1 Hence, with probability 1, for all n sufficiently large, n1 nj=1 I(Zj 6∈ U ) ≤ β/2 holds. P When it holds, we have that for any A ⊂ {1, . . . , n} such that |A| ≥ βn, n 1 X 1 X 1 X I(Zj ∈ U ) = 1 − I(Zj 6∈ U ) ≥ 1 − I(Zj 6∈ U ) ≥ 1/2, |A| j∈A |A| j∈A βn j=1 a.s. i.e. when it holds, we have Iβ (Y1:n , [ 12 , 1]) = 1. Hence, Iβ (Y1:n , [ 21 , 1]) −−→ 1. In the following, µx = (ξ + s(x))/(ν + 1), as in Equation 3.7.1. Proposition 3.8.3. Consider a well-behaved exponential family with conjugate prior. Let P be a probability measure on X such that P(s(X) ∈ M) = 1 when iid X ∼ P . Let X1 , X2 , . . . ∼ P . Then for any β ∈ (0, 1] there exists U ⊂ M compact a.s. such that Iβ (Y1:n , [ 12 , 1]) −−→ 1 as n → ∞, where Yj = I(µXj ∈ U ). Proof. Since M is open and convex, then for any y ∈ M, z ∈ M, and ρ ∈ (0, 1), we have ρy + (1 − ρ)z ∈ M (by e.g. Rockafellar (1970) 6.1). Taking z = ξ/ν and ρ = 1/(ν + 1), this implies that the set U0 = {(ξ + y)/(ν + 1) : y ∈ M} is contained in M. Note that U0 is closed and P(µX ∈ U0 ) = P(s(X) ∈ M) = 1. Let β ∈ (0, 1], and choose r ∈ (0, ∞) such that P(|µX | > r) < β/2. Letting U = {y ∈ U0 : |y| ≤ r}, we have that U ⊂ M, and U is compact. Further, P(µX 6∈ U ) < β/2, so by applying a.s. Proposition 3.8.2 with Zj = µXj , we have Iβ (Y1:n , [ 21 , 1]) −−→ 1. 54 Lemma 3.8.4. Consider a well-behaved exponential family with conjugate prior, and the resulting collection of single-cluster marginals m(·). Let P be a probability measure on X satisfying Conditions 3.6.1 (for the s and M from the exponential iid family under consideration), and let X1 , X2 , . . . ∼ P . Then for any β ∈ (0, 1] there exists c ∈ (0, ∞) such that with probability 1, for all n sufficiently large, the following event holds: for every subset J ⊂ {1, . . . , n} such that |J| ≥ βn, there exists K ⊂ J such that |K| ≥ 12 |J| and for any j ∈ K, m(XJ ) ≤ c m(XJrj ) m(Xj ). Proof. Let β ∈ (0, 1]. Since M is open and convex, and Conditions 3.6.1 hold by assumption, then by Lemma 3.8.1 (with V = M) there exists U1 ⊂ M compact a.s. such that Iβ/2 (s(X1:n ), U1 ) −−→ 1 as n → ∞, where s(X1:n ) = (s(X1 ), . . . , s(Xn )). By a.s. Proposition 3.8.3 above, there exists U2 ⊂ M compact such that Iβ (Y1:n , [ 21 , 1]) −−→ 1 as n → ∞, where Yj = I(µXj ∈ U2 ). Hence, a.s. Iβ/2 (s(X1:n ), U1 ) Iβ (Y1:n , [ 21 , 1]) −−−→ 1. n→∞ Choose C ∈ (0, ∞) according to Proposition 3.7.2 applied to U := U1 ∪ U2 . We will prove the result with c = (ν +1)k/2 C. (Recall that k is the dimension of s : X → Rk .) Let n large enough that βn ≥ 2, and suppose that Iβ/2 (s(X1:n ), U1 ) = 1 and Iβ (Y1:n , [ 21 , 1]) = 1. Let J ⊂ {1, . . . , n} such that |J| ≥ βn. Then for any j ∈ J, 1 X s(Xi ) ∈ U1 ⊂ U |J r j| i∈Jrj since Iβ/2 (s(X1:n ), U1 ) = 1 and |J r j| ≥ |J|/2 ≥ (β/2)n. Hence, for any j ∈ K, 55 where K = {j ∈ J : µXj ∈ U }, we have m(XJ )  (ν + |J| − 1)(ν + 1) k/2 ≤C ≤ C (ν + 1)k/2 = c m(XJrj ) m(Xj ) ν + |J| by our choice of C above, and |K| 1 X 1 X ≥ I(µXj ∈ U2 ) = Yj ≥ 1/2 |J| |J| j∈J |J| j∈J since Iβ (Y1:n , [ 12 , 1]) = 1 and |J| ≥ βn. 3.9 Miscellaneous proofs Proof of Proposition 3.3.3. There are two cases: (A) σ ∈ [0, 1) and ϑ > −σ, or (B) σ < 0 and ϑ = N |σ|. In either case, σ < 1, so wn (a) 1−σ+a−2 1−σ = ≤ +1 awn (a − 1)wn (1) a 2 whenever n ≥ 2 and a ∈ {2, . . . , n}, and hence lim supn cwn < ∞. For any n > t ≥ 1, in case (A) we have vn (t) t+1 = , vn (t + 1) ϑ + tσ and the same holds in case (B) if also t < N . Meanwhile, whenever N < t < n in case (B), vn (t)/vn (t + 1) = 0/0 = 0 by convention. Therefore, lim supn cvn (t) < ∞ in either case, for any t ∈ {1, 2, . . . } except t = N in case (B). 56 Proof of Theorem 3.3.4. Let t ∈ {1, 2, . . . }, and assume Conditions 3.3.1 and 3.3.2 hold. Let x1 , x2 , . . . ∈ X , and suppose supc∈[0,∞) lim inf n ϕt (x1:n , c) > 0 (which occurs with probability 1). We show that this implies lim supn p(Tn = t | x1:n ) < 1, proving the theorem. Let α ∈ (0, ∞) such that lim supn cwn < α and lim supn cvn (t) < α. Choose c ∈ [0, ∞) and ε ∈ (0, 1) such that lim inf n ϕt (x1:n , c) > ε. Choose N > 2t/ε large enough that for any n > N we have cwn < α, cvn (t) < α, and ϕt (x1:n , c) > ε. Then by Lemma 3.3.5, for any n > N , Ct (x1:n , c) 2tcα2 /ε p(Tn = t | x1:n ) ≤ ≤ , 1 + Ct (x1:n , c) 1 + 2tcα2 /ε since ϕt (x1:n , c) − t/n > ε − ε/2 = ε/2 (and y 7→ y/(1 + y) is monotone increasing on [0, ∞)). Since this upper bound does not depend on n (and is less than 1), then lim supn p(Tn = t | x1:n ) < 1. Proof of Proposition 3.6.3. (1) For any θ ∈ Θ and any j ∈ {1, . . . , k}, ∂2 Z Z 2 sj (x) pθ (x) dλ(x) = exp(−κ(θ)) 2 exp(θT s(x)) dλ(x) < ∞ X ∂θj X P (Hoffmann-Jørgensen, 1994, 8.36.1). Since P has density f = πi pθ(i) with respect to λ, then Z t X Z 2 2 Esj (X) = sj (x) f (x) dλ(x) = πi sj (x)2 pθ(i) (x) dλ(x) < ∞, X i=1 X and hence (E|s(X)|)2 ≤ E|s(X)|2 = Es1 (X)2 + · · · + Esk (X)2 < ∞. 57 (2) Note that SP (s) ⊂ Sλ (s) (in fact, they are equal since Pθ and λ are mutually absolutely continuous for any θ ∈ Θ), and therefore SP (s) ⊂ Sλ (s) ⊂ Cλ (s) = M by Proposition 3.10.1(8). Hence, P(s(X) ∈ M) ≥ P(s(X) ∈ SP (s)) = P s−1 (support(P s−1 )) = 1. (3) Suppose λ is absolutely continuous with respect to Lebesgue measure, X is open and connected, and s is real analytic. Let L ⊂ Rk be a hyperplane, and write L = {z ∈ Rk : z T y = b} where y ∈ Rk r {0}, b ∈ R. Define g : X → R by g(x) = s(x)T y − b. Then g is real analytic on X , since a finite sum of real analytic functions is real analytic. Since X is connected, it follows that either g is identically zero, or the set V = {x ∈ X : g(x) = 0} has Lebesgue measure zero (Krantz, 1992). Now, g cannot be identically zero, since for any θ ∈ Θ, letting Z ∼ Pθ , we have 0 < y T κ00 (θ)y = y T (Cov s(Z))y = Var(y T s(Z)) = Var g(Z) by Proposition 3.10.1(2) and (3). Consequently, V must have Lebesgue measure zero. Hence, P (V ) = 0, since P is absolutely continuous with respect to λ, and thus, with respect to Lebesgue measure. Therefore, P(s(X) ∈ L) = P(g(X) = 0) = P (V ) = 0. 58 3.10 Exponential family properties We note some well-known properties of exponential families satisfying Conditions 3.5.1. For a general reference on this material, see Hoffmann-Jørgensen (1994). Let Sλ (s) = support(λs−1 ), that is, Sλ (s) = z ∈ Rk : λ(s−1 (U )) 6= 0 for every neighborhood U of z .  Let Cλ (s) be the closed convex hull of Sλ (s) (that is, the intersection of all closed convex sets containing it). Given U ⊂ Rk , let U ◦ denote its interior. Given a (sufficiently smooth) function f : Rk → R, we use f 0 to denote its gradient, that is, ∂f ∂2f f 0 (x)i = ∂xi (x), and f 00 (x) to denote its Hessian matrix, that is, f 00 (x)ij = ∂xi ∂xj (x). Proposition 3.10.1. If Conditions 3.5.1 are satisfied, then: (1) κ is C ∞ smooth and strictly convex on Θ, (2) κ0 (θ) = Es(X) and κ00 (θ) = Cov s(X) when θ ∈ Θ and X ∼ Pθ , (3) κ00 (θ) is symmetric positive definite for all θ ∈ Θ, (4) κ0 : Θ → M is a C ∞ smooth bijection, (5) κ0−1 : M → Θ is C ∞ smooth, (6) Θ is open and convex, (7) M is open and convex, (8) M = Cλ (s)◦ and M = Cλ (s), and (9) κ0−1 (µ) = argmaxθ∈Θ (θT µ−κ(θ)) for all µ ∈ M. The maximizing θ ∈ Θ always exists and is unique. 59 Proof. These properties are all well-known. Let us abbreviate Hoffmann-Jørgensen (1994) as HJ. For (1), see HJ 8.36(1) and HJ 12.7.5. For (6),(2),(3), and (4), see HJ 8.36, 8.36.2-3, 12.7(2), and 12.7.11, respectively. Item (5) and openness in (7) follow, using the inverse function theorem (Knapp, 2005, 3.21). Item (8) and convexity in (7) follow, using HJ 8.36.15 and Rockafellar (1970) 6.2-3. Item (9) follows from HJ 8.36.15 and item (4). Given an exponential family with conjugate prior as in Section 3.5.2, the joint density of x1 , . . . , xn ∈ X and θ ∈ Rk is pθ (x1 ) · · · pθ (xn )πξ,ν (θ) (3.10.1)   = exp (ν + n) θT µx1:n − κ(θ) exp(−ψ(ξ, ν)) I(θ ∈ Θ) Pn where µx1:n = (ξ + j=1 s(xj ))/(ν +n). The marginal density, defined as in Equation 3.2.5, is  P   m(x1:n ) = exp ψ ξ + s(xj ), ν + n − ψ(ξ, ν) (3.10.2) when this quantity is well-defined. Proposition 3.10.2. If Conditions 3.5.1 are satisfied, then: (1) ψ(ξ, ν) is finite and C ∞ smooth on Ξ, P (2) if s(x1 ), . . . , s(xn ) ∈ Sλ (s) and (ξ, ν) ∈ Ξ, then (ξ + s(xj ), ν + n) ∈ Ξ, (3) {πξ,ν : (ξ, ν) ∈ Ξ} is a conjugate family for {pθ : θ ∈ Θ}, and (4) if s : X → Rk is continuous, (ξ, ν) ∈ Ξ, and λ(U ) 6= 0 for any nonempty U ⊂ X that is open in X , then m(x1:n ) < ∞ for any x1 , . . . , xn ∈ X . 60 Proof. (1) For finiteness, see Diaconis and Ylvisaker (1979), Theorem 1. Smoothness holds for the same reason that κ is smooth (Hoffmann-Jørgensen, 1994, 8.36(1)). (Note that Ξ is open in Rk+1 , since M is open in Rk .) 1 P (2) Since Cλ (s) is convex, n s(xj ) ∈ Cλ (s). Since Cλ (s) = M and M is open P and convex (3.10.1(7) and (8)), then (ξ + s(xj ))/(ν + n) ∈ M, as a (strict) convex combination of n1 P s(xj ) ∈ M and ξ/ν ∈ M (Rockafellar, 1970, 6.1). iid (3) Let (ξ, ν) ∈ Ξ, θ ∈ Θ. If X1 , . . . , Xn ∼ Pθ then s(X1 ), . . . , s(Xn ) ∈ Sλ (s) P almost surely, and thus (ξ + s(Xj ), ν + n) ∈ Ξ (a.s.) by (2). By Equations 3.10.1 and 3.10.2, the posterior is πξ+P s(Xj ), ν+n . (4) The assumptions imply {s(x) : x ∈ X } ⊂ Sλ (s), and therefore, for any P x1 , . . . , xn ∈ X , we have (ξ + s(xj ), ν + n) ∈ Ξ by (2). Thus, by (1) and Equation 3.10.2, m(x1:n ) < ∞.  It is worth mentioning that while Ξ ⊂ (ξ, ν) ∈ Rk+1 : ψ(ξ, ν) < ∞ , it may be a strict subset — often, Ξ is not quite the full set of parameters on which πξ,ν can be defined. 3.11 Bounds on the Laplace approximation Our proof uses the following simple bounds on the Laplace approximation. These bounds are not fundamentally new, but the precise formulation we require does not seem to appear in the literature, so we have included it for the reader’s convenience. Lemma 3.11.1 is simply a multivariate version of the bounds given by De Bruijn (1970), and Corollary 3.11.2 is a straightforward consequence, putting the lemma in 61 a form most convenient for our purposes. Given symmetric matrices A and B, let us write A  B to mean that B − A is positive semidefinite. Given A ∈ Rk×k symmetric positive definite and ε, t ∈ (0, ∞), define √ C(t, ε, A) = P(|A−1/2 Z| ≤ ε t) where Z ∼ Normal(0, Ik×k ). Note that C(t, ε, A) → 1 as t → ∞. Let Bε (x0 ) = {x ∈ Rk : |x − x0 | ≤ ε} denote the closed ball of radius ε > 0 at x0 ∈ Rk . Lemma 3.11.1. Let E ⊂ Rk be open. Let f : E → R be C 2 smooth with f 0 (x0 ) = 0 for some x0 ∈ E. Define Z g(t) = exp(−tf (x)) dx E for t ∈ (0, ∞). Suppose ε ∈ (0, ∞) such that Bε (x0 ) ⊂ E, 0 < δ ≤ inf{f (x) − f (x0 ) : x ∈ E r Bε (x0 )}, and A, B are symmetric positive definite matrices such that A  f 00 (x)  B for all x ∈ Bε (x0 ). Then for any 0 < s ≤ t we have C(t, ε, B) g(t) C(t, ε, A)  t k/2 −(t−s)δ sf (x0 ) ≤ ≤ + e e g(s) |B|1/2 (2π/t)k/2 e−tf (x0 ) |A|1/2 2π where |A| = | det A|. Remark. In particular, these assumptions imply f is strictly convex on Bε (x0 ) with unique global minimum at x0 . Note that the upper bound is trivial unless g(s) < ∞. Proof. By Taylor’s theorem, for any x ∈ Bε (x0 ) there exists zx on the line between x0 and x such that, letting y = x − x0 , f (x) = f (x0 ) + y T f 0 (x0 ) + 21 y T f 00 (zx )y = f (x0 ) + 12 y T f 00 (zx )y. 62 Since zx ∈ Bε (x0 ), and thus A  f 00 (zx )  B, 1 T 2 y Ay ≤ f (x) − f (x0 ) ≤ 21 y T By. Hence, Z Z tf (x0 ) e exp(−tf (x)) dx ≤ exp(− 12 (x − x0 )T (tA)(x − x0 )) dx Bε (x0 ) Bε (x0 ) = (2π)k/2 |(tA)−1 |1/2 P |(tA)−1/2 Z| ≤ ε .  Along with a similar argument for the lower bound, this implies  2π k/2 C(t, ε, B) Z  2π k/2 C(t, ε, A) tf (x0 ) ≤e exp(−tf (x)) dx ≤ . t |B|1/2 Bε (x0 ) t |A|1/2 Considering the rest of the integral, outside of Bε (x0 ), we have Z  0≤ exp(−tf (x)) dx ≤ exp − (t − s)(f (x0 ) + δ) g(s). ErBε (x0 ) Combining the preceding four inequalities yields the result. The following corollary tailors the lemma to our purposes. Given a symmetric positive definite matrix A ∈ Rk×k , let λ∗ (A) and λ∗ (A) be the minimal and maximal eigenvalues, respectively. By diagonalizing A, it is easy to check that λ∗ (A)Ik×k  A  λ∗ (A)Ik×k and λ∗ (A)k ≤ |A| ≤ λ∗ (A)k . Corollary 3.11.2. For any α, β, γ, δ, ε ∈ (0, ∞) there exist c1 = c1 (β, ε) ∈ (0, ∞) and c2 = c2 (α, γ, δ) ∈ (0, ∞) such that if E, f, x0 , A, B satisfy all the conditions of Lemma 3.11.1 (for this choice of δ, ε) and additionally, α ≤ λ∗ (A), β ≥ λ∗ (B), and γ ≥ ef (x0 ) g(1), then R E exp(−tf (x)) dx c1 ≤ ≤ c2 t−k/2 exp(−tf (x0 )) 63 for all t ≥ 1. Proof. The first term in the upper bound of the lemma is C(t, ε, A)/|A|1/2 ≤ 1/αk/2 , and with s = 1 the second term is less or equal to (t/2π)k/2 e−(t−1)δ γ, which is bounded above for t ∈ [1, ∞). For the lower bound, a straightforward calculation (using z T Bz ≤ λ∗ (B)z T z ≤ βz T z in the exponent inside the integral) shows that √ C(t, ε, B)/|B|1/2 ≥ P(|Z| ≤ ε β)/β k/2 for t ≥ 1. Although we do not need it (and thus, we omit the proof), the following corollary gives the well-known asymptotic form of the Laplace approximation. (As usual, g(t) ∼ h(t) as t → ∞ means that g(t)/h(t) → 1.) Corollary 3.11.3. Let E ⊂ Rk be open. Let f : E → R be C 2 smooth such that for some x0 ∈ E we have that f 0 (x0 ) = 0, f 00 (x0 ) is positive definite, and f (x) > f (x0 ) for all x ∈ E r {x0 }. Suppose there exists ε > 0 such that Bε (x0 ) ⊂ E and inf{f (x) − f (x0 ) : x ∈ E r Bε (x0 )} is positive, and suppose there is some s > 0 such that E e−sf (x) dx < ∞. Then R Z  2π k/2 exp(−tf (x )) 0 exp(−tf (x)) dx ∼ 00 1/2 E t |f (x0 )| as t → ∞. 3.12 Capture lemma In this section, we prove Lemma 3.8.1, which is restated here for the reader’s conve- nience. 64 The following definitions are standard. Let S denote the unit sphere in Rk , that is, S = {x ∈ Rk : |x| = 1}. We say that H ⊂ Rk is a halfspace if H = {x ∈ Rk : xT u ≺ b}, where ≺ is either < or ≤, for some u ∈ S, b ∈ R. We say that L ⊂ Rk is a hyperplane if L = {x ∈ Rk : xT u = b} for some u ∈ S, b ∈ R. Given U ⊂ Rk , let ∂U denote the boundary of U , that is, ∂U = U r U ◦ . So, for example, if H is a halfspace, then ∂H is a hyperplane. The following notation is also useful: given x ∈ Rk , we call the set Rx = {ax : a > 0} the ray through x. We give the central part of the proof first, postponing some plausible intermediate results for the moment. Lemma 3.12.1 (Capture lemma). Let V ⊂ Rk be open and convex. Let P be a probability measure on Rk such that: (1) E|X| < ∞ when X ∼ P , (2) P (V ) = 1, and (3) P (L) = 0 for any hyperplane L that does not intersect V . iid If X1 , X2 , . . . ∼ P , then for any β ∈ (0, 1] there exists U ⊂ V compact such that a.s. Iβ (X1:n , U ) −−→ 1 as n → ∞. Proof. Without loss of generality, we may assume 0 ∈ V (since otherwise we can translate to make it so, obtain U , and translate back). Let β ∈ (0, 1]. By Proposition 3.12.3 below, for each u ∈ S there is a closed halfspace Hu such that 0 ∈ Hu◦ , Ru a.s. intersects V ∩ ∂Hu , and Iβ (X1:n , Hu ) −−→ 1 as n → ∞. By Proposition 3.12.5 below, there exist u1 , . . . , ur ∈ S (for some r > 0) such that the set U = ri=1 Hui is compact T 65 and U ⊂ V . Finally, r a.s. Y Iβ (X1:n , U ) = Iβ (X1:n , Hui ) −−−→ 1. n→∞ i=1 The main idea of the lemma is exhibited in the following simpler case, which we will use to prove Proposition 3.12.3. Proposition 3.12.2. Let V = (−∞, c), where −∞ < c ≤ ∞. Let P be a probability measure on R such that: (1) E|X| < ∞ when X ∼ P , and (2) P (V ) = 1. iid If X1 , X2 , . . . ∼ P , then for any β ∈ (0, 1] there exists b < c such that a.s. Iβ (X1:n , (−∞, b]) −−→ 1 as n → ∞. Proof. Let β ∈ (0, 1]. By continuity from above, there exists a < c such that P(X > a) < β. If P(X > a) = 0 then the result is trivial, taking b = a. Suppose P(X > a) > 0. Let b such that E(X | X > a) < b < c, which is always possible, by a straightforward argument (using E|X| < ∞ in the c = ∞ case). Let Bn = Bn (X1 , . . . , Xn ) = {i ∈ {1, . . . , n} : Xi > a}. Then n 1 X 1 1X Xi = 1 Xi I(Xi > a) |Bn | i∈B n |Bn | n i=1 n a.s. E(X I(X > a)) −−−→ = E(X | X > a) < b. n→∞ P(X > a) 1 P Now, fix n ∈ {1, 2, . . . }, and suppose 0 < |Bn | < βn and |Bn | i∈Bn Xi < b, noting that with probability 1, this happens for all n sufficiently large. We show that this 66 implies Iβ (X1:n , (−∞, b]) = 1. This will prove the result. Let A ⊂ {1, . . . , n} such that |A| ≥ βn. Let M = {π1 , . . . , π|A| } where π is a permutation of {1, . . . , n} such that Xπ1 ≥ · · · ≥ Xπn (that is, M ⊂ {1, . . . , n} consists of the indices of |A| of the largest entries of (X1 , . . . , Xn )). Then |M | = |A| ≥ βn ≥ |Bn |, and it follows that Bn ⊂ M . Therefore, 1 X 1 X 1 X Xi ≤ Xi ≤ Xi ≤ b, |A| i∈A |M | i∈M |Bn | i∈B n as desired. The first of the two propositions used in Lemma 3.12.1 is the following. Proposition 3.12.3. Let V and P satisfy the conditions of Lemma 3.12.1, and also iid assume 0 ∈ V . If X1 , X2 , . . . ∼ P then for any β ∈ (0, 1] and any u ∈ S there is a closed halfspace H ⊂ Rk such that (1) 0 ∈ H ◦ , (2) Ru intersects V ∩ ∂H, and a.s. (3) Iβ (X1:n , H) −−→ 1 as n → ∞. Proof. Let β ∈ (0, 1] and u ∈ S. Either (a) Ru ⊂ V , or (b) Ru intersects ∂V . (Case (a)) Suppose Ru ⊂ V . Let Yi = XiT u for i = 1, 2, . . . . Then E|Yi | ≤ E|Xi ||u| = E|Xi | < ∞, and thus, by Proposition 3.12.2 (with c = ∞) there exists a.s. b ∈ R such that Iβ (Y1:n , (−∞, b]) −−→ 1. Let us choose this b to be positive, which is always possible since Iβ (Y1:n , (−∞, b]) is nondecreasing as a function of b. Let H = {x ∈ Rk : xT u ≤ b}. Then 0 ∈ H ◦ , since b > 0, and Ru intersects V ∩ ∂H at bu, 67 1 1 P P since Ru ⊂ V and buT u = b. And since |A| i∈A Yi ≤ b if and only if |A| i∈A Xi ∈ H, a.s. we have Iβ (X1:n , H) −−→ 1. (Case (b)) Suppose Ru intersects ∂V at some point z ∈ Rk . Note that z 6= 0 since 0 6∈ Ru . Since V is convex, it has a supporting hyperplane at z, and thus, there exist v ∈ S and c ∈ R such that G = {x ∈ Rk : xT v ≤ c} satisfies V ⊂ G and z ∈ ∂G (Rockafellar, 1970, 11.2). Note that c > 0 and V ∩ ∂G = ∅ since 0 ∈ V and V is open. Letting Yi = XiT v for i = 1, 2, . . . , we have P(Yi ≤ c) = P(XiT v ≤ c) = P(Xi ∈ G) ≥ P(Xi ∈ V ) = P (V ) = 1, and hence, P(Yi ≥ c) = P(Yi = c) = P(XiT v = c) = P(Xi ∈ ∂G) = P (∂G) = 0, by our assumptions on P , since ∂G is a hyperplane that does not intersect V . Con- sequently, P(Yi < c) = 1. Also, as before, E|Yi | < ∞. Thus, by Proposition 3.12.2, a.s. there exists b < c such that Iβ (Y1:n , (−∞, b]) −−→ 1. Since c > 0, we may choose this b to be positive (as before). Letting H = {x ∈ Rk : xT v ≤ b}, we have a.s. Iβ (X1:n , H) −−→ 1. Also, 0 ∈ H ◦ since b > 0. Now, we must show that Ru intersects V ∩ ∂H. First, since z ∈ Ru means z = au for some a > 0, and since z ∈ ∂G means z T v = c > 0, we find that uT v > 0 and z = cu/uT v. Therefore, letting y = bu/uT v, we have y ∈ Ru ∩ V ∩ ∂H, since (i) b/uT v > 0, and thus y ∈ Ru , (ii) y T v = b, and thus y ∈ ∂H, (iii) 0 < b/uT v < c/uT v, and thus y is a (strict) convex combination of 0 ∈ V and 68 z ∈ V , hence y ∈ V (Rockafellar, 1970, 6.1). To prove Proposition 3.12.5, we need the following geometrically intuitive facts. Proposition 3.12.4. Let V ⊂ Rk be open and convex, with 0 ∈ V . Let H be a closed halfspace such that 0 ∈ H ◦ . Let T = {x/|x| : x ∈ V ∩ ∂H}. Then (1) T is open in S, (2) T = {u ∈ S : Ru intersects V ∩ ∂H}, and (3) if x ∈ H, x 6= 0, and x/|x| ∈ T , then x ∈ V . Proof. Write H = {x ∈ Rk : xT v ≤ b}, with v ∈ S, b > 0. Let S+ = {u ∈ S : uT v > 0}. (1) Define f : ∂H → S+ by f (x) = x/|x|, noting that 0 6∈ ∂H. It is easy to see that f is a homeomorphism. Since V is open in Rk , then V ∩ ∂H is open in ∂H. Hence, T = f (V ∩ ∂H) is open in S+ , and since S+ is open in S, then T is also open in S. Items (2) and (3) are easily checked. Proposition 3.12.5. Let V ⊂ Rk be open and convex, with 0 ∈ V . If (Hu : u ∈ S) is a collection of closed halfspaces such that for all u ∈ S, (1) 0 ∈ Hu◦ and (2) Ru intersects V ∩ ∂Hu , Tr then there exist u1 , . . . , ur ∈ S (for some r > 0) such that the set U = i=1 Hui is compact and U ⊂ V . 69 Proof. For u ∈ S, define Tu = {x/|x| : x ∈ V ∩ ∂Hu }. By part (1) of Proposition 3.12.4, Tu is open in S, and by part (2), u ∈ Tu , since Ru intersects V ∩ ∂Hu . Thus, (Tu : u ∈ S) is an open cover of S. Since S is compact, there is a finite subcover: there exist u1 , . . . , ur ∈ S (for some r > 0) such that ri=1 Tui ⊃ S, and in fact, S Sr Tr i=1 Tui = S. Let U = i=1 Hui . Then U is closed and convex (as an intersection of closed, convex sets). Further, U ⊂ V since for any x ∈ U , if x = 0 then x ∈ V by assumption, while if x 6= 0 then x/|x| ∈ Tui for some i ∈ {1, . . . , r} and x ∈ U ⊂ Hui , so x ∈ V by Proposition 3.12.4(3). In order to show that U is compact, we just need to show it is bounded, since we already know it is closed. Suppose not, and let x1 , x2 , . . . ∈ U r {0} such that |xn | → ∞ as n → ∞. Let vn = xn /|xn |. Since S is compact, then (vn ) has a convergent subsequence such that vni → u for some u ∈ S. Then for any a > 0, we have avni ∈ U for all i sufficiently large (since avni is a convex combination of 0 ∈ U and |xni |vni = xni ∈ U whenever |xni | ≥ a). Since avni → au, and U is closed, then au ∈ U . Thus, au ∈ U for all a > 0, i.e. Ru ⊂ U . But u ∈ Tuj for some j ∈ {1, . . . , r}, so Ru intersects ∂Huj (by Proposition 3.12.4(2)), and thus au 6∈ Huj ⊃ U for all a > 0 sufficiently large. This is a contradiction. Therefore, U is bounded. Chapter Four The mixture of finite mixtures model 71 4.1 Introduction In this chapter, we explore the properties of a class of variable-dimension mixture models that we refer to as mixtures of finite mixtures (MFMs), finding that in fact they share many of the same attractive properties as Dirichlet process mixtures (DPMs) — a simple partition distribution, restaurant process, random discrete mea- sure interpretation, stick-breaking representation in some cases, and exchangeability properties — and that many of the inference algorithms developed for DPMs can be directly applied to MFMs as well. Also, as with DPMs, in some cases there are theoretical guarantees on the posterior rate of convergence for density estimation with MFMs. Further, in Chapter 5, we empirically observe that the posterior properties of MFMs and DPMs are remarkably similar in many ways. For instance, density es- timates under the two models are usually nearly indistinguishable (which we might expect since both models are consistent for the density). Despite these similarities, there are at least two major differences between MFMs and the usual nonparametric mixture models (e.g., DPMs and Pitman–Yor process mixtures (PYMs)): (1) The prior distribution of the number of clusters is very different. In an MFM, one has complete control over the distribution of the number of components — leading to control over the distribution of the number of clusters — and (under the prior) as the sample size grows, the number of clusters converges to a finite value with probability 1 (but has no a priori bound). In a DPM or PYM (with discount parameter σ ∈ [0, 1)), the distribution of the number of 72 clusters follows a particular parametric form, and (under the prior) the number of clusters grows to infinity with probability 1. (2) Given the number of clusters, the prior distribution of the cluster sizes is very different. In particular (under the prior), in an MFM the sizes of the clusters are all the same order of magnitude (asymptotically with respect to the number of samples), while in a DPM or PYM (with σ ∈ [0, 1)) many of the clusters are negligible in size relative to the largest clusters. These differences in the priors carry over to differences in certain aspects of the posteriors. Partitions sampled from the DPM posterior tend to have multiple small transient clusters, which are not seen in samples from the MFM posterior; conse- quently, the DPM posterior on the number of clusters tends to put more probability mass on higher numbers of clusters than the MFM. Further, given the number of clusters, the entropy of a partition sampled from the DPM posterior tends be lower than one sampled from the MFM posterior. One could argue that, compared to the DPM, the MFM is a more natural Bayesian approach for a data distribution of unknown complexity (if something is unknown, put a prior on it), and consequently, that it enjoys greater interpretability and conceptual simplicity. Further, in an MFM, the parameter space is a countable union of finite-dimensional spaces, while in a DPM or PYM, the parameter space is infinite-dimensional. There are also some disadvantages to using MFMs. A minor inconvenience is that the coefficients of the partition distribution need to be precomputed when doing in- ference. A more significant issue is that the mixing time of incremental Markov chain Monte Carlo (MCMC) samplers can be worse than for the DPM, since the MFM 73 dislikes small clusters, making it difficult to create or destroy substantial clusters by moving one element at a time. A natural solution to this mixing issue would be to use a split-merge sampler. In general, we would not say that one model is uniformly better than the other — the choice of whether to use an MFM or a DPM or PYM (or something else) should be made based on the appropriateness of the model to the application at hand. 4.1.1 Mixture of finite mixtures (MFM) Consider the following variable-dimension mixture model: K ∼ p(k), where p(k) is a p.m.f. on {1, 2, . . . } (π1 , . . . , πk ) ∼ Dirichletk (γ, . . . , γ), given K = k iid Z1 , . . . , Zn ∼ π, given π (4.1.1) iid θ1 , . . . , θk ∼ H, given K = k Xj ∼ fθZj independently for j = 1, . . . , n, given θ1:K , Z1:n . Here, H is a prior or “base measure” on Θ ⊂ R` , and {fθ : θ ∈ Θ} is a family of probability densities with respect to a sigma-finite measure λ on X ⊂ Rd . (As usual, we give Θ and X the Borel sigma-algebra.) In a typical application, the values X1 , . . . , Xn would be observed, and all other variables would be hidden/latent. See Figure 4.1 for a graphical model representation. We refer to this as a mixture of finite mixtures (MFM) model, since there does not seem to be a well-established, unambiguous name. This type of model has been studied by many authors (Nobile, 1994, Phillips and Smith, 1996, Richardson 74 K π θ Zj Xj n Figure 4.1: Graphical model for the MFM. and Green, 1997, Green and Richardson, 2001, Stephens, 2000, Nobile and Fearnside, 2007) and has found many applications, including neuron spike classification (Nguyen et al., 2003), survival analysis (Marin et al., 2005), DNA restriction mapping (Lee et al., 1998), and gene expression profiling (Vogl et al., 2005). The most common approach to inference in such a model is reversible jump Markov chain Monte Carlo (RJMCMC) (Richardson and Green, 1997, Green and Richardson, 2001), although other strategies have been proposed as well (Nobile, 1994, Phillips and Smith, 1996, Stephens, 2000, Nobile and Fearnside, 2007). Re- versible jump is a very general and widely-applicable method, but it is not a “black box”. In contrast, a nice aspect of many DPM samplers is that they are fairly generic. In their paper developing a reversible jump sampler for the DPM, Green and Richardson (2001) pointed out that it would be interesting if, conversely, DPM- style samplers could be used for the MFM model: In view of the intimate correspondence between DP and DMA models dis- cussed above, it is interesting to examine the possibilities of using either class of MCMC methods for the other model class. We have been unsuc- cessful in our search for incremental Gibbs samplers for the DMA models, but it turns out to be reasonably straightforward to implement reversible 75 jump split/merge methods for DP models. (Green and Richardson, 2001) Note: They seem to use the term Dirichlet-multinomial allocation (DMA) to refer both to the MFM model above with k random, and to the corresponding model with k fixed (the latter of which is the more standard usage, (Chen, 1980)); to avoid confusion, we say MFM instead of DMA for the k random case. The key to many DPM samplers is that the model can be characterized by a nice distribution on combinatorial structures — namely, the Chinese restaurant pro- cess. It turns out, as we observe below, that MFMs have a similar characterization, enabling the application of the same sampling algorithms to them. It is important to note that in the MFM, we assume a symmetric Dirichlet with a single parameter γ > 0 not depending on k. This assumption is key to deriving a simple form for the partition probability distribution and the resulting restaurant process. Assuming symmetry in the distribution of π is quite natural, since the distri- bution of X1 , . . . , Xn under any asymmetric distribution on π would be the same as if this were replaced by its symmetrized version (e.g., if the entries of π were uniformly permuted). Assuming the same γ for all k is a genuine restriction, albeit a fairly natural one, often made in such models even when not strictly necessary (Nobile, 1994, Richardson and Green, 1997, Green and Richardson, 2001, Stephens, 2000, Nobile and Fearnside, 2007). Note that prior information about the relative sizes of the mixing weights π1 , . . . , πk can be introduced through γ — roughly speaking, small γ favors lower entropy π’s, while large γ favors higher entropy π’s. On the other hand, we will make very few assumptions on p(k), the distribution of the number of components. For practical purposes, we need the infinite series P∞ k=1 p(k) to converge to 1 reasonably quickly, although any choice of p(k) arising in 76 practice should not be a problem. For certain theoretical purposes — in particular, consistency for the number of components — it is desirable to have p(k) > 0 for all {1, 2, . . . }. This chapter is organized as follows. In Section 4.2, we derive a number of prop- erties of the MFM model, including the partition distribution, restaurant process, random discrete measure formulation and other equivalent models, stick-breaking representation in a special case, density estimates, and various basic properties. In Section 4.3, we discuss existing results on the posterior asymptotics of the density, the mixing distribution, and the number of components, and we derive some asymp- totic properties of the MFM model. In Section 4.4, we show how the properties derived in Section 4.2 facilitate inference with Markov chain Monte Carlo (MCMC) algorithms; in particular, we show how partition-based samplers for the DPM — in the case of both conjugate and non-conjugate priors — can be directly adapted to the MFM. 4.2 Properties of MFMs 4.2.1 Partition distribution The primary observation on which our development relies is that the distribution on partitions induced by an MFM takes a form which is simple enough that it can be easily computed. Consider the MFM model defined above. Let C denote the unordered partition of [n] := {1, . . . , n} induced by Z1 , . . . , Zn ; in other words, C = {Ei : |Ei | > 0} where Ei = {j : Zj = i} for i = 1, 2, . . . . Then, as we show 77 below, Y p(C) = Vn (t) γ (|c|) (4.2.1) c∈C where t = |C| is the number of parts in the partition, ∞ X k(t) Vn (t) = p(k), (4.2.2) k=1 (γk)(n) x(m) = x(x + 1) · · · (x + m − 1), and x(m) = x(x − 1) · · · (x − m + 1). (By convention, x(0) = 1 and x(0) = 1.) Viewed as a function of the part sizes (|c| : c ∈ C), Equation 4.2.1 is an “ex- changeable partition probability function” (EPPF) in the terminology of Pitman (1995, 2006) — that is, it is a symmetric function of the part sizes (and depends on C only through t = |C| and the part sizes). Consequently, C is an exchangeable random partition of [n]; that is, its distribution is invariant under permutations of [n] (alter- natively, this can be seen directly from the definition of the model, since Z1 , . . . , Zn are exchangeable). More specifically, Equation 4.2.1 is a member of the family of Gibbs partition distributions (Pitman, 2006). This can be viewed as a special case of the general results of Gnedin and Pitman (2006) characterizing the extremal points of the space of Gibbs partition distributions; also see Lijoi et al. (2008), Ho et al. (2007), Cerquetti (2008, 2011), Gnedin (2010), and Lijoi and Pr¨ unster (2010) for further results on Gibbs partitions. However, the utility of this representation for inference in a wide variety of variable-dimension mixture models does not seem to have been previously explored in the literature. The derivation of Equation 4.2.1 is straightforward, and we provide it here. Let- ting Ei = {j : zj = i} as above, and writing C(z) for the partition induced by 78 z = (z1 , . . . , zn ), by Dirichlet-multinomial conjugacy we have Qk Γ(|Ei | + γ) Z Γ(kγ) 1 Y p(z|k) = p(z|π)p(π|k)dπ = i=1 = γ (|c|) , Γ(γ)k Γ(n + kγ) (kγ)(n) c∈C(z) provided that p(k) > 0. It follows that for any partition C of [n], X p(C|k) = p(z|k) z∈[k]n : C(z)=C n o 1 Y (|c|) = # z ∈ [k]n : C(z) = C γ (γk)(n) c∈C k(t) Y (|c|) = γ , (4.2.3) (γk)(n) c∈C k   where t = |C|, since # z ∈ [k]n : C(z) = C = t t! = k(t) . Finally, ∞ ∞ X Y X k(t) Y p(C) = p(C|k)p(k) = γ (|c|) p(k) = Vn (t) γ (|c|) , k=1 c∈C k=1 (γk)(n) c∈C with Vn (t) as in Equation 4.2.2 above. 4.2.2 Equivalent models For various purposes, it useful to write the model (Equation 4.1.1) in various equiv- iid alent ways. First, note that instead of sampling only θ1 , . . . , θk ∼ H given K = k, iid we could simply sample θ1 , θ2 , . . . ∼ H independently of K, and the distribution of X1:n would be the same; the graphical model for this version is in Figure 4.2(a). Now, Z1:n determines which subset of the i.i.d. variables θ1 , θ2 , . . . will actually be used, and the indices of this subset are independent of θ1:∞ ; hence, denoting these random indices I1 < · · · < IT , we have that θI1 , . . . , θIT |Z1:n are i.i.d. from H. 79 K π θ1:∞ Zj Xj n (a) K C φ π C φ Xj n Zj Xj n (c) (b) Figure 4.2: Alternative graphical models for the MFM. From this we see that we could equivalently write the model in terms of C as follows, resulting in the same distribution for C and X1:n (see Figure 4.2(b,c)): C ∼ p(C), with p(C) as in Equation 4.2.1 iid φc ∼ H for c ∈ C, given C (4.2.4) Xj ∼ fφc independently for j ∈ c, c ∈ C, given φ, C. Note that for notational convenience, here, φ = (φc : c ∈ C) is a tuple of t = |C| parameters φc ∈ Θ, one for each part c ∈ C. This representation of the model is particularly useful for doing inference, since one does not have to deal with cluster labels or empty components. The formulation of models starting from a partition distribution has been a fruitful approach, exem- 80 plified by the development of product partition models (Quintana and Iglesias, 2003, Quintana, 2006, Dahl, 2009, Park and Dunson, 2010, M¨ uller et al., 2011, Barry and Hartigan, 1992); also see M¨ uller and Quintana (2010) for a review. 4.2.3 Various basic properties We list here some basic properties of the MFM model above. The derivations are all straightforward, but for completeness we provide them in Section 4.2.10.1. Denoting R Q  xc = (xj : j ∈ c) and m(xc ) = Θ j∈c fθ (xj ) H(dθ) (with the convention that m(x∅ ) = 1), we have Y p(x1:n |C) = m(xc ). (4.2.5) c∈C The number of components K and the number of clusters T = |C| are related by k(t) X Y (|c|) p(t|k) = γ , (4.2.6) (γk)(n) c∈C C:|C|=t 1 k(t) p(k|t) = p(k), (4.2.7) Vn (t) (γk)(n) where in p(t|k), the sum is over partitions C of [n] such that |C| = t. The formula for p(k|t) is required for doing inference about the number of components K based on posterior samples of C; fortunately, it is easy to compute. We have the conditional independence relations C ⊥ K | T, (4.2.8) X1:n ⊥ K | T. (4.2.9) 81 4.2.4 The coefficients Vn (t) γ (|c|) , where t = |C| and Q Recall that p(C) = Vn (t) c∈C ∞ X k(t) Vn (t) = p(k). k=1 (γk)(n) The numbers Vn (t) satisfy the recursion Vn+1 (t + 1) = Vn (t)/γ − (n/γ + t)Vn+1 (t) (4.2.10) for any 0 ≤ t ≤ n and γ > 0; this is easily seen by plugging the identity k(t+1) = (γk + n)k(t) /γ − (n/γ + t)k(t) into the expression for Vn+1 (t + 1). This recursion is a special case of a more general recursion for a large class of partition distributions (Gnedin and Pitman, 2006). In the case of γ = 1, Gnedin (2010) has given a nice example of a distribution on K for which both p(k) and Vn (t) have closed-form expressions. In previous work on the MFM model, it has been common for p(k) to be chosen to be proportional to a Poisson distribution restricted to a subset of the positive integers (Phillips and Smith, 1996, Stephens, 2000, Nobile and Fearnside, 2007), and Nobile (2005) has proposed a theoretical justification for this choice. Interestingly, the model has some nice mathematical properties if one instead chooses K − 1 to be given a Poisson distribution, that is, p(k) = Poisson(k − 1|λ) for some λ > 0. One example of this arises here (for another example, see Section 4.2.9): it turns out that 82 if p(k) = Poisson(k − 1|λ) and γ = 1 then (see Section 4.2.10.2 for details) n 1 X  Vn (0) = n 1 − p(k) . (4.2.11) λ k=1 However, to do inference, it is not necessary to choose p(k) to have any particular form. To do inference, we just need to be able to compute p(C), and in turn, we need to be able to compute Vn (t). To this end, note that k(t) /(γk)(n) ≤ k t /(γk)n , and thus the infinite series for Vn (t) converges rapidly when t  n. It always converges to a finite value when 1 ≤ t ≤ n; this is clear from the fact that p(C) ∈ [0, 1]. This finiteness can also be seen directly from the series since k t /(γk)n ≤ 1/γ n , and in fact, this shows that the series for Vn (t) converges at least as rapidly (up to a constant) as the series ∞ P k=1 p(k) converges to 1. Hence, for any reasonable choice of p(k) (i.e., not having an extraordinarily heavy tail), Vn (t) can easily be numerically approximated to a high level of precision. For instance, in Chapter 5, we apply the model using   1 if k ∈ {1, ..., 30}  p(k) ∝  1/(k − 30)2 if k > 30,  and computing the coefficients Vn (t) presents no problems, even though this has a fairly heavy tail (the mean is infinite) and does not take an analytic form. Note that Vn (t) is often exceedingly small, making it necessary to use a standard technique such as the “log-sum-exp trick” in order to avoid numerical underflow (that is, represent each term in log space and use log(exp(a)+exp(b)) = log(exp(a−m)+exp(b−m))+m, where m = max{a, b}, to compute the log of the sum). Although we find it easier to use this numerical approach for all Vn (t) values, it 83 would also be possible to use the recursion in Equation 4.2.10 (taking care, again, regarding underflow) and numerically approximate only the required values of Vn (0) (or use the expression for Vn (0) in Equation 4.2.11 if K − 1 ∼ Poisson(λ) and γ = 1). In the course of doing inference, n is usually fixed and the MCMC sampling methods we describe below typically only require one to compute Vn (t) for t = 1, . . . , tmax for some relatively small tmax , since higher values of t have such low probability that in practice they are never visited by the sampler. Consequently, it suffices to precompute Vn (1), . . . , Vn (tmax ) for some suitably chosen tmax , and this takes a negligible amount of time compared to running the sampler. Also, on a practical note, the samplers described below involve moving one ele- ment at a time between clusters, so in order to achieve irreducibility of the Markov chain, it is necessary to have {t ∈ {1, 2, . . . } : Vn (t) > 0} be a block of consecutive integers. In fact, it turns out that this is always the case (and this block includes t = 1), since for any k such that p(k) > 0, we have Vn (t) > 0 for all t = 1, . . . , k. When we come to the restaurant process in Section 4.2.6 below, it will become clear that precomputing the Vn (t)’s effectively amounts to precomputing the prob- abilities for “choosing a new table”. The idea of precomputing these probabilities has previously been applied to DPMs with a prior on the concentration parameter α (MacEachern, 1998), in which the pre-computation consists of numerically integrat- ing α out. In the MFM model as described above (Equation 4.1.1), the Dirichlet parameter γ is fixed, but in some cases it might be interesting to put a prior on γ, to allow greater flexibility regarding the relative sizes of the components. (Thanks to Vinayak Rao for this suggestion.) Although we have not tried this, it should be straightforward to 84 periodically update γ with a Metropolis-Hastings step and recompute the required Vn (t) coefficients, possibly using a simple caching scheme. 4.2.5 Self-consistent marginals For each n = 1, 2, . . . , let qn (C) denote the distribution on partitions of [n] as defined above (Equation 4.2.1). This family of partition distributions has the property that qm coincides with the “marginal” distribution on partitions of [m] induced by qn when n ≥ m; in other words, drawing a sample from qn and removing elements m + 1, . . . , n from it yields a sample from qm . This can be seen directly from the model definition (Equation 4.1.1), since C is the partition induced by the Z’s, and the distribution of Z1:m is the same when the model is defined with any n ≥ m. This property is sometimes referred to as consistency in distribution (Pitman, 2006). Using Kolmogorov’s extension theorem (e.g., Durrett, 1996), one can show that this implies the existence of a unique probability distribution on partitions of the positive integers Z>0 = {1, 2, . . . } such that the marginal distribution on partitions of [n] is qn for all n = 1, 2, . . . . A random partition of Z>0 from such a distribution is a combinatorial stochastic process; for background, see Pitman (2006). It is instructive to also verify this self-consistency property for MFMs by the more tedious approach of directly marginalizing. Suppose C is a partition of [n], denote by C∗ the partition of [n + 1] obtained from C by making n + 1 a singleton, and for c ∈ C denote by Cc the partition of [n + 1] obtained from C by adding n + 1 to P c. We would like to show that qn+1 (C∗ ) + c∈C qn+1 (Cc ) = qn (C), and indeed, letting 85 t = |C|, 0 0 Y X Y Vn+1 (t + 1)γ γ (|c |) + Vn+1 (t)(γ + |c|) γ (|c |) c0 ∈C c∈C c0 ∈C   Y  Y = Vn+1 (t + 1)γ + Vn+1 (t)(tγ + n) γ (|c|) = Vn (t) γ (|c|) , c∈C c∈C the last step following from the recursion, Equation 4.2.10. The property follows by induction. 4.2.6 Restaurant process / P´ olya urn scheme Pitman (1996) considered a general class of “restaurant processes”, or P´olya urn schemes, corresponding to exchangeable partition probability functions (EPPFs). Following this general rule, we can derive a restaurant process for the MFM. For any vector of random variables (Y1 , . . . , Yn ), we can obtain a sample by first drawing Y1 , then Y2 |Y1 , and so on, up to Yn |Yn−1 , . . . , Y1 . In the same way, if we write C1 , . . . , Cn for the partitions of [1], . . . , [n], respectively, induced by a partition Cn of [n], then for any distribution on partitions of [n], we can obtain a sample Cn by sampling C1 , C2 |C1 , and so on. So, in a crude sense, any distribution on partitions has a “restaurant process” describing this sequence of conditional distributions, but it may depend on n and it may be quite complicated. However, since we are in the fortunate circumstance that our partition distri- bution has the self-consistency property described in the previous section, when Cn ∼ qn , we have Cm ∼ qm for m = 1, . . . , n (where qn denotes the distribution on partitions as in Equation 4.2.1). Consequently, there is a single restaurant process 86 that works for all n, and it takes a simple form involving qm at step m, since p(Cm |Cm−1 , . . . , C1 ) = p(Cm |Cm−1 ) ∝ qm (Cm ) I(Cm \ m = Cm−1 ), where C \ m denotes C with element m removed. γ (|c|) , we have, letting t = |Cm−1 |, Q Recalling that qm (Cm ) = Vm (|Cm |) c∈Cm   Vm (t + 1)γ if m is a singleton in Cm , i.e., {m} ∈ Cm  p(Cm |Cm−1 ) ∝  Vm (t)(γ + |c|) if c ∈ Cm−1 and c ∪ {m} ∈ Cm ,  for Cm such that Cm \ m = Cm−1 (and p(Cm |Cm−1 ) = 0 otherwise). In other words, we have the following restaurant process for the MFM: • The first customer sits at a table: C1 = {{1}}. • For n = 2, 3, . . . , the nth customer sits . . . at an existing table c ∈ Cn−1 with probability ∝ |c| + γ Vn (t + 1) at a new table with probability ∝ γ Vn (t) where t = |Cn−1 |. Clearly, this bears a resemblance to the Chinese restaurant process, in which the nth customer sits at an existing table c with probability ∝ |c| or at a new table with probability ∝ α (the concentration parameter) (Blackwell and MacQueen, 1973). 87 G C φ Xj βj n n (a) (b) Figure 4.3: Graphical models for (a) the random discrete measure formulation, and (b) the equivalent construction of the distribution of β1 , . . . , βn . 4.2.7 Random discrete measure formulation The MFM can also be formulated starting from a distribution on discrete measures that is analogous to the Dirichlet process. Let K ∼ pK (k) (π1 , . . . , πk ) ∼ Dirichletk (γ, . . . , γ), given K = k iid θ1 , . . . , θk ∼ H, given K = k K X G= πi δθi i=1 and denote the distribution of G by M(pK , γ, H). Then the distribution of X1:n is the same as before if we take X1 , . . . , Xn |G i.i.d. from the resulting mixture, namely, Z K X fG (x) := fθ (x)G(dθ) = πi fθi (x). i=1 So, in this notation, the model is (Figure 4.3(a)): G ∼ M(pK , γ, H) X1 , . . . , Xn ∼ fG , given G. 88 When H is continuous, the random discrete measure G belongs to the class of “species sampling” models studied by Pitman (1996), who proved many properties regarding this type of distribution. For completeness, we derive the prediction rule for G, along with some other properties. We refer to Pitman (1996), Hansen and Pitman (2000), Ishwaran and James (2003, 2001), Lijoi et al. (2005a, 2007), and Lijoi et al. (2008) for more background on species sampling models and further examples. iid Let G ∼ M(pK , γ, H) and let β1 , . . . , βn ∼ G, given G. Then, by construc- tion, the joint distribution of (β1 , . . . , βn ) (with G marginalized out) is the same as (θZ1 , . . . , θZn ) in the original model (Equation 4.1.1). Further, letting C denote the partition induced by Z1 , . . . , Zn as usual, we have (θZ1 , . . . , θZn ) = (φc1 , . . . , φcn ) where cj is defined to be the c ∈ C such that j ∈ c, and φc is defined to be equal to the θi such that Zj = i for all j ∈ c. By the same argument as in Section 4.2.2, (φc : c ∈ C) are i.i.d. from H given C. Therefore, we have the following equivalent construction for (β1 , . . . , βn ): Cn ∼ qn , with qn as in Section 4.2.5 iid φc ∼ H for c ∈ Cn , given Cn βj = φc for j ∈ c, c ∈ Cn , given Cn , φ. See Figure 4.3(b). Due to the self-consistency property of q1 , q2 , . . . , we can sample Cn , (φc : c ∈ Cn ), β1:n sequentially for n = 1, 2, . . . by sampling from the restaurant process for Cn |Cn−1 , sampling φ{n} from H if n is seated at a new table (or setting φc ∪{n} = φc if n is seated at c ∈ Cn−1 ), and setting βn accordingly. In particular, if the base measure H is continuous (that is, if H({θ}) = 0 for all θ ∈ Θ), then conditioning on β1:n−1 is the same as conditioning on Cn−1 , (φc : 89 c ∈ Cn−1 ), β1:n−1 , so we can sample βn |β1:n−1 in the same way as was just described. Therefore, when H is continuous, the distribution of βn given β1 , . . . , βn−1 is propor- tional to Vn (t + 1) X γ H+ (|c| + γ)δφc , (4.2.12) Vn (t) c∈C n−1 where t = |Cn−1 |, with the Cn−1 and (φc : c ∈ Cn−1 ) induced by β1 , . . . , βn−1 , or equivalently, n−1 t Vn (t + 1) X X γ H+ δβ j + γ δβj∗ , Vn (t) j=1 j=1 where β1∗ , . . . , βt∗ are the distinct values taken by β1 , . . . , βn−1 . For comparison, if G ∼ DP(α, H) instead, the distribution of βn given β1 , . . . , βn−1 is proportional to n−1 X αH + δβ j . j=1 This is the classic P´olya urn scheme for the Dirichlet process, described by Blackwell and MacQueen (1973). 4.2.8 Density estimates In this section, we derive formulas for density estimation with the MFM using sam- ples of C, φ | x1:n , or, if the single-cluster marginals can be computed, using samples of C | x1:n . Using the restaurant process (at the end of Section 4.2.6 above), it is straightforward to show that, if C is a partition of [n] and φ = (φc : c ∈ C) then Vn+1 (t + 1) X p(xn+1 | C, φ, x1:n ) ∝ γ m(xn+1 ) + (|c| + γ)fφc (xn+1 ) (4.2.13) Vn+1 (t) c∈C 90 where t = |C|, and, using the recursion for Vn (t) (Equation 4.2.10), this is normalized when multiplied by Vn+1 (t)/Vn (t). Further, Vn+1 (t + 1) X m(xc ∪{n+1} ) p(xn+1 | C, x1:n ) ∝ γ m(xn+1 ) + (|c| + γ) , (4.2.14) Vn+1 (t) c∈C m(x c) with the same normalization constant. Therefore, when the single-cluster marginals m(xc ) can be easily computed, Equation 4.2.14 can be used to estimate the posterior predictive density p(xn+1 |x1:n ) based on samples from C | x1:n . When m(xc ) cannot be easily computed, Equation 4.2.13 can be used to estimate p(xn+1 |x1:n ) based iid on samples from C, φ | x1:n , along with samples θ1 , . . . , θN ∼ H to approximate m(xn+1 ) ≈ N1 N P i=1 fθi (xn+1 ). The posterior predictive density is, perhaps, the most natural estimate of the density. However, following Green and Richardson (2001), another way to obtain a fairly natural estimate is by assuming that element n+1 is added to an existing clus- ter; this will be very similar to the posterior predictive density when n is sufficiently large. To this end, we define p∗ (xn+1 | C, φ, x1:n ) = p(xn+1 | C, φ, x1:n , |Cn+1 | = |C|), where Cn+1 is the partition of [n + 1], and observe that X |c| + γ p∗ (xn+1 | C, φ, x1:n ) = fφc (xn+1 ) c∈C n + γt where t = |C| (Green and Richardson, 2001). Using this, we can estimate the density by N 1 X p∗ (xn+1 | C (i) , φ(i) , x1:n ), (4.2.15) N i=1 where (C (1) , φ(1) ), . . . , (C (N ) , φ(N ) ) are samples from C, φ | x1:n . 91 The corresponding expressions for the DPM are all very similar, with the obvious changes, using its restaurant process instead. These formulas are conditional on additional parameters such as γ for the MFM, and α for the DPM. If priors are placed on such parameters and they are sampled along with C and φ given x1:n , then the posterior predictive density can be estimated using the same formulas as above, but also using the posterior samples of these additional parameters. 4.2.9 Stick-breaking representation The Dirichlet process has an elegant stick-breaking representation for the mixture weights π1 , π2 , . . . (Sethuraman, 1994, Sethuraman and Tiwari, 1981). This extraor- dinarily clarifying perspective has inspired a number of other nonparametric models (MacEachern et al., 1999, MacEachern, 2000, Hjort, 2000, Ishwaran and Zarepour, 2000, Ishwaran and James, 2001, Griffin and Steel, 2006, Dunson and Park, 2008, Chung and Dunson, 2009, Rodriguez and Dunson, 2011, Broderick et al., 2012), has provided insight into the properties of other models (Favaro et al., 2012, Teh et al., 2007, Thibaux and Jordan, 2007, Paisley et al., 2010), and has been used to develop efficient inference algorithms (Ishwaran and James, 2001, Blei and Jordan, 2006, Papaspiliopoulos and Roberts, 2008, Walker et al., 2007, Kalli et al., 2011). In a certain special case — namely, when p(k) = Poisson(k − 1|λ) and γ = 1 — we have noticed that the MFM also has a nice stick-breaking representation. (There may also be other cases of which we are not yet aware.) This is another example of the nice mathematical properties resulting from this choice of p(k) and γ. Consider the following “stick-breaking” procedure: 92 Take a unit-length stick, and break off i.i.d. Exponential(λ) pieces until you run out of stick. iid e = min{j : Pj εi > 1}, In other words, let ε1 , ε2 , . . . ∼ Exponential(λ), define K i=1 eKe = 1 − K−1 e − 1 and π P and set π ei = εi for i = 1, . . . , K i=1 π ei . Then the resulting e stick lengths π e have the same distribution as the mixture weights π in the MFM model when p(k) = Poisson(k − 1|λ) and γ = 1. This is a consequence of a standard construction of a Poisson process; for details, see Section 4.2.10.3. This suggests an alternative way of constructing a variable-dimension mixture: take any sequence of nonnegative random variables (ε1 , ε2 , . . . ) (not necessarily in- dependent or identically distributed) such that ∞ P i=1 εi > 1 with probability 1, and define K e and π e as above. Although the distribution of K e and π e may be complicated, in some cases it might still be possible to do inference based on the stick-breaking representation. This might be an interesting way to introduce different kinds of prior information on the mixture weights, however, we have not explored this possibility. 93 4.2.10 Proofs and details 4.2.10.1 Derivation of various basic properties Here, we derive the properties listed in Section 4.2.3. Abbreviate x = x1:n , z = z1:n , and θ = θ1:k , and assume p(z, k) > 0. Then p(x|θ, z, k) = ki=1 j∈Ei fθi (xj ) and Q Q Z k Z h Y Y i p(x|z, k) = p(x|θ, z, k)p(dθ|k) = fθi (xj ) H(dθi ) Θk i=1 Θ j∈Ei k Y Y = m(xEi ) = m(xc ). i=1 c∈C(z) Since this last expression depends only on z, k through C = C(z), we have Y m(xc ) = p(x|z, k) = p(x|z) = p(x|C), c∈C k(t) γ (|c|) (where t = Q establishing Equation 4.2.5. Next, recall that p(C|k) = (γk)(n) c∈C |C|) from Equation 4.2.3, and thus X k(t) X Y (|c|) p(t|k) = p(C|k) = γ , (γk)(n) c∈C C:|C|=t C:|C|=t (where the sum is over partitions C of [n] such that |C| = t) establishing Equation 4.2.6. Equation 4.2.7 follows, since k(t) p(k|t) ∝ p(t|k)p(k) ∝ p(k), (γk)(n) 94 (provided p(t) > 0) and the normalizing constant is precisely Vn (t). To see that C ⊥ K | T (Equation 4.2.8), note that if t = |C| then p(C, t|k) p(C|k) p(C|t, k) = = , p(t|k) p(t|k) (provided p(t, k) > 0) and due to the form of p(C|k) and p(t|k) just above, this quantity does not depend on k; hence, p(C|t, k) = p(C|t). To see that X ⊥ K | T (Equation 4.2.9), note that from the graphical model in Figure 4.2(b), X ⊥ K | C; using this in addition to C ⊥ K | T , we have X X p(x|t, k) = p(x|C, t, k)p(C|t, k) = p(x|C, t)p(C|t) = p(x|t). C:|C|=t C:|C|=t 4.2.10.2 Vn (0) in the Poisson case with γ = 1 To see Equation 4.2.11, observe that when p(k) = Poisson(k − 1|λ) and γ = 1, we have ∞ ∞ ∞ X 1 e−λ λk−1 −n X e−λ λk+n−1 −n X Vn (0) = (n) (k − 1)! = λ = λ p(k), k=1 k k=1 (k + n − 1)! k=n+1 P∞ Pn and k=n+1 p(k) = 1 − k=1 p(k). 4.2.10.3 Proof of the stick-breaking representation iid e = min{j : Pj εi > 1}, π Let ε1 , ε2 , . . . ∼ Exponential(λ), K i=1 ei = εi for eKe = 1 − K−1 e − 1 and π Pe i = 1, . . . , K i=1 π ei . Pj Let Sj = i=1 εi and define N (τ ) = #{j : Sj ≤ τ }. Then N (τ ) is a Poisson pro- 95 cess with rate λ (e.g., Durrett, 1996). Hence, N (τ ) ∼ Poisson(λτ ), and in particular, since K e − 1 ∼ Poisson(λ). e = N (1) + 1 a.s., we have K It is a standard result that (S1 , . . . , Sm )|N (τ ) = m has the same joint dis- tribution as the order statistics (U(1) , . . . , U(m) ) of i.i.d. uniform random variables U1 , . . . , Um ∼ Uniform[0, τ ]. Considering the case of τ = 1, the conditional density of (S1 , . . . , Sm )|N (1) = m is therefore m! I(0 < s1 < · · · < sm < 1), and since the change of variables from (S1 , . . . , Sm ) to (ε1 , . . . , εm ) has Jacobian equal to 1, it fol- lows that (ε1 , . . . , εm )|N (1) = m has density m! I( m P i=1 ε < 1, εi > 0 ∀i). Therefore, e|K taking k = m + 1, we see that π e = k has the uniform k-dimensional Dirichlet distribution. 4.3 Asymptotics In this section, we first summarize and compare asymptotic results for posterior infer- ence in the mixture of finite mixtures (MFM) and Dirichlet process mixture (DPM) models with respect to three objects of possible interest: the density, the mixing distribution, and the number of components. We then provide some basic asymp- totic results for various aspects of the MFM: the asymptotic relationship between the number of components and the number of clusters, the asymptotic behavior of Vn (t), and the form of the conditional distribution on part sizes given the number of parts. 96 4.3.1 Posterior asymptotics To set things up, we first pose the questions of interest from a general perspective. Suppose {fθ : θ ∈ Θ} is a parametric family of densities with respect to a measure on X ⊂ Rd . Given a discrete probability measure G = ∞ P i=1 πi δθi (where θ1 , θ2 , . . . ∈ Θ, P δθ is the unit point mass at θ ∈ Θ, and π1 , π2 , . . . ≥ 0 such that i πi = 1), let fG denote the density of the resulting mixture, that is, Z ∞ X fG (x) = fθ (x) G(dθ) = πi fθi (x), Θ i=1 and let s(G) denote the number of elements in the support of G, that is, s(G) = | support(G)| ∈ {1, 2, . . . } ∪ {∞}. Suppose we observe data X1 , . . . , Xn ∼ f∗ for some true density f∗ . The Bayesian approach is to put a prior on f∗ , and given X1 , . . . , Xn , make posterior inferences about f∗ . Mixtures of the form fG above tend to generate a very flexible class of densities, and considerable success has been obtained by constructing Bayesian models following this form, Dirichlet process mixtures being the prime example. In the density estimation setting, the true density f∗ does not have to take the form fG in order to obtain posterior guarantees — only general regularity conditions are required; meanwhile, when estimating the mixing distribution and the number of clusters, it is typically necessary to assume f∗ is of the form fG . Further, although it is irrelevant for density estimation, to have any hope of consistency for the mixing distribution, it is necessary to assume that the family {fθ } is mixture identifiable in the following sense: if s(G), s(G0 ) < ∞ and fG = fG0 a.e., then G = G0 . Note that we are not constraining the θ’s to be distinct and 97 ordered (nor constraining the π’s to be positive), so the π’s and θ’s will be non- identifiable; however, the assumption here is identifiability of the measure G, and this is satisfied for many commonly-used families {fθ }, such as multivariate Gaussian, Gamma, Poisson, Exponential, Cauchy (Teicher, 1963, Yakowitz and Spragins, 1968) in the continuous case, and Poisson, Geometric, Negative Binomial, or any power- series distribution (Sapatinas, 1995) in the discrete case. It should be said, perhaps, at this point, that inferences about the mixing distri- bution or the number of components will typically be sensitive to misspecification of the family {fθ }; for instance, any sufficiently regular density can be approximated arbitrarily well by a mixture of Gaussians, so even if the true density is very close to but not exactly a finite mixture Gaussians, the model will introduce “extra” com- ponents in order to fit the data distribution. On the other hand, when making inferences about the density, there is no assumption that the true density f∗ has any particular form (other than general regularity conditions), so misspecification is not a serious issue. Nonetheless, it is desirable to obtain consistency results for the mixing distribu- tion and the number of components when the model is well-specified, since lack of consistency under ideal conditions would be a red flag. Given a prior on G, then, consider the following questions. (1) Density estimation. Does the posterior on the density concentrate at the true density, and if so, at what rate of convergence? That is, for an appropriate metric dist(·, ·), does P Pmodel (dist(fG , f∗ ) < ε | X1:n ) −−data −→ 1 n→∞ 98 for all ε > 0, and if so, how rapidly can we let εn go to 0 and still have P Pmodel (dist(fG , f∗ ) < Cεn | X1:n ) −−data −→ 1 n→∞ for some C > 0? The standard choices of metric here are L1 and Hellinger. (2) Mixing distribution. Does the posterior on the mixing distribution concen- trate at the true mixing distribution (assuming there is one)? That is, for an appropriate metric dist(·, ·), does P Pmodel (dist(G, G∗ ) < ε | X1:n ) −−data −→ 1 n→∞ for all ε > 0? The Wasserstein metric has been used in this context (Nguyen, 2013). (3) Number of components. Does the posterior on the number of components concentrate at the true number of components, for data from a finite mixture? That is, does P Pmodel (s(G) = s(G∗ ) | X1:n ) −−data −→ 1? n→∞ Also of interest, due to its (albeit, questionable) use in practice for inference about the number of components, is the number of clusters T , defined — when the prior on G is introduced via priors on π and θ — as the number of distinct values of the latent allocation variables Z1 , . . . , Zn ∼ π, where Xi ∼ fθZi . Does the posterior on the number of clusters T concentrate at the true number of components? That is, does P Pmodel (T = s(G∗ ) | X1:n ) −−data −→ 1? n→∞ In varying degrees of generality, answers to these questions for both MFMs and 99 DPMs have been provided in the literature, as summarized in Table 4.1. We provide references for these results in the sections below. It should be emphasized that most of these results — especially regarding rates of convergence — have only been shown for certain types of models; the table is simply meant to give a rough idea of the general picture. Table 4.1: Summary of known posterior consistency results (Yes = consistent, No = not consistent); see text for explanation. MFMs DPMs Density estimation Yes (optimal rate) Yes (optimal rate) Mixing distribution Yes Yes Number of components Yes No Here, by “optimal rate”, we mean optimal up to a logarithmic factor in n. Also, note that if the true number of components is finite, then the DPM posterior on the number of components is trivially inconsistent, since the number of components is infinite with probability 1; the inconsistency in the table refers to inconsistency of the number of clusters T . 4.3.1.1 Density estimation Density estimation is perhaps the primary type of application for a flexible Bayesian mixture model, such as a DPM or MFM. In this setting, the true density f∗ is assumed only to satisfy certain regularity conditions; it does not have to be of the form assumed in the model. Posterior inferences for f∗ enable, for example, predictive inference regarding future data points Xm for m > n, inference about functionals of f∗ , or decision-theoretic inferences to minimize a loss of interest. Since posterior asymptotics of the density are not our main focus, here we simply reference the primary results from the literature. 100 For DPMs, there has been great progress in establishing posterior consistency and rates of concentration for the density, in a number of cases. For a general overview, see Ghosal (2010), and for background, see Ghosh and Ramamoorthi (2003). For some of the many contributions to this area, see Ghosal et al. (1999), Barron et al. (1999), Ghosal and Van der Vaart (2001), Lijoi et al. (2005b), Tokdar (2006), Ghosh and Ghosal (2006), Tang and Ghosal (2007), Ghosal and Van der Vaart (2007), Walker et al. (2007), James (2008), Wu and Ghosal (2010), Bhattacharya and Dunson (2010), Khazaei et al. (2012), Scricciolo (2012), Pati et al. (2013), and references therein. More recently, these efforts have been extended to variable-dimension mixtures, and similar results have been obtained. Under fairly general conditions, Kruijer (2008) and Kruijer et al. (2010) have shown that for a large class of variable- dimension location-scale mixtures, the posterior on the density concentrates at the true density at the minimax-optimal rate (up to a logarithmic factor). It should be noted that these optimality proofs apply to a model that is not exactly an MFM as defined in Equation 4.1.1, since there is a single scale parameter shared among all mixture components, rather than a separate scale parameter for each component; the same assumption is made in the DPM optimality proofs. However, this can be viewed as extending the hierarchical model by first sampling the common scale from some prior, and given this, defining the MFM as before. These results suggest that for density estimation, variable-dimension mixtures may perform as well as Dirichlet process mixtures. 101 4.3.1.2 Mixing distribution Now, suppose the true density f∗ is in fact a finite mixture from the assumed family {fθ }, that is, f∗ = fG∗ for some G∗ with s(G∗ ) < ∞. Under this (strong) assumption, MFMs can be shown to be consistent for the mixing distribution (and the density, and the number of components), under very general conditions, by appealing to Doob’s theorem (Nobile, 1994). Roughly speak- ing, Doob’s theorem says that for a correctly-specified, identifiable model Pθ , with probability 1, if the parameter θ0 is drawn from the prior, and the data is i.i.d. from Pθ0 (given θ0 ), then the posterior concentrates at θ0 (Doob, 1949, Ghosh and Ramamoorthi, 2003). The theorem applies under very general conditions, however, a slightly subtle point is that the condition “with probability 1, if the parameter θ0 is drawn from the prior” allows for an exceptional set of prior probability zero, and in an infinite-dimensional space, this exceptional set can be quite large even for seemingly reasonable priors (Freedman, 1963, Diaconis and Freedman, 1986). Fortunately, in an MFM, the parameter space is a countable union of finite- dimensional spaces, since for each k = 1, 2, . . . , a mixture with k components has a finite-dimensional parameter space. Consequently, one can choose a prior such that any set of prior probability zero also has Lebesgue measure zero (where Lebesgue measure is extended in the natural way to this disjoint union of spaces), and indeed, an MFM will have this property whenever p(k) > 0 for all k = 1, 2, . . . and Lebesgue measure is absolutely continuous with respect to the base measure H (for instance, if H has a density with respect to Lebesgue measure that is strictly positive on Θ). Consequently, Doob’s theorem can be applied to prove consistency of MFMs for Lebesgue almost-all true mixing distributions (Nobile, 1994). (A minor techni- 102 cal issue is that Doob’s theorem also requires identifiability, and mixtures are not identifiable in the usual sense; however, as long as the family {fθ } is mixture iden- tifiable, the theorem can be applied by defining a new parameter space containing one representative from each equivalence class.) Regarding DPM consistency for the mixing distribution, Nguyen (2013) has shown that the DPM posterior on the mixing measure G concentrates at the true mixing measure G∗ (in the Wasserstein metric), and has provided a rate of conver- gence. 4.3.1.3 Number of components As mentioned in the preceding section, when the true density f∗ is a finite mixture from the assumed family {fθ }, under very general conditions the MFM posterior on the number of components will concentrate at the true number of components, by Doob’s theorem (Nobile, 1994). (Other approaches have also been proposed for estimating the number of components (Henna, 1985, Keribin, 2000, Leroux, 1992, Ishwaran et al., 2001, James et al., 2001, Henna, 2005), including methods with claims of robustness (Woo and Sriram, 2006, 2007).) On the other hand, DPMs are not well-suited to estimating the number of components in a finite mixture, as discussed at length in Chapters 2 and 3. In Chapter 3, we gave conditions under which the number of clusters in a Gibbs partition mixture model is inconsistent for the number of components, and observed that DPMs satisfy the conditions on the partition distribution. Since MFMs are also Gibbs partition mixtures, the fact that they are consistent for the number of components implies that these conditions are not satisfied for MFMs. This is easy to 103 verify directly, using the following asymptotic expression for Vn (t), which is proved in Section 4.3.4 below: for any t ∈ {1, 2, . . . }, if pK (t) > 0 then t! t! Γ(γt) Vn (t) ∼ (n) pK (t) ∼ pK (t) (4.3.1) (γt) n! nγt−1 as n → ∞. (Here, pK is the distribution of the number of components K, and γ > 0 is the Dirichlet parameter.) Let C be the random partition of [n] as in Section 4.2.1, and let A = (A1 , . . . , AT ) be the ordered partition of [n] obtained by randomly ordering the parts of C, uniformly among the T ! possible choices, where T = |C|. Note that given C, the T ! such choices all yield distinct A’s. Therefore, t p(C) 1 Y 1 Y p(A) = = Vn (t) γ (|c|) = Vn (t) γ (|Ai |) , |C|! t! c∈C t! i=1 where t = |C|. The inconsistency conditions are in terms of a distribution on A of the form p(A) = vn (t) ti=1 wn (|Ai |), and we can put the MFM distribution on A in Q this form by choosing vn (t) = Vn (t)/t! and wn (a) = γ (a) . One of the conditions for inconsistency is that lim supn vn (t)/vn (t + 1) < ∞ for some t. This condition is not satisfied for MFMs — provided that pK (k) > 0 for all k = 1, 2, . . . — since for any t ∈ {1, 2, . . . }, by Equation 4.3.1, vn (t) Vn (t) Γ(γt) pK (t) = (t + 1) ∼ nγ −→ ∞ vn (t + 1) Vn (t + 1) Γ(γ(t + 1)) pK (t + 1) as n → ∞. 104 4.3.2 Relationship between the number of clusters and num- ber of components In the MFM, it is perhaps intuitively clear that, under the prior at least, the number of clusters T = |C| should behave very similarly to the number of components K when n is large. It turns out that under the posterior they also behave very similarly for large n. Assume pK (k) > 0 for all k. We show that for any x1 , x2 , . . . ∈ X , k ∈ {1, 2, . . . }, if lim inf n p(K = k|x1:n ) > 0 then p(T = k | x1:n ) ∼ p(K = k | x1:n ) as n → ∞. See Nobile (2005) for related results. To show this, first we observe that for any t ∈ {1, 2, . . . }, 1 t(t) pn (K = t | T = t) = pK (t) −−−→ 1 (4.3.2) Vn (t) (γt)(n) n→∞ (where pn denotes the MFM distribution with n samples), by Equations 4.2.7 and 4.3.1. Next, write k X p(K = k | x1:n ) = p(K = k | T = t, x1:n ) p(T = t|x1:n ), (4.3.3) t=1 and note that by Equations 4.2.9 and 4.3.2, p(K = k | T = t, x1:n ) = pn (K = k | T = t) → I(k = t) (also see Nobile (2005)). If lim inf n p(K = k | x1:n ) > 0, then dividing Equation 4.3.3 by p(K = k | x1:n ) and taking n → ∞, the terms for t = 1, . . . , k − 1 go to 0, and we are left with 1 = limn p(T = k | x1:n )/p(K = k | x1:n ). 105 4.3.3 Distribution of the cluster sizes under the prior Here, we illustrate one of the major differences between the MFM and DPM priors, namely, the distribution on partitions given the number of clusters/parts. (See Green and Richardson (2001) for similar observations along these lines.) Under the prior, the MFM prefers all parts to be roughly the same order of magnitude, while the DPM prefers having some parts very small. Interestingly, these prior influences remain visible in certain aspects of the posterior, even in the limit as n goes to infinity, as shown by our results on inconsistency of DPMs for the number of components in a finite mixture; see Chapters 2 and 3. In Chapter 5 below, we also empirically observe differences in the posteriors. Let C be the random partition of [n] as in Section 4.2.1, and let A = (A1 , . . . , AT ) be the ordered partition of [n] obtained by randomly ordering the parts of C, uni- formly among the T ! such choices, where T = |C|. Recall that t p(C) 1 Y p(A) = = Vn (t) γ (|Ai |) , |C|! t! i=1 where t = |C|. Now, let S = (S1 , . . . , ST ) be the vector of part sizes of A, that is, Si = |Ai |. Then t t X n! 1 Y n! Y γ (si ) p(S = s) = p(A) = Vn (t) γ (si ) = Vn (t) s1 ! · · · st ! t! i=1 t! i=1 si ! A:S(A)=s for s ∈ ∆t , where ∆t = {s ∈ Zt : P i si = n, si ≥ 1 ∀i} (i.e., the t-part compositions of n). For any x > 0, writing x(m) /m! = Γ(x + m)/(m! Γ(x)) and using Stirling’s approximation, we have x(m) mx−1 ∼ m! Γ(x) 106 as m → ∞. This yields the approximations t Vn (t) n! Y γ−1 p(S = s) ≈ s Γ(γ)t t! i=1 i and t Y p(S = s | T = t) ≈ κ sγ−1 i i=1 for s ∈ ∆t , where κ is a normalization constant. Note that p(s|t), although a discrete distribution, has approximately the same shape as a symmetric t-dimensional Dirichlet distribution. This would be obvious if we were conditioning on the number of components K, and it makes intuitive sense when conditioning on T , since K and T are essentially the same for large n. It is interesting to compare this to the corresponding distributions for the Chinese t Q restaurant process. In the CRP, we have pCRP (C) = αα(n) c∈C (|c|−1)!, and pCRP (A) = pCRP (C)/|C|! as before, so for s ∈ ∆t , t n! 1 αt Y n! αt −1 pCRP (S = s) = (s i − 1)! = s · · · s−1 s1 ! · · · st ! t! α(n) i=1 α(n) t! 1 t and pCRP (S = s | T = t) ∝ s−1 −1 1 · · · st , (4.3.4) which has the same shape as a t-dimensional Dirichlet distribution with all the parameters taken to 0 (noting that this is normalizable since ∆t is finite). Asymp- totically in n, the distribution pCRP (s|t) puts all of its mass in the “corners” of the discrete simplex ∆t , while under the MFM, p(s|t) remains more evenly dispersed. 107 4.3.4 Asymptotics of Vn (t) Recall that (Equation 4.2.2) ∞ X k(t) Vn (t) = pK (k), k=1 (γk)(n) for 1 ≤ t ≤ n, with γ > 0 and pK a p.m.f. on {1, 2, . . . }. Here, we consider the asymptotics of Vn (t): we show that for any t ∈ {1, 2, . . . }, if pK (t) > 0 then t! t! Γ(γt) Vn (t) ∼ (n) pK (t) ∼ pK (t) (4.3.5) (γt) n! nγt−1 as n → ∞. Hence, asymptotically, Vn (t) has a simple interpretation — it behaves like the k = t term in the series. To show this, we use the following elementary result; this is a special case of the dominated convergence theorem, but since the proof is so simple we provide it here. Proposition 4.3.1. For j = 1, 2, . . . , let a1j ≥ a2j ≥ · · · ≥ 0 such that aij → 0 as i → ∞. If ∞ P P∞ j=1 a1j < ∞ then j=1 aij → 0 as i → ∞. P Proof. Letting εR = j>R a1j , we have εR → 0 as R → ∞. For any R > 0, P∞ lim supi j=1 aij ≤ lim supi R P j=1 aij + εR = εR . Taking R → ∞ gives the result. For any x > 0, writing x(n) /n! = Γ(x + n)/(n! Γ(x)) and using Stirling’s approx- imation, we have x(n) nx−1 ∼ n! Γ(x) as n → ∞. Therefore, the k = t term of Vn (t) is t(t) t! Γ(γt) (n) pK (t) ∼ pK (t). (γt) n! nγt−1 108 The first t − 1 terms of Vn (t) are 0, so to establish Equation 4.3.5, we need to show that the rest of the series, divided by the k = t term, goes to 0. (Recall that we have assumed pK (t) > 0.) To this end, let k(t) bnk = (γt)(n) pK (k). (γk)(n) P∞ We must show that k=t+1 bnk → 0 as n → ∞. We apply Proposition 4.3.1 with aij = bt+i,t+j . Clearly, for any k > t, b1k ≥ b2k ≥ · · · ≥ 0. Further, for any k > t, (γt)(n) nγt−1 Γ(γk) ∼ −→ 0 (γk)(n) Γ(γt) nγk−1 as n → ∞, hence, for any k > t, bnk → 0 as n → ∞. Finally, observe that P∞ (n) k=t+1 bnk ≤ (γt) Vn (t) < ∞ for any n ≥ t. Therefore, by Proposition 4.3.1, P∞ k=t+1 bnk → 0 as n → ∞. This proves Equation 4.3.5. 4.4 Inference algorithms The usual approach to inference in variable-dimension mixture models is reversible jump Markov chain Monte Carlo (Richardson and Green, 1997). However, now that we have established that MFMs (as defined in Equation 4.1.1) have many of the same attractive properties as DPMs, much of the extensive body of work on Markov chain Monte Carlo (MCMC) samplers for the DPM can be directly applied to MFMs. One advantage of this is that these samplers tend to be applicable to a wide range of component distribution families {fθ }, without requiring one to design specialized moves. We show how this works for two well-known MCMC sampler algorithms: (1) “Algorithm 3”, for conjugate priors (MacEachern, 1994, Neal, 1992, 2000), and 109 (2) “Algorithm 8”, for non-conjugate priors (Neal, 2000, MacEachern and M¨ uller, 1998). The names Algorithm 3 and Algorithm 8 come from Neal (2000); it has become fairly common to refer to these algorithms by the numbers he gave them. Algorithm 3 is a collapsed sampler (i.e., the component parameters have been integrated out) involv- ing the partition C, while Algorithm 8 is an auxiliary variable algorithm involving C and the component parameters. In fact, we describe a generalization of Algorithm 8, allowing a distribution other than the prior to be used for the auxiliary variables. This allows the algorithm to be used even when one cannot sample from the prior, as well as allowing for the possibility of improved mixing time with a better auxiliary variable distribution. A well-known issue with incremental samplers such as these, when applied to DPMs, is that the mixing time can be somewhat slow, since it may take a long time to create or destroy substantial clusters by moving one element at a time. With MFMs, this issue seems to be exacerbated, since MFMs tend to put small probability (compared with DPMs) on partitions with tiny clusters (see Section 4.3.3), making it difficult for the sampler to move through these regions of the sample space; see Section 5.1.7 for examples of this in experiments. To circumvent this issue, split-merge samplers for DPMs have been proposed in which a large number of elements can be reassigned in a single move (Dahl, 2003, 2005, Jain and Neal, 2004, 2007). In the same way as the incremental samplers, it should be possible to directly apply such split-merge samplers to MFMs, using the properties of the partition distribution described in Section 4.2. More generally, it seems likely that any partition-based MCMC sampler for DPMs could be applied to 110 MFMs as well. (We should point out that split-merge samplers are often used with reversible jump MCMC also (Richardson and Green, 1997, Green and Richardson, 2001); however, as before, there might be an advantage to using the DPM-style split-merge samplers in terms of the ease of application to new models.) 4.4.1 Inference with conjugate priors When H is a conjugate prior for the family of component distributions {fθ : θ ∈ Θ}, often there is an easy-to-compute expression for the “single-cluster marginal likelihood”, Z Y m(xc ) = fθ (xj )H(dθ), j∈c of a subset c ⊂ [n] of the data. To do posterior inference, we would like to draw samples from the posterior on the partition C, Y  p(C|x1:n ) ∝ p(x1:n |C)p(C) = m(xc ) p(C), c∈C (using Equation 4.2.5). The following MCMC algorithm provides an easy way to sample (approximately) from p(C|x1:n ). The algorithm was described by MacEachern (1994), Neal (1992), and Neal (2000) (Algorithm 3) for DPMs; here, we show how it is immediately adapted to MFMs. For ease of comparison, we describe both the MFM and DPM versions, side-by-side. The DPM concentration parameter is denoted by α. Given a partition C, denote by C \ j the partition obtained by removing element j from C. We use restaurant terminology, by way of analogy with the restaurant process described in Section 4.2.6. 111 Algorithm 4.4.1 (Collapsed sampler for MFM and DPM). Initialize C = {[n]} (i.e., a single cluster). Repeat the following N times, to obtain N samples: For j = 1, . . . , n: (1) Choose to reseat customer j . . . MFM DPM m(xc ∪j ) m(xc ∪j ) at table c ∈ C \ j with probability ∝ (|c| + γ) |c| m(xc ) m(xc ) Vn (t + 1) at a new table with probability ∝ γ m(xj ) α m(xj ) Vn (t) where t = |C \ j|. (2) Update C to reflect this reassignment. It is reasonable to assume that m(x1:n ) > 0 (since otherwise something is probably wrong with the model); this ensures that m(xc ) > 0 for all c ⊂ [n]. Further, Vn (t) > 0 if and only if t < t∗ + 1 where t∗ = sup{k : pK (k) > 0} (allowing, possibly, t∗ = ∞); see Section 4.2.4. Hence, the transition probabilities are all well-defined and positive up to t∗ (and there is zero probability of going beyond t∗ ). This algorithm is usually described as a Gibbs sampler in terms of the allocation variables Z1 , . . . , Zn , however, in that formulation one has to accommodate the fact that these variables must take particular values. Our description of the algorithm is directly in terms of the partition C. This can be seen as a blocked Gibbs sampler in terms of the n × n binary matrix in which entry (i, j) is 1 when i and j belong to the same cluster, taking the lth block (l = 1, . . . , n) to be all the entries (i, j) such that i = l or j = l. Another easy way to see that it has the correct invariant 112 distribution (the posterior, p(C|x1:n )) is as a special case of the auxiliary variable algorithm described in Section 4.4.3.2 below; see that section for this argument. Irreducibility is satisfied, since, for example, from any state C 0 with p(C 0 |x1:n ) > 0, there is a positive probability of reaching the state C = {[n]} (i.e., one cluster) in a single sweep, by moving each element to the cluster containing element 1. (Note that this requires Vn (t) > 0 for all t = 1, . . . , |C 0 |, and this is indeed the case; see the discussion above and observe that p(C 0 |x1:n ) > 0 implies Vn (|C 0 |) > 0.) The chain is aperiodic, since there is positive probability of remaining in the same state. 4.4.2 Inference with non-conjugate priors Often, the need arises to use a non-conjugate prior — for instance, we may wish to use a family of component distributions {fθ : θ ∈ Θ} for which there does not exist a conjugate prior. Neal (2000) described a clever auxiliary variable method (Algorithm 8), based on a very similar algorithm of MacEachern and M¨ uller (1998), for inference in DPMs with a non-conjugate prior. The method is particularly attractive since it does not require approximation of the marginal likelihoods m(xc ) to compute the transition probabilities, and consequently, is exact in the sense that its invariant distribution is exactly equal to the posterior, and not an approximation thereof. We show how this algorithm can be adapted to MFMs as well, in the same manner as the previous section. In fact, we describe a generalization of the algorithm allowing for a distribution other than the prior to be used for the auxiliary variables. This permits the algorithm to be used when it is not possible to sample exactly from the prior, and also introduces the possibility of improved mixing time via a better auxiliary variable distribution. We show that when the “single-cluster posterior” 113 (defined below) is used instead of the prior, this reduces — with one minor exception — to the sampler referred to as Algorithm 2 in Neal (2000) (for earlier descriptions of the algorithm, see MacEachern (1994), West et al. (1994), MacEachern (1998)). As before, we present the algorithm in terms of the partition C (instead of the allocation variables Z1 , . . . , Zn as is usually done). Recall from Section 4.2.2 that the model can be equivalently formulated as: C ∼ p(C) iid φc ∼ H for c ∈ C, given C Xj ∼ fφc independently for j ∈ c, c ∈ C, given φ, C. The algorithm constructs a Markov chain with state (C, φ), where φ = (φc : c ∈ C), producing samples from the posterior C, φ|x1:n . (Note that for mathematical convenience, the parameters φc are indexed by subsets c, rather than numbers, but of course any indexing scheme could be used in a software implementation.) The algorithm has two “parameters”: s is the number of auxiliary variables to be used, and (Rx : x ∈ X ) is a family of probability distributions on Θ, parametrized by points x ∈ X in the data space. (On first reading, one can think of the special case in which Rx = H for all x.) Assume H and Rx , for all x ∈ X , have densities h and rx with respect to a common measure ν. Further, assume fθ (x), h(θ), and rx (θ) are strictly positive for all x ∈ X , θ ∈ Θ. We use π(θ|xc ) to denote the density of the single-cluster posterior with respect to ν, Y  π(θ|xc ) = fθ (xj ) h(θ)/m(xc ). j∈c 114 Also, we denote φ \ j = (φ0c0 : c0 ∈ C \ j) where φ0c0 = φc if c0 ⊂ c. Algorithm 4.4.2 (Auxiliary variable sampler for MFM and DPM). Initialize C = {[n]}, and initialize φ[n] ∈ Θ. Repeat the following N times, to obtain N samples: • For each c ∈ C, sample φc ∼ π(φc |xc ), or move φc according to a Markov chain guaranteed to converge to this distribution. • For j = 1, . . . , n: iid (1) If j is seated alone, set η1 ← φ{j} and sample η2 , . . . , ηs ∼ Rxj ; otherwise, iid sample η1 , . . . , ηs ∼ Rxj . (2) Choose to reseat customer j . . . at table c ∈ C \ j with probability ∝ MFM DPM (|c| + γ) fφ0c (xj ) |c| fφ0c (xj ) at a new table with parameter ηi (i = 1, . . . , s) with probability ∝ MFM DPM γ Vn (t + 1) fηi (xj )h(ηi ) α fηi (xj )h(ηi ) s Vn (t) rxj (ηi ) s rxj (ηi ) where t = |C \ j| and φ0 = φ \ j. (3) Update C and φ to reflect this reassignment. To be clear, when reseating customer j, we sample from among |C \ j| + s choices: |C \ j| existing tables or a new table with one of the s parameters η1 , . . . , ηs . Note that for the DPM, the original Algorithm 8 of Neal (2000) is the special 115 case obtained when rx = h for all x ∈ X . Further, we recover a slight variation of Algorithm 2 (Neal, 2000, MacEachern, 1994, West et al., 1994, MacEachern, 1998) (also see Ishwaran and James (2001) for a general formulation) when s = 1 and rx is the single-cluster posterior, rx (θ) = π(θ|x). For, in this case, in step (2), fη1 (xj )h(η1 ) = m(xj ), rxj (η1 ) and thus, for the DPM, the probability of choosing a new table with parameter η1 is proportional to α m(xj ). If j was not seated alone in step (1), then η1 was sampled from π(θ|xj ), and this is equivalent to Algorithm 2; the only difference is that here, if j was seated alone then η1 was set equal to φ{j} , rather than being sampled from π(θ|xj ). This suggests that mixing time could be improved by choosing rx (θ) to approx- imate π(θ|x) (rather than using rx = h). Note that if rx (θ) can be chosen to be exactly π(θ|x), then there seems to be (essentially) no point in choosing s > 1, since it ends up being equivalent to s = 1, except on the occasion that j is seated alone. On the other hand, if rx (θ) is not equal to π(θ|x), then choosing s large can be seen as approximating this choice (as noted by Neal (2000) in the case of rx = h), since then s Z 1 X fηi (xj )h(ηi ) ≈ fη (xj )H(dη) = m(xj ), s i=1 rxj (ηi ) making the probability of choosing a new table approximately the same as in Algo- rithm 2, and given that we choose a new table, choosing among η1 , . . . , ηs propor- tionally to fηi (xj )h(ηi )/rxj (ηi ) amounts to an approximate sample from π(θ|xj ). It should be emphasized, however, that the algorithm has exactly the correct invariant distribution (see Section 4.4.3), and the preceding remarks are simply intended to 116 aid the intuition. From this perspective, the original Algorithm 8 can be interpreted as employing a Monte Carlo approximation to π(θ|xj ) (based on samples from the prior), while the generalization above can be interpreted as using importance sampling. Consequently, if there are many points xj for which π(θ|xj ) is very different than h(θ), and we have a decent approximation to π(θ|x) that is easy to evaluate and sample from, then we would expect the algorithm above to mix significantly better than the original Algorithm 8. In Section 4.4.3, we show the correctness of the algorithm above. In Chapter 5, we apply this algorithm in experiments. 4.4.3 Justification of the non-conjugate sampler In this section, we show that the non-conjugate sampler in Section 4.4.2 produces a Markov chain converging to the posterior distribution C, φ | x1:n . For a Markov chain (Ui ) = (U0 , U1 , U2 , . . . ) in a measurable space (U, A) with a countably-generated sigma-algebra A, it is sufficient to have an invariant distri- bution µ, irreducibility, and aperiodicity to guarantee convergence to the invariant distribution µ, for µ-almost all initial points u0 (Tierney, 1994). If U0 is drawn from a distribution that is absolutely continuous with respect to µ, this guarantees convergence with probability 1. Before establishing these properties, we recall some standard definitions; for general background on Markov chains, we refer to Tierney (1994), Robert and 117 Casella (2004), and Meyn and Tweedie (2009). A transition kernel is a function Q : U × A → [0, 1] such that u 7→ Q(u, A) is measurable for each fixed A ∈ A, and Q(u, ·) is a probability measure on (U, A) for each fixed u ∈ U. The Markov chain (Ui ) is said to have transition kernel Q if P(Ui ∈ A | Ui−1 = u) = Q(u, A). A probability measure µ on (U, A) is invariant with respect to Q if Z µ(A) = µ(du)Q(u, A) (4.4.1) U for all A ∈ A; this is conveniently abbreviated as µ = µQ. For A ∈ A, let T (A) = inf{i > 0 : Ui ∈ A} denote the first time the chain hits A after time 0. Given a sigma-finite measure ϕ on (U, A), the chain (Ui ) is ϕ-irreducible if, for any A ∈ A such that ϕ(A) > 0, we have P(T (A) < ∞ | U0 = u) > 0 for all u ∈ U. The chain is irreducible if it is ϕ-irreducible for some sigma-finite measure ϕ such that ϕ(U) > 0. Assuming Q has invariant distribution µ, it is periodic if there exist disjoint sets A1 , . . . , Ad ∈ A for some d ≥ 2, such that (a) µ(A1 ) > 0, (b) Q(u, Ai+1 ) = 1 for all u ∈ Ai , for i = 1, . . . , d − 1, and (c) Q(u, A1 ) = 1 for all u ∈ Ad . Otherwise, it is aperiodic. 4.4.3.1 Irreducibility and aperiodicity {C} × Θ|C| where the union is over partitions C of [n] such that p(C) > 0. S Let U = Note that the sigma-algebra A associated to U is countably generated, since U is the disjoint union of finitely many spaces of the form Θt , and Θ ⊂ R` has been given the Borel sigma-algebra (which is countably generated). 118 Consider the Markov chain with state U = (C, φ) ∈ U resulting from the non- conjugate sampler; for definiteness, let U0 be the state after the first φ move, and each successive Ui be the state after a complete sweep through the reassignments of j = 1, . . . , n and the following φ move. Let Q be the corresponding transition kernel, and take µ to be the posterior distribution of C, φ | x1:n . Here, we show that the chain is irreducible and aperiodic; in Section 4.4.3.2, we verify that µ is an invariant distribution. To show irreducibility, choose ϕ = δC ∗ ×Π(·|x[n] ), where C ∗ = {[n]} (i.e., one clus- R ter) and Π(·|x[n] ) is the single-cluster posterior, that is, Π(B|x[n] ) = B π(θ|x[n] )ν(dθ). Let A ∈ A such that ϕ(A) > 0, let u ∈ U, and write u = (C 0 , φ0 ). First, observe that similarly to the case of the conjugate prior algorithm, since p(C 0 ) > 0, we have Vn (|C 0 |) > 0, and there is positive probability of reaching C ∗ in a single pass through the reassignments of j = 1, . . . , n, since fθ (x) > 0 for all x ∈ X , θ ∈ Θ by as- sumption. Define B = {θ ∈ Θ : (C ∗ , θ) ∈ A}. When C = C ∗ , the φ move consists of moving φ[n] according to a Markov chain guaranteed to converge to Π(·|x[n] ), so since Π(B|x[n] ) = ϕ(A) > 0, it is guaranteed that φ[n] will hit B in finitely many steps, with probability 1. Since there is positive probability of staying at C ∗ during successive reassignments of j = 1, . . . , n, we have P(T (A) < ∞ | U0 = u) > 0. To show aperiodicity, suppose d ≥ 2 and A1 , . . . , Ad ∈ A are disjoint sets such that µ(A1 ) > 0, Q(u, Ai+1 ) = 1 for all u ∈ Ai , for i = 1, . . . , d − 1, and Q(u, A1 ) = 1 for all u ∈ Ad . Pick C ∗ such that the set B1 = {φ : (C ∗ , φ) ∈ A1 } satisfies µ({C ∗ } × B1 ) > 0, i.e., P(C = C ∗ , φ ∈ B1 | X1:n = x1:n ) > 0. Then P(φ ∈ B1 | C = C ∗ , X1:n = x1:n ) > 0 also. Let Bi = {φ : (C ∗ , φ) ∈ Ai } for i = 2, . . . , d. Then B1 , . . . , Bd are disjoint sets, and since there is a positive probability of staying at C ∗ (with no modifications to φ) during the reassignments of j = 1, . . . , n, the transition kernel K for the φ move must satisfy K(φ, Bi+1 ) = 1 for all φ ∈ Bi , 119 for i = 1, . . . , d − 1, and K(φ, B1 ) = 1 for all φ ∈ Bd . However, this cannot be, since by assumption, K is a transition kernel for a Markov chain that converges to the posterior distribution φ | C ∗ , x1:n , and therefore K must be aperiodic (and have invariant distribution φ | C ∗ , x1:n ) (Tierney, 1994). Hence, a contradiction is obtained, and the chain (Ui ) must be aperiodic. 4.4.3.2 Invariant distribution A surprisingly useful way of constructing a transition kernel with a particular invari- ant distribution µ is via auxiliary variables (Besag and Green, 1993, Edwards and Sokal, 1988): let U ∼ µ and introduce a random variable V that is dependent on U ; if the current state of the chain is ui , then sample V = v | U = ui , and make a move from ui to ui+1 according to a transition kernel Rv for which the conditional distribution U |V = v is invariant. A transition kernel Q constructed in this way (for the overall move from ui to ui+1 ) has µ as an invariant distribution. Often, such a Q will not yield an irreducible chain, but the composition Q1 Q2 · · · Qm of a well-chosen sequence of finitely-many such kernels will yield irreducibility (and will always still have the same invariant distribution). An important special case arises when V = f (U ) for some function f (Tierney, 1994); this can be thought of as partitioning u-space according to the level sets of f , and making a move according to the conditional distribution of the level set containing the current state. The Gibbs sampler can be thought of as a further special case in which U consists of multiple variables and f is a projection onto a subset of variables (typically, projecting onto all but one). The approach of the preceding paragraph can be used to see that the conjugate 120 prior algorithm in Section 4.4.1 has the correct invariant distribution: let U = C, take µ to be the posterior distribution p(C|x1:n ), define fj (C) = C \ j for j = 1, . . . , n, and define the transition matrices Qj (C 0 , C) = p C | fj (C) = fj (C 0 ), X1:n = x1:n  for j = 1, . . . , n. Then, the posterior p(C | x1:n ) is an invariant distribution for each Qj . Applying Q1 Q2 · · · Qn is equivalent to the conjugate prior algorithm described in Section 4.4.1. To show that the non-conjugate prior algorithm in Section 4.4.2 has the correct invariant distribution, we actually use the auxiliary variable technique in two ways: first in a way resembling the preceding paragraph, and again to sample from the conditional distribution. We focus on the MFM case; the DPM case is essentially the same. Let U = (C, φ), take µ to be the posterior distribution of C, φ | x1:n , define fj (C, φ) = (C \ j, φ \ j) for j = 1, . . . , n, and define the transition kernels Qj ((C 0 , φ0 ), A) = P (C, φ) ∈ A | fj (C, φ) = fj (C 0 , φ0 ), X1:n = x1:n .  Clearly, this conditional distribution has a discrete part and a “continuous” part: j can be placed in one of the clusters c ∈ C\j, in which case φ is completely determined, or j can be placed in a separate cluster, in which case φ{j} is free. (Strictly speaking, it is continuous only if H is continuous.) Recall that H has density h with respect to ν. With respect to a measure on U (that is, on the space of (C, φ) pairs) consisting of a unit point mass at these discrete points, and ν on the “continuous” part, the 121 density of the conditional distribution is     fφ0c (xj )(|c| + γ) if c ∈ C \ j and c ∪ j ∈ C   Vn (t + 1) p(C, φ | C \ j, φ \ j, x1:n ) ∝ fη (xj ) γ h(η) if {j} ∈ C and φ{j} = η ∈ Θ   Vn (t)    0 otherwise, where t = |C \ j| and φ0 = φ \ j. (Note that sampling from this distribution would R be easy if we could compute m(xj ) = fη (xj )h(η)ν(dη) and sample from π(η|xj ) ∝ fη (xj )h(η), however, since we are in the non-conjugate case, this is not possible.) In order to construct a move for this conditional distribution, the algorithm uses another auxiliary variable technique, resulting in the probabilities as given in the algorithm; see the next section for details. 4.4.3.3 A simple variable-dimension move Here, we consider — in an abstract setting — a particular auxiliary variable technique for a simple variable-dimension move; then we apply it to the situation at hand. R Let a1 , . . . , at ≥ 0 and g : Y → [0, ∞), Y ⊂ R` , such that κ = Y g(y)ν(dy) < ∞, P where ν is a sigma-finite Borel measure on Y, and define c = i ai + κ. Define X ∈ {0, 1, . . . , t} and Y ∈ Y ∪ {∗}, where ∗ is an arbitrary element not in Y, such that P (X = i) = ai /c for i = 1, . . . , t, P (X = 0) = κ/c, P (Y = ∗ | X > 0) = 1, Y | X = 0 has density g(y)/κ. Suppose we know a1 , . . . , at and can evaluate g(y) for any y ∈ Y, and we would like to sample (X, Y ), but we cannot easily evaluate κ or c. A transition kernel for which 122 the distribution of (X, Y ) is invariant can be constructed using auxiliary variables, as follows. Choose a distribution R on Y that we can easily sample from, and that has a density r with respect to ν that is strictly positive and can be evaluated at any y ∈ Y. Fix s ∈ {1, 2, . . . }. Let M ∼ Uniform{1, . . . , s} independently of X and Y . iid Given X, Y, M : if X > 0 then let Z1 , . . . , Zs ∼ R, and if X = 0 then let ZM = Y iid and Zi ∼ R for i ∈ {1, . . . , s} \ {M }. Denote Z = (Z1 , . . . , Zs ). The following move preserves the distribution of (X, Y ): from state (x, y), sample Z = z | X = x, Y = y, then sample X, Y | Z = z, and discard z. To sample Z|X, Y , we can just sample M, Z|X, Y as described above, and discard M . It turns out that sampling X, Y |Z is also easy, since it can be shown that this is a discrete 1 g(zi ) distribution placing mass ax /C at (x, ∗) for x = 1, . . . , t, and mass at (0, zi ) Cs r(zi ) Pt 1 Ps g(zi ) for i = 1, . . . , s, where C = x=1 ax + s i=1 r(zi ) is the normalizing constant. (Note: if H is not continuous, then multiple zi ’s may coincide, and in this case their mass contributions accumulate.) In fact, from this we can see that, algorithmically, there is actually no need to sample M : by symmetry, the same transition kernel results from always choosing M = 1. We apply this to finish showing that the non-conjugate algorithm has the correct invariant distribution, continuing the discussion from the end of Section 4.4.3.2. Replace a1 , . . . , at by (ac : c ∈ C \ j) where ac = fφ0c (xj )(|c| + γ). Take Y = Θ and Vn (t+1) g(η) = fη (xj ) γ Vn (t) h(η) for η ∈ Θ, and take R = Rxj . Then, the move above is precisely the reassignment move made for each j = 1, . . . , n in Algorithm 4.4.2. Chapter Five Experiments with the MFM 124 In this chapter, we apply the MFM model in various ways, and compare it to the corresponding DPM. Our objectives are to empirically demonstrate: (1) posterior consistency properties, (2) similarities and differences between the MFM and DPM, (3) correctness of the inference algorithms and agreement with previously pub- lished results, and (4) broad applicability and ease of application. To this end, we consider two families of mixtures: univariate normal and bivariate skew-normal. We apply these mixture models to a variety of datasets — simulated and real — and consider the posterior behavior of the models with respect to density estimation, clustering, and the number of components and clusters. We also compare the mixing time of the MCMC samplers for the MFM and DPM. For a number of other experiments with the MFM, see Nobile (1994), Phillips and Smith (1996), Richardson and Green (1997), Stephens (2000), Green and Richardson (2001), and Nobile and Fearnside (2007). 5.1 Univariate normal mixtures Here, we apply the MFM and DPM models in the most standard setting: mixtures of univariate normals. Partly in order to demonstrate agreement with previously published results, some of the experiments are modeled after those of Richardson and Green (1997) and Green and Richardson (2001), who developed reversible jump 125 MCMC samplers for both the MFM and DPM models (in the univariate normal setting) and compared the two models empirically, by considering density estimates, deviance, mean posterior cluster sizes, the entropy of posterior partitions, and the posteriors on the number of components and clusters. Rather than using reversible jump for inference, we use the (non-conjugate) incre- mental Gibbs sampler from Section 4.4.2. In addition to comparing density estimates, mean posterior cluster sizes, and the posteriors on the number of components and clusters, we also consider Hellinger distance to the true density (if known), likelihood on a held-out test set, sample clusterings, and pairwise probability matrices. The characteristics we observe here appear to be true for other families of compo- nent distributions as well; for instance, in Section 5.2, we observe similar properties with mixtures of bivariate skew-normal distributions. 5.1.1 Data We consider the following data distributions and datasets. iid • Standard normal. The data is X1 , . . . , Xn ∼ N (0, 1). This can be interpreted as a mixture with a single component. P4 • Four components. The data is X1 , . . . , Xn i.i.d. from i=1 πi N (x|µi , σi2 ) where π = (0.44, 0.25, 0.3, 0.01), µ = (5, 5, 8, 10), and σ = (1.2, 0.2, 0.6, 0.2). See Figure 5.1. Note that two of the components have the same mean, and one of the components has fairly small weight: 0.01. • Classic galaxy dataset. The “galaxy dataset” (Roeder, 1990) is a standard 126 Figure 5.1: Datasets used, in addition to standard normal. Histograms are shown for the classic galaxy, Shapley galaxy, and SLC datasets, accompanied by rug plots for classic galaxy and SLC. Note: The Shapley dataset also has a small amount of additional data above 45,000 km/s, extending in a long tail up to nearly 80,000 km/s. benchmark dataset for mixture models. It consists of measurements of the velocities of 82 galaxies in the Corona Borealis region. See Figure 5.1. • Shapley galaxy dataset. This is a larger dataset of the same type as the classic galaxy dataset, containing measurements of the velocity of 4215 galaxies in the Shapley supercluster, a large concentration of gravitationally-interacting galaxies (Drinkwater et al., 2004). See Figure 5.1. The clustering tendency of galaxies continues to be a subject of interest in astronomy. Since the as- sumption of normal components appears to be incorrect, this dataset provides the opportunity to see how the MFM and DPM behave on a relatively large amount of real data with a potentially misspecified model. • Sodium-lithium countertransport (SLC) dataset. The “SLC dataset” (Roeder, 1994) consists of the SLC activity level of 190 individuals; it is another bench- mark dataset for mixture models, particularly for estimating the number of 127 components. See Figure 5.1. 5.1.2 Model description To enable comparison, we use the exact same model as in Richardson and Green (1997) and Green and Richardson (2001). The component densities are univariate normal, r −1 λ λ exp − (x − µ)2 ,  fθ (x) = fµ,λ (x) = N (x|µ, λ ) = 2π 2 and the base measure (prior) H on θ = (µ, λ) is µ ∼ N (µ0 , σ02 ), λ ∼ Γ(a, b), independently (where the gamma distribution Γ(a, b) is parametrized to have density ba Γ(a) xa−1 e−bx ). Further, a hyperprior is placed on b, by taking b ∼ Γ(a0 , b0 ). The remaining parameters are set to µ0 = (max{xi } + min{xi })/2, σ0 = max{xi } − min{xi }, a = 2, a0 = 0.2, and b0 = 10/σ02 . Note that the parameters µ0 , σ0 , and b0 are functions of the observed data x1 , . . . , xn . See Richardson and Green (1997) for the rationale behind these parameter choices. (Note: This choice of σ0 may be a bit too large, affecting the posteriors on the number of clusters and components — see Section 5.1.6 — however, we stick with it to enable comparisons to Richardson and Green (1997).) For the MFM, following Richardson and Green (1997), we take K ∼ Uniform{1, . . . , 30} and γ = 1 for the finite-dimensional Dirichlet parameters, as a default. The Shapley galaxy data, however, is highly heterogeneous and K ≤ 30 is too restrictive, 128 so on it we use   c if k ∈ {1, ..., 30}  p(k) =  c/(k − 30)2 if k > 30  where c = 1/(30 + π 2 /6), and γ = 1. For the DPM, we consider two versions: a fixed concentration parameter of α = 1 (DPM-fixed), and a random concentration parameter with α ∼ Exponential(1) (DPM-random). Note that taking µ and λ to be independent results in a non-conjugate prior. Such a prior is appropriate when the location of the data is not informative about the scale (and vice versa). 5.1.3 Approximate inference The hyperprior on b could be handled in at least two different ways. First, integrating b out would give rise to a distribution on (µ, λ) that we could have called H, and in fact, in the case above b could be analytically integrated out since it has been given a conjugate prior. Alternatively, we could use Gibbs sampling (i.e., adding b to the state of the Markov chain, running the sampler given b as usual, and periodically sampling b given everything else). We use the latter approach in order to illustrate that there is no difficulty in handling more general hyperprior structures. Since the prior is non-conjugate, we use the non-conjugate sampler described in Section 4.4.2. For simplicity, we use a single auxiliary variable (s = 1) and use the prior H (given the current value of the hyperparameter b) for the auxiliary variable distributions Rx . For the φ move (that is, moving φc = (µc , λc ) according to the single-cluster posterior π(φc |xc ), for c ∈ C), we use partial conjugacy to sample from π(µc |λc , xc ) and then π(λc |µc , xc ). 129 In contrast to Richardson and Green (1997), we do not restrict the parameter space in any way (e.g., forcing the component means to be ordered to obtain identifi- ability). All of the quantities we consider are invariant to the labeling of the clusters. See Jasra et al. (2005) for an interesting discussion on this point. For the DPM with random α (DPM-random), we numerically integrate out the prior on α in a pre-computation step, as described by MacEachern (1998). For some of the datasets, multiple runs were performed using different sets of samples; for others, a single run was performed with the whole dataset. In each run, the sampler executed 100,000 burn-in sweeps, and 200,000 sample sweeps (for a total of 300,000). Judging by traceplots and running averages of various statistics, this appeared to be sufficient for mixing. The number t = |C| and sizes (|c| : c ∈ C) of the clusters were recorded after each sweep. To reduce memory storage requirements, the full state (C, φ) of the chain was recorded once every 10 sweeps. For each run, the seed of the random number generator was initialized to the same value for both the MFM and DPM. For a dataset of size n, the sampler used for these experiments took approximately 2.5 × 10−6 n seconds per sweep, using a 2.80 GHz processor with 6 GB of RAM. 5.1.4 Density estimation Following Green and Richardson (2001), we use Equation 4.2.15 to estimate the density based on samples from the posterior of C, φ | x1:n . The 20,000 recorded samples (out of the 200,000 sample sweeps) are used. 130 Overall, the MFM and DPM yield highly similar density estimates. For large n, this makes sense, since both models are consistent for the density; the degree of similarity for smaller n as well is interesting. The estimated densities using the DPM with fixed α and random α were very similar, so the density estimation results displayed here are for fixed α, unless mentioned otherwise. 5.1.4.1 Estimated densities Figures 5.2 and 5.3 show the estimated densities for increasing amounts of data (n ∈ {50, 200, 500, 2000}) from the four component distribution and Shapley dataset, respectively; for the Shapley dataset, random subsets of size n were used. Figure 5.4 shows the estimated densities for the classic galaxy, (full) Shapley galaxy, and SLC datasets. 5.1.4.2 Hellinger distance For the simulated data sets (i.e., standard normal and four component), the true density is known, so we can evaluate the performance of a given density estimate by considering its distance from the true density. Hellinger distance is one of the standard metrics used in theoretical studies of the convergence of density estimates, so it makes for a natural choice in this context. The Hellinger distance between densities f and g is  Z 1 p p 2 1/2 H(f, g) = f (x) − g(x) dx . 2 131 Figure 5.2: Estimated densities for MFM (left) and DPM (right) on data from the four component distribution (bottom plot). Results for 5 datasets are displayed for each n ∈ {50, 200, 500, 2000}. For each dataset, the MFM and DPM estimated den- sities are nearly indistinguishable. As n increases, the estimates appear to converge to the true density, as expected. 132 Figure 5.3: Estimated densities for MFM (left) and DPM (right) on velocity data (units: 1000 km/s) from the Shapley galaxy dataset (bottom plot, histogram nor- malized to a density). For each n ∈ {50, 200, 500, 2000}, five random subsets of size n were used. Again, for each set of data, the MFM and DPM estimated densities are nearly indistinguishable. 133 Figure 5.4: Estimated densities using MFM and DPM on the classic galaxy, Shapley galaxy, and SLC datasets. Once again, the MFM and DPM results are visually very similar. For the classic galaxy dataset, we can compare the MFM estimated density to the results of Richardson and Green (1997), and they appear to be similar. 134 Figure 5.5: Estimated Hellinger distances for MFM (red, left) and DPM (blue, right) density estimates of the four component distribution. For each n ∈ {50, 200, 500, 2000}, five independent datasets of size n were used, and the lines connect the averages of the estimated distances. If we can sample from g, and g(x) > 0 whenever f (x) > 0, then we can estimate the squared Hellinger distance by using a simple Monte Carlo approximation, s Z p N 1 p 2 1 X  f (Yi ) 2 H(f, g)2 = f (x) − g(x) dx ≈ −1 2 2N i=1 g(Yi ) iid where Y1 , . . . , YN ∼ g. For the simulated data sets, the true density is simply a mixture of normals, so it is easy to sample from (and it is strictly positive). Figure 5.5 shows the estimated Hellinger distance between the true density and the estimated density of the four component distribution, for increasing amounts of data (n ∈ {50, 200, 500, 2000}). Each distance was estimated with N = 103 independent samples from the true density; the same samples were used to evaluate both MFM and DPM. The MFM and DPM results are very similar. 135 5.1.4.3 Log-likelihood on a test set For a real dataset, of course, we cannot compute the Hellinger distance. In this situation, it is common to evaluate performance by measuring how well the estimated density fits a held-out subset of the data. A standard measure of fit is the log- likelihood on the held-out data. To this end, we compute the density estimate fˆn as before using a “training set” of n points from a dataset, and compute the log-likelihood m X `(fˆn ; y1:m ) = log fˆ(yj ) j=1 where y1 , . . . , ym is a “test set” of m points from the dataset, disjoint from the training set. For this experiment, we used the Shapley galaxy dataset, and also, for com- parison, the four component distribution. Figure 5.6 shows 1 m `(fˆn ; y1:m ) for the MFM and DPM estimated densities, using training sets of increasing size (n ∈ {50, 200, 500, 2000}), and test sets of size m = 1000. For each n, five different training sets (each with a different test set) were used; the same training and test sets were used for the MFM and DPM. The MFM and DPM results are very similar. 5.1.5 Clustering In addition to density estimation, mixture models are very frequently used for cluster- ing. In this section, we examine various properties of the MFM and DPM posteriors that are relevant to clustering. In contrast with density estimation, some significant 136 (a) Four component distribution (b) Shapley galaxy dataset Figure 5.6: Test set log-likelihood (divided by m) for MFM (red, left) and DPM (blue, right) density estimates. For each n ∈ {50, 200, 500, 2000}, five training sets of size n were used, each with a test set of size m = 1000. 137 differences emerge between the MFM and DPM with respect to clustering. It is well-known that partitions, or “clusterings”, C sampled from the DPM pos- terior tend to have small transient clusters (e.g., West et al. (1994), Lartillot and Philippe (2004), Onogi et al. (2011), and others). Here, we empirically observe that MFMs do not exhibit this behavior. When using DPMs for clustering, these small ex- tra clusters can be a nuisance. This makes MFMs a potentially attractive alternative to DPMs for model-based clustering. The differences we observe are what one might expect, based on the rather ex- treme differences in the conditional distribution of the partition a priori, given the number of clusters (as discussed in Section 4.3.3). Roughly speaking, for a given number of clusters, the DPM has a strong preference for lower entropy partitions. In addition to comparing the MFM’s clustering behavior to the DPM with fixed α, we also compare to the DPM with random α. When the data is well-modeled by a finite mixture, it appears that putting a prior on α enables the DPM to adapt somewhat and reduce the tendency to make extra clusters, however, empirically it does seem to still have this tendency to some degree. 5.1.5.1 Sample clusterings First, we look at some samples of clusterings from the posterior. To visualize the clusters, we make a scatter plot in which the x-values are the data points and the y-values are the cluster assignments; the y-values are then slightly vertically jittered in order to distinguish nearby points. Figure 5.7 shows three representative samples for each of the three different mod- 138 els — the MFM, the DPM with random α (DPM-random), and the DPM with fixed α (DPM-fixed) — on n = 500 points from the standard normal distribution (i.e., a mixture with one component). The MFM samples typically have a single cluster. The DPM-random samples usually have a single cluster, but occasionally have one or more small extra clusters. The DPM-fixed samples almost always have multiple small extra clusters. It should be noted that the extra clusters are “transient”, in the sense that they do not consistently contain the same points. Figure 5.8 shows samples from the three models, on n = 2000 points from the four component distribution. All three models fairly consistently have 4 dominant clusters corresponding to the 4 components, accounting for the majority of the data. With respect to extra clusters, we observe the same trend as before: the MFM usually has 4 or 5 clusters, DPM-random usually has 4 to 6, and DPM-fixed usually has 4 to 8. When the MFM does have extra clusters, they are typically substantial in size, whereas the DPM’s extra clusters range in size and usually some are very small. Also, even though one of the true components has small weight (0.01, so it accounts for only around 20 data points out of 2000), all three models clearly recognize it by creating a separate cluster. For the classic galaxy dataset (Figure 5.9), all the models tended to have small clusters. This seems attributable to fact that the data set is fairly heterogeneous and rather small (n = 82). For the Shapley galaxy dataset (Figure 5.10), all the models have a fairly large number of clusters (10 to 20). Nonetheless, as before, the DPM models usually have a few tiny clusters in addition. Interestingly, the behavior seems quite different on the long flat tail at the upper range of values — the MFM often puts all of these values together in one cluster, while the DPMs often break it up into several clusters. 139 Figure 5.7: Typical sample clusterings from the posterior, on standard normal data. See the text for discussion. Also, in this case, DPM-fixed tends to have fewer clusters than DPM-random (see also Section 5.1.6), apparently due to the prior on t. Figure 5.11 shows samples from the three models on the SLC dataset. We observe similar trends as on the one component and four component data. The preceding heuristic observations are based on looking at individual samples. In the subsequent sections, we consider more principled estimates of various posterior properties. 140 Figure 5.8: Typical sample clusterings from the posterior, on four component data. The true component assignments are shown at bottom. See the text for discussion. 141 Figure 5.9: Typical sample clusterings from the posterior, on the classic galaxy dataset. See the text for discussion. 142 Figure 5.10: Typical sample clusterings from the posterior, on the Shapley galaxy dataset. See the text for discussion. 143 Figure 5.11: Typical sample clusterings from the posterior, on the SLC dataset. See the text for discussion. 144 5.1.5.2 Cluster sizes To summarize the relative sizes of the clusters c ∈ C for clusterings C sampled from the posterior, we consider the posterior means of the sorted cluster sizes. In other words, let N1 ≥ N2 ≥ · · · such that (Ni : i ∈ {1, 2, . . . }, Ni > 0) is a permutation of (|c| : c ∈ C), and consider the posterior means E(Ni | x1:n ). Following Green and Richardson (2001), we also consider the mean sizes given the number of clusters t = |C|, that is, E(Ni | t, x1:n ). Figure 5.12 shows estimates of these conditional and unconditional mean sizes for the standard normal distribution with n = 500, the four component distribution with n = 2000, the classic galaxy dataset, and the (full) Shapley galaxy dataset. To visualize these quantities, we draw a box for each i = 1, 2, . . . , stacked up from largest to smallest. The estimates shown are for values of t with sufficiently large posterior probability that we have enough samples to compute the estimates. Note that conditional on t, the DPM with fixed α (DPM-fixed) and with random α (DPM- random) are formally equivalent. Similarly to Green and Richardson (2001), we very consistently observe that under the MFM, the mean cluster sizes given t are “more equally-sized” than under the DPM; more precisely, the discrete distribution corresponding to the sizes (i.e., with probabilities pi = Ni /n) has higher entropy under the MFM. 5.1.5.3 Pairwise probability matrix A useful summary of the posterior distribution of the cluster assignments C is the matrix containing, for each pair of points, the posterior probability that they belong 145 (a) Standard normal distribution, n = 500 (b) Four component distribution, n = 2000 (c) Classic galaxy dataset (d) Shapley galaxy dataset Figure 5.12: Left: Means of the sorted cluster sizes, given t, for MFM (red, left) and DPM (blue, right). Right: Means of the sorted cluster sizes, unconditionally, for MFM (red, left), DPM-random (green, middle), and DPM-fixed (blue, right). 146 to the same cluster. In other words, entry (i, j) is pij = P(∃c ∈ C s.t. i, j ∈ c | x1:n ). This matrix is often used in clustering applications (Medvedovic and Sivaganesan, 2002, Kim et al., 2006, Dahl, 2006, Murino et al., 2010). For each dataset and each model, we estimate this matrix using 5000 samples of C. Figure 5.13 displays estimates of these matrices for the four component distribution (n = 500), the classic galaxy dataset (n = 82), and a random subset of the Shapley galaxy dataset (n = 500). The matrices shown are for the DPM with fixed α; the random α case is nearly identical. The indices of the matrix for the four component distribution are ordered by true component origin, then by the value of the data point; the other matrices are ordered by the value of the data point. The MFM and DPM matrices are visually similar in all three cases, however we observe — particularly for the galaxy dataset — that some of the blocks are darker in the DPM matrices (i.e., these pairs are more likely to be together); this is a real effect, not simply a plotting artifact. This seems to be attributable to the higher entropy of MFM clusterings compared to DPM clusterings. Table 5.1: RMS of the differences between estimated matrices. k = 1 k = 4 galaxy Shapley SLC MFM vs. DPM-fixed 0.09 0.03 0.14 0.04 0.04 MFM vs. DPM-random 0.03 0.03 0.16 0.03 0.06 DPM-fixed vs. DPM-random 0.12 0.02 0.03 0.03 0.04 Table 5.1 displays the root-mean-square (RMS) of the differences between the estimated matrices for the different models, i.e., for matrices (pij ) and (qij ), we 1/2 compute n12 i,j (pij − qij )2 P . For the one component data we use n = 500, and for SLC, n = 190. To get a sense of the accuracy of these estimated matrices, the RMS distance between two estimates of the MFM matrix based on independent sets of 5000 samples was computed for each dataset, and in each case it was approximately 0.01. For the most part, the differences between the models are small, but there are 147 some noticeable disparities: the MFM differs from DPM-fixed and DPM-random on the galaxy dataset, meanwhile, DPM-fixed differs from the others on the standard normal data. 5.1.6 Number of components and clusters It is tempting to use the posterior of the number of components k or of the number of clusters t as a measure of the heterogeneity of the data, however, there are some subtle issues with this that one should bear in mind. (See Aitkin (2001) for a detailed discussion in the context of the galaxy dataset.) In particular, (1) these posteriors can be strongly affected by the base measure H, and (2) these posteriors can be highly sensitive to misspecification of the family of component distributions {fθ }. Issue (1) can be seen, for instance, in the case of normal mixtures: it might seem desirable to choose the prior on the component means to have large variance in order to be less informative, however, this causes the posterior of t to favor smaller values of t. In the MFM, the same is true for the posterior of k (Richardson and Green, 1997, Stephens, 2000, Jasra et al., 2005). With some care, this issue can be dealt with by varying the base measure H and observing the effect on the posterior; for instance, see Richardson and Green (1997). Issue (2) is more serious; in practice, we typically cannot expect our choice of {fθ : θ ∈ Θ} to contain the true component densities (assuming the data is even from a mixture). When the model is misspecified in this way, the posteriors of k 148 (a) Four component distribution, n = 500 (b) Classic galaxy dataset (c) Shapley galaxy dataset (random subset of size n = 500) Figure 5.13: Matrices showing, for each pair of data points, the posterior probability of belonging to the same cluster (white is probability 0, black is 1). 149 and t will typically be severely affected and depend strongly on n. For instance, if the model uses mixtures of normals, and the true data distribution is not a finite mixture of normals, then these posteriors will diverge as n increases. Consequently, the effects of misspecification need to be carefully considered if these posteriors are to be used as measures of heterogeneity. Here, we consider the posteriors of k and t, not for the purpose of measuring heterogeneity, but simply to empirically demonstrate: (1) posterior consistency (and inconsistency), assuming correct specification, (2) differences between the MFM and DPM posteriors, and (3) correctness of the inference algorithms. The posterior of t = |C| is easily estimated from posterior samples of C. Fig- ures 5.14 (standard normal data), 5.15 (four component data), 5.16 (classic galaxy data), 5.17 (Shapley galaxy data), and 5.18 (SLC) show estimates of p(t|x1:n ) for the MFM, DPM-fixed, and DPM-random. For the standard normal, four component, and Shapley galaxy data, for each n ∈ {50, 200, 500, 2000} (as well as n = 5000 for the four component data) we show the average over 5 datasets (using random subsets for the Shapley galaxy data). Under the DPM, the prior on t diverges as n grows, and this influences the posterior. Following Green and Richardson (2001), we can eliminate this influence by “tilting” the posterior in postprocessing; that is, since p(t|x1:n ) ∝ p(x1:n |t)p(t), we can switch to any other prior q(t) on t simply by multiplying the estimated posterior 150 Figure 5.14: Posterior distribution of t on standard normal data. of t by q(t)/p(t) and renormalizing. To this end, define p(t | x1:n ) p∗ (t | x1:n ) ∝ pn (t) where pn (t) is the prior on t when there are n data points; we refer to this as “DPM- flat”. This can be interpreted as corresponding to an (improper) uniform prior on t; in practice it is numerically equivalent to a uniform prior on the range of t values sampled. In the figures mentioned above, we also display p∗ (t|x1:n ) (DPM-flat). For the DPM, the posterior of the number of components k is always trivially a point mass at infinity. For the MFM, to compute the posterior of k, note that ∞ X ∞ X p(k|x1:n ) = p(k|t, x1:n )p(t|x1:n ) = p(k|t)p(t|x1:n ), t=1 t=1 by Equation 4.2.9, and the formula for p(k|t) is given by Equation 4.2.7; therefore, it is easy to transform our estimate of the posterior of t into an estimate of the posterior of k. Figure 5.19 shows the posterior of k alongside that of t for the 151 Figure 5.15: Posterior distribution of t on four component data. Figure 5.16: Posterior distribution of t on the classic galaxy dataset. 152 Figure 5.17: Posterior distribution of t on the Shapley galaxy dataset. Figure 5.18: Posterior distribution of t on the SLC dataset. 153 MFM, for each dataset. We know that the difference between p(t|x1:n ) and p(k|x1:n ) becomes negligible as n grows (see Section 4.3.2); empirically, in these experiments the difference is modest for n = 50 and becomes indistinguishable for larger values of n. Note that the mass of p(k|x1:n ) is always shifted “to the right” (i.e. to higher values), compared to p(t|x1:n ). For the classic galaxy dataset, we can compare the MFM posterior estimated here with the estimate of Richardson and Green (1997); since the exact same model is used, the results should be very similar. Indeed, as Table 5.2 shows, they are. Table 5.2: MFM posterior on k for the classic galaxy dataset. k 1 2 3 4 5 6 7 Here 0.000 0.000 0.067 0.139 0.188 0.194 0.156 R&G 0.000 0.000 0.061 0.128 0.182 0.199 0.160 8 9 10 11 12 13 14 15 0.107 0.066 0.038 0.021 0.011 0.006 0.003 0.002 0.109 0.071 0.040 0.023 0.013 0.006 0.003 0.002 On the finite mixture data, the MFM posteriors on k and t appear to be con- centrating at the true number of components (although perhaps slowly on the four component data), and DPM-fixed certainly does not appear to be concentrating. DPM-random and DPM-flat are somewhere in-between — it is not apparent whether they are concentrating at the true number of components; we would conjecture that they will not concentrate. These empirical observations are consistent with what we know from theory: the MFM is consistent for the number of components (see Section 4.3.1.3) and DPM- fixed is not (see Chapters 2 and 3). Note that the four component data is rather challenging for this task, since one of the components has a weight of only 0.01; this may contribute to the slow 154 (a) Standard normal data (b) Four component data (c) Classic galaxy dataset (n = 82) (d) Shapley galaxy dataset (e) SLC dataset (n = 190) Figure 5.19: Posterior distributions of t (left) and k (right) for the MFM. 155 concentration of the MFM posterior at k = 4. For comparison, in Section 5.2.6, the MFM posterior concentrates much more quickly for data from a three-component mixture of bivariate skew-normals. The SLC dataset is often used as a benchmark for methods estimating the num- ber of components in a finite mixture. The sodium-lithium countertransport (SLC) activity of red blood cells is associated with hypertension, and is thought to have a genetic basis. Assuming a normal model, the objective is to determine whether a “simple dominance model” (two components, with weights p2 + 2pq and q 2 ) or an “additive model” (three components, with weights p2 , 2pq, and q 2 ) is more appro- priate; see (Roeder, 1994) for details. Some methods find more evidence for two components, while others indicate three components — there is no clear consensus (Woo and Sriram, 2006, Roeder, 1994, Ishwaran et al., 2001, Chen et al., 2012). The MFM results shown here (Figures 5.18 and 5.19) are consistent with either two or three components (as are the DPM-random and DPM-flat results), however, as expected, using the DPM-fixed posterior on t would appear to be rather misleading. It is interesting to note that on the Shapley galaxy dataset, DPM-fixed favors a smaller number of clusters than the others; by comparing it with DPM-flat, we can see that this is due to the prior on the number of clusters. 5.1.7 Mixing issues As described in Section 4.4 (and shown empirically in this chapter), it is straight- forward to adapt partition-based DPM samplers to the MFM. An issue that arises with the incremental MCMC samplers is that when n is large, the mixing time for the MFM can be significantly worse than for the DPM. 156 To illustrate, in Figure 5.20, we display examples of traceplots of the number of clusters t for posterior samples using the MFM and DPM, on n = 50 and n = 2000 data points from the four component distribution. For visualization purposes, we uniformly jitter the t values within t ± 0.25, and we display only every 10th sample. Only the burn-in period (105 sweeps) is displayed. DPM-fixed is shown; DPM- random is similar. When n = 50, the mixing is fine for both MFM and DPM, however, when n = 2000 the MFM mixing is substantially worse. Out of the 5 runs with n = 2000, the example shown was a particularly bad run in this respect, but this behavior is not too uncommon, based on previous experiments. The explanation for this seems to be that — as we know from Section 4.3.3 — the DPM likes having small clusters, while the MFM does not. Consequently, the MFM takes a longer time to create or destroy substantial clusters by reassigning one element at a time, since in order to do so it must move through the regions of low probability where there are small clusters. We have observed similar behavior for the conjugate prior algorithm also (Section 4.4.1). While traceplots of t are particularly useful for illustrating this issue, it is worth mentioning that they are a bit misleading with regard to DPM mixing. Since the DPM eagerly creates and destroys small transient clusters, t is significantly more variable than it would be if these tiny clusters were ignored, and consequently, the trace of t indicates a level of mixing that is artificially high. (Adding independent noise to any trace will appear to improve mixing, in an artificial way.) It seems that mixing would be more appropriately assessed by using a statistic that is not strongly affected by tiny clusters, for instance, the entropy of the clustering, − c∈C |c| log |c| P n n , or the sum of the squares of the sizes, c∈C |c|2 . Nonetheless, by considering the P 157 Figure 5.20: Examples of traceplots of the number of clusters t for the first 105 sweeps (the burn-in period), on four component data. 158 traces of various other statistics, it seems that the DPM does indeed mix significantly better than the MFM when n is large. A natural solution would be to use a split-merge sampler (Dahl, 2003, 2005, Jain and Neal, 2004, 2007); although we have not explored this, we expect it to be straightforward to adapt existing DPM split-merge samplers to the MFM. 5.2 Multivariate skew-normal mixtures In order to illustrate the flexibility with which the MFM model can be used, in this section we apply it to mixtures of multivariate skew-normals. This example is fairly general, in the following respects. • Fully non-conjugate. For the univariate normal mixtures, we used a non- conjugate prior, but one in which each parameter separately had a conjugate prior; here, we use a prior with no such conjugacy properties. • Unbounded K. We use a prior on the number of components K under which K has no a priori upper bound (and indeed, pK (k) > 0 for all k = 1, 2, . . . ). • Multivariate. We use a multivariate family of component distributions. Going from univariate to multivariate is routine, due to the convenient form of the samplers. For visualization purposes, we show results for the bivariate case. Overall, the results are quite similar to the case of univariate normal mixtures, so we discuss only a subset of the posterior properties considered in the previous section. 159 Introduced by Azzalini and Dalla Valle (1996) and Azzalini and Capitanio (1999), the multivariate skew-normal distribution SN (ξ, Q, w) with parameters ξ ∈ Rd (lo- cation), Q ∈ Rd×d positive definite (scale and correlation), and w ∈ Rd (skew), has density SN (x | ξ, Q, w) = 2 N (x | ξ, Q) Φ(wT S −1 (x − ξ)) √ for x ∈ Rd , where S ∈ Rd×d is the diagonal matrix with Sii = Qii , and Φ is the univariate standard normal CDF. This is a generalization of the multivariate normal family (which can be recovered by setting w = 0) which retains some of its nice properties. Intuitively, it is a multivariate normal which has been skewed, according to the skew parameter w ∈ Rd , in the direction of w/|w| by a magnitude of |w|. Figure 5.21 shows a contour plot and a scatterplot of n = 2000 points from SN (ξ, Q, w) with ξ = (−1.5, −1)T , Q = 6.25 1.5  1.5 1 , and w = (1, 3)T . A convenient  2  σ ρσ σ parametrization of Q in the bivariate case (d = 2), is as Q = ρσ11σ2 σ12 2 , where 2 σ1 , σ2 > 0 and ρ ∈ (−1, 1) is the correlation coefficient; for instance, in the preceding example, σ1 = 2.5, σ2 = 1, and ρ = 0.6. Due to its increased flexibility over the multivariate normal, the skew-normal may be useful for applications involving mixture models, particularly clustering. 5.2.1 Data We consider data simulated from a three-component mixture of skew-normals P3 i=1 πi SN (x|ξi , Qi , wi ), with parameters as shown in Table 5.3. Figure 5.22 displays contour plots and samples from each of the three mixture components. 160 Figure 5.21: Contour plot and scatterplot of a bivariate skew-normal distribution. Table 5.3: Parameters of the skew-normal mixture used for data simulation i πi ξi1 ξi2 σi1 σi2 ρi wi1 wi1 1 0.45 2 −2 2.5 1 −0.6 −3 6 2 0.3 0 0 2 1 0.95 6 0 3 0.25 −1 2 1 1 0 0 −2 161 Figure 5.22: Mixture components used for simulations (top: contour plots, bottom: scatterplots). The scatterplots shows 450, 300, and 250 samples from the three components, respectively. 162 5.2.2 Model description We take the family of component densities to be bivariate skew-normal, fθ (x) = SN (x | ξ, Q, w), with the parametrization θ = (ξ1 , ξ2 , log σ1 , log σ2 , logit ρ+1 2 , w1 , w2 ) (where logit p = p log 1−p ). For the base measure H, for convenience we simply choose θ ∼ N (0, C) where C 1/2 = diag(5, 5, 2, 2, 2, 4, 4). Note that this parametrization covers all of R7 . This parametrization and choice of H was chosen rather arbitrarily, without any particular regard to the structure of the skew-normal distribution; a more principled choice of H is almost certainly possible. The location and scale of H was chosen to be roughly appropriate for the data at hand. For the MFM parameters, we use γ = 1 and   c if k ∈ {1, ..., 30}  p(k) =  c/(k − 30)2 if k > 30  where c = 1/(30+π 2 /6). For the DPM, we use only a fixed concentration parameter, α = 1. 5.2.3 Approximate inference Since H is not a conjugate prior for {fθ }, we again use the non-conjugate sam- pler (Section 4.4.2). For the auxiliary variable distribution, we again use H. For 163 the φ move (that is, sampling the component parameters), we use the Metropolis algorithm with a normal proposal distribution. Specifically, for c ∈ C, we pro- iid pose φ0c = φc + 0.1 ε where ε1 , . . . , ε7 ∼ N (0, 1) (and accept φ0c with probability min{1, π(φ0c |xc )/π(φc |xc )}); this is repeated 20 times between each pass through the reassignments of j = 1, . . . , n. Note that the acceptance probability is ( ) h(φ0c ) j∈c fφ0c (xj ) Q min 1, Q , h(φc ) j∈c fφc (xj ) which does not actually require computing the posterior density. Although a small amount of experimentation was done with these parameters, this move could al- most certainly be improved upon (e.g., perhaps by scaling the proposal distribution separately for each dimension). For each n ∈ {50, 200, 500, 2000}, five independent sets of data from the true distribution were considered. As before, for each such set, the sampler was run for 100,000 burn-in sweeps and 200,000 samples sweeps (for a total of 300,000); this appeared to be sufficient for mixing. For a dataset of size n, the sampler used for these experiments took approximately 1.5 × 10−5 n seconds per sweep, using a 2.80 GHz processor with 6 GB of RAM. 5.2.4 Density estimation As before, Equation 4.2.15 was employed to estimate the density based on samples from the posterior of C, φ | x1:n . The 20,000 recorded samples (out of the 200,000 sample sweeps) were used. Figure 5.23 shows representative examples of the esti- mated densities for increasing amounts of data (n ∈ {50, 200, 500}). Although slight 164 differences between the MFM and DPM can be detected, overall, the densities ap- pear to be quite similar. The wiggliness of the contour lines for the DPM is probably due to the presence of tiny clusters in the posterior samples; see below. Since the true density is known, we can assess the performance of the density estimates using Hellinger distance (see Section 5.1.4). Figure 5.24 shows the esti- mated Hellinger distance between the true density and the estimated densities, for increasing amounts of data. Each distance was estimated with N = 103 independent samples from the true density; the same samples were used to evaluate both MFM and DPM. As usual, the MFM and DPM results are very similar. 5.2.5 Clustering We display some representative clusterings sampled from the posterior. Figure 5.25 shows (Top) three samples from the MFM and the DPM, and (Bottom) the true component assignments of the points, for n = 500 data points. The markers for points in small clusters have been enlarged for visualization purposes. All of the data points are visible within the window shown. Similarly to before, we observe that both models consistently have three dominant clusters corresponding to the three true components, and that the DPM tends to have small transient extra clusters, while the MFM does not. The behavior of the means of the sorted cluster sizes is similar to the univariate normal case. Also, the pairwise probability matrices for the MFM and DPM are very similar to one another. 165 Figure 5.23: Contour plots of typical density estimates for the MFM (left) and DPM (right); compare with the true density (bottom). 166 Figure 5.24: Estimated Hellinger distances for MFM (red, left) and DPM (blue, right) density estimates. For each n ∈ {50, 200, 500, 2000}, five independent datasets of size n were used, and the lines connect the averages of the estimated distances. 5.2.6 Number of components and clusters We estimate the posteriors on the number of components and clusters, as described in Section 5.1.6. As explained there, we do not consider these for the purpose of measuring heterogeneity, but simply to demonstrate posterior consistency properties and differences between the MFM and DPM. Figures 5.26 shows estimates of p(t|x1:n ) and p(k|x1:n ) for the MFM, p(t|x1:n ) for the DPM (DPM-fixed), and p∗ (t|x1:n ) for the DPM (DPM-flat). For each n ∈ {50, 200, 500, 2000}, we show the average over 5 datasets. As before, the MFM appears to be concentrating at the true number of components, the DPM does not appear to be, and for DPM-flat it is not clear (although it does seem to be “stalled”, suggesting that it might not concentrate). This empirical observation for the MFM is actually rather interesting, since skew- normal mixtures are not mixture identifiable in general, and therefore we have no 167 Figure 5.25: Typical sample clusterings from the posterior; see the text. 168 Figure 5.26: Posterior distributions of k and t. guarantee that the posterior will concentrate at the true number of components. This lack of identifiability can be easily exhibited (as pointed out to me by Dan Klein) by observing that 1 1 SN (x | ξ, Q, w) + SN (x | ξ, Q, −w) = SN (x | ξ, Q, 0) 2 2 for any x, ξ, Q, w, since Φ(z) + Φ(−z) = 2Φ(0) for any z ∈ R. Chapter Six Combinatorial stochastic processes for other variable-dimension models 170 6.1 Introduction Like the Dirichlet process mixture (DPM), a number of other nonparametric models are derived as a certain infinite-dimensional limit of a family of finite-dimensional models. For each such model, we could alternatively construct a variable-dimension model by putting a prior on the dimension, in the same way that we constructed the MFM. As we have seen with the MFM, a variable-dimension model can have many of the same characteristics as an infinite-dimensional model, as well as having certain advantages. Two key observations underlying this development were that the MFM has a nice partition distribution, and that this distribution can be generated by a simple combinatorial stochastic process (a restaurant process, in this case). In this chapter, we show that for two other popular nonparametric models — the hierarchical Dirichlet process (HDP) and the Indian buffet process (IBP) — the natural variable-dimensional counterparts also have many of the same attractive properties as the nonparametric, infinite-dimensional models. In particular, they give rise to nice distributions on discrete structures that can be generated by simple combinatorial stochastic processes (a “franchise process” and a “buffet process”, respectively). 6.2 Hierarchical mixture of finite mixtures In many applications, the data comes in groups that may be related in some ways but may also have significant differences. For instance, in modeling a set of documents, 171 the words in each document may be considered as a group of data points, and documents may be related by sharing various topics. One way to model such data would be to combine the data from all groups and use a single mixture model, however, this does not allow for possible differences between the groups. At the other extreme, one could use a separate mixture model for each group, but this does not take advantage of possible relationships between the groups. A third way — the one we pursue here — is to use a mixture model for each group and allow mixture components to be shared among groups. The hierarchical Dirichlet process (HDP) (Teh et al., 2004, 2006) is a model of this type, in which each group is modeled using a Dirichlet process mixture for which the base measure is itself a draw from a Dirichlet process, and is shared among all groups. Due to the discreteness of random measures drawn from a DP, this allows subsets of the component parameters of these different mixtures to be shared. The HDP is a special case of the Analysis of Densities model of Tomlinson and Escobar (1999), and can also be viewed as an example of a dependent Dirichlet process (MacEachern et al., 1999, MacEachern, 2000, MacEachern et al., 2001). The HDP has seen a wide range of applications, including document modeling (Teh et al., 2004), natural language modeling (Liang et al., 2007), object tracking (Fox et al., 2007), haplotype inference (Xing et al., 2006), natural image processing (Kivinen et al., 2007), cognitive science (Griffiths et al., 2007), and measuring similarity between musical pieces (Hoffman et al., 2008), and a dynamic hierarchical Dirichlet process (dHDP) has been proposed by Ren et al. (2008). (For other nonparametric approaches to modeling data in multiple related groups, see M¨ uller et al. (2004), Dunson (2006), and Rodriguez et al. (2008).) In the same way that the HDP is constructed from the Dirichlet process, one can construct a hierarchical mixture model using the corresponding distribution on 172 discrete measures for the MFM. We show that the resulting model, which we refer to as a hierarchical mixture of finite mixtures (HMFM), has some of the same attractive properties as the HDP mixture model. In particular, it has the same exchangeabil- ity properties, it can be represented in terms of a “franchise process” that closely parallels the Chinese restaurant franchise process of the HDP, and MCMC inference algorithms for the HDP that are based on the franchise process representation can be adapted to the HMFM. 6.2.1 A hierarchical variable-dimension mixture model First, we give a succinct description of the model, employing the distribution M(pK , γ, H) on discrete measures G associated with the MFM; see Section 4.2.7. As before, we use fG to denote the mixture obtained from the discrete measure G, that R is, fG (x) = fθ (x)G(dθ). Assume we have a family {fθ : θ ∈ Θ} and base measure H on Θ just as before. Suppose we have data in m groups, with nr data points in the rth group. For instance, in the example of documents mentioned above, we have m documents with nr words in document r. The hierarchical mixture of finite mixtures (HMFM) model is as follows: G0 ∼ M(q0 , γ0 , H) For r ∈ {1, . . . , m} independently, given G0 : Gr ∼ M(qr , γr , G0 ) (6.2.1) Xrj ∼ fGr (x) independently for j ∈ {1, . . . , nr }, where q0 , q1 , . . . , qm are p.m.f.’s on {1, 2, . . . } and γ0 , γ1 , . . . , γm > 0. See Figure 173 G0 Kr K0 Gr πr Wr π0 Xrj Zrj Xrj θ nr nr m m (a) (b) Figure 6.1: Graphical models for (a) the random discrete measure formulation, and (b) the equivalent representation including additional latent variables. 6.1(a) for a graphical model representation. The only difference between this and an HDP mixture model is that in the HDP, the M distributions are replaced by Dirichlet processes: specifically, M(q0 , γ0 , H) is replaced by DP(α0 , H), and M(qr , γr , G0 ) is replaced by DP(αr , G0 ) for r = 1, . . . , m. This model has the same exchangeability properties as the HDP mixture model: the datapoints within each group are exchangeable (being conditionally i.i.d.), and if q1 = · · · = qm and γ1 = · · · = γm , then the order of the groups does not matter, as long as the group sizes n1 , . . . , nm are handled appropriately; for instance, the distribution is the same if the sizes are permuted, the data is sampled, and then the groups are permuted back to the original order of the sizes (or alternatively, if the sizes are taken to be infinite, then the groups are exchangeable). This model can also be extended to multi-level hierarchies in the same way as the HDP. It is also useful to represent the model in an equivalent way, using a number of additional latent variables, as follows. Despite the flurry of symbols, this is a very 174 natural construction. See Figure 6.1(b) for a graphical model representation. K0 ∼ q0 π0 ∼ Dirichletk0 (γ0 , . . . , γ0 ), given K0 = k0 iid θ1 , θ2 , . . . ∼ H For r = 1, . . . , m: (6.2.2) Kr ∼ qr πr ∼ Dirichletkr (γr , . . . , γr ), given Kr = kr iid Wr1 , Wr2 , . . . ∼ π0 , given π0 iid Zr1 , . . . , Zrnr ∼ πr , given πr Xrj ∼ fθYrj , where Yrj = WrZrj , indep. for j ∈ {1, . . . , nr }, given (θi ), Wr , Zr . The representation in terms of discrete measures can be recovered by setting k0 X kr X G0 = π0i δθi and Gr = πri δθWri i=1 i=1 for r = 1, . . . , m, and the distribution on the data (Xrj ) is the same in both of these representations. 6.2.2 Franchise process In the same way as an HDP mixture, this model can also be described in terms of a simple “franchise process” and an associated distribution on combinatorial struc- tures. We use the same restaurant/franchise analogy as the HDP. We envision m different restaurants all belonging to the same franchise, and imagine the m groups of datapoints as groups of customers visiting these restaurants, respectively. Label the 175 jth customer at restaurant r by (r, j). For each restaurant r = 1, . . . , m, we partition the customers according to which table they are seated at, resulting in a partition Cr of {(r, 1), (r, 2), . . . , (r, nr )}. There are tr := |Cr | occupied tables at restaurant r, and L := m P r=1 tr occupied tables altogether. Now, we imagine that each occupied table is served a single dish, chosen from a franchise-wide menu. We partition the L occupied tables according to which dish they are served, and represent this as a partition C0 of m S r=1 Cr . (In this notation, C0 is a set of sets of sets of customers.) Define ∞ X k(t) Vnr (t) = qr (k) k=1 (γr k)(n) for r = 0, 1, . . . , m and 0 ≤ t ≤ n. Consider the following process for generating C1 , . . . , Cr and C0 : • Initialize C0 ← {} and L ← 0. • For each restaurant r = 1, . . . , m: – Initialize Cr ← {}. – For j = 1, . . . , nr : (1) The jth customer at restaurant r sits . . . ∗ at an existing table c ∈ Cr with probability ∝ |c| + γr Vjr (tr + 1) ∗ at a new table with probability ∝ γr Vjr (tr ) where tr = |Cr | is the number of tables occupied so far at restaurant r. (2) If a new table is chosen in step (1), set L ← L + 1 and select a dish for it from the franchise-wide menu, choosing . . . ∗ an existing “dish” d ∈ C0 with probability ∝ |d| + γ0 (note that |d| is the total number of tables having this dish so far, in all restaurants) 176 VL0 (t0 + 1) ∗ a new dish with probability ∝ γ0 VL0 (t0 ) where t0 = |C0 | is the number of dishes tried so far in all restaurants. (3) Update Cr and C0 to reflect these choices. Note that the first customer at each restaurant always sits at a new table since there are no existing tables, and similarly, the first table in the whole process always gets a new dish. The only difference between this and the Chinese restaurant franchise process of the HDP is that for the HDP, in step (1) table c is chosen with probability ∝ |c| or a new table is chosen with probability ∝ αr , and in step (2) dish d is chosen with probability ∝ |d| or a new dish is chosen with probability ∝ α0 . For each r ∈ {1, . . . , m}, the seating of customers in restaurant r simply follows the MFM restaurant process (Section 4.2.6), independently of the other restaurants. And given the seating arrangements C1 , . . . , Cm , the serving of dishes also follows the MFM restaurant process, this time with occupied tables playing the role of customers. Note that it is not necessary to immediately choose a dish for each new table; in particular, it would be equivalent if we waited until all the customers at all the restaurants had been seated before choosing dishes. Therefore, the probability of obtaining C0:m := (C0 , C1 , . . . , Cm ) is m Y p(C0:m ) = PL0 (C0 ) Pnrr (Cr ) (6.2.3) r=1 where Y Pnr (C) = Vnr (|C|) γr(|c|) (6.2.4) c∈C 177 Pm is the MFM distribution on partitions (Equation 4.2.1), and L = r=1 |Cr |. Note that C0 depends on Cr through L. (Of course, p(C0:m ) is as above only when C0:m is valid, and p(C0:m ) = 0 otherwise; C0:m is valid when Cr is a partition of {(r, 1), (r, 2), . . . , (r, nr )} for r = 1, . . . , m and C0 is a partition of m S r=1 Cr .) For comparison, the formula for the probability of C0:m under the Chinese restau- rant franchise process is the same, except that Pnr (C) is replaced by the DP distri- |C| (n) Q bution on partitions, (αr /αr ) c∈C (|c| − 1)!. By exchangeability of the restaurant process, the distribution of C0:m is invariant to the order in which the customers enter each restaurant, and further, it is not nec- essary for all the customers at restaurant r to be seated and served before customers at r0 > r start sitting and being served — the same distribution of C0:m will result if steps (1), (2), and (3) of the process are followed for any temporal ordering of the customers. We can use this distribution on combinatorial structures C0:m to construct a hierarchical mixture model, as follows: Cr ∼ Pnrr independently for r = 1, . . . , m C0 ∼ PL0 given C1:m , where L = m P r=1 |Cr | (6.2.5) iid φd ∼ H for d ∈ C0 , given C0 Xrj ∼ fφd independently for (r, j) ∈ c, c ∈ d, d ∈ C0 , given C0 . See Figure 6.2 for a graphical model representation. We claim that this is equivalent to the models described above; see the next section for this argument. 178 C0 φ Cr Xrj nr m Figure 6.2: Graphical model for combinatorial representation of the model. 6.2.3 Equivalence of the models Here, we show that the combinatorial model in Equation 6.2.5 is equivalent to the model in Equation 6.2.2, in the sense that it simply represents different variables in the latent structure, and gives rise to the same distribution for the data (Xrj ). To show this, we will start with the model as described in Equation 6.2.2, then define C0 , C1 , . . . , Cm and (φd : d ∈ C0 ) in this context, and finally, observe that we obtain the same distributions as in Equation 6.2.5. For r = 1, . . . , m, let Cr be the partition of {(r, 1), . . . , (r, nr )} induced by Zr1 , . . . , Zrnr . Note that C1 , . . . , Cm are independent, and by construction, Cr simply has the MFM distribution on partitions, Pnrr , in the notation of Equation 6.2.4 above. Sm Now, let T = r=1 Cr . Recalling that Yrj = WrZrj , note that for any c ∈ T all of the elements (r, j) ∈ c have the same value of Yrj ; denote this value by Yc . Let C0 be the partition of T induced by (Yc : c ∈ T ). Given C1:m and π0 , observe that (Yc : c ∈ T ) are i.i.d. from π0 (by applying the argument of Section 4.2.2 to each r). Therefore, given C1:m (marginalizing out π0 ), the partition induced by (Yc : c ∈ T ) — namely C0 — simply has an MFM partition distribution; specifically, C0 | C1:m has p.m.f. PL0 in the notation of Equation 6.2.4, where L = |T | = m P r=1 |Cr |. 179 For d ∈ C0 , let φd equal the θi such that Yc = i for all c ∈ d. Again, by the same argument as in Section 4.2.2, (φd : d ∈ C0 ) are i.i.d. from H, given C0 . Finally, it is easily verified that the distribution of (Xrj ) | C0 , (φd ) is Xrj ∼ fφd independently for (r, j) ∈ c, c ∈ d, d ∈ C0 . 6.2.4 Inference Since the HMFM model has a franchise process which closely parallels that of the HDP mixture model, the Gibbs sampling algorithm for the HDP mixture (Teh et al., 2004) can be easily adapted to the HMFM. It might be interesting to see if the other inference algorithms proposed for the HDP mixture, such as the augmented sampler of Teh et al. (2006), can also be adapted to the HMFM. 6.3 Mixture of finite feature models In a sense, a mixture model can be viewed as having a single discrete latent feature, controlling which component each data point comes from. Often, it is more realistic to allow for multiple latent features, enabling the model to parsimoniously represent more complex relationships. The latent features are sometimes chosen to be binary- valued (zero or one), and in this case, the feature values can be represented in a binary matrix Z such that Zij = 1 when item i has feature j. Given such a matrix Z, the observed data X may be modeled in a variety of ways. A simple example is to take X Xi = Zij µj + εi j 180 where the µj ’s and εi ’s are multivariate normal; this can be viewed as a binary h i factor analysis model, where M := µ1 · · · µk is the “factor loading matrix” and Zi1 , . . . , Zik are the “factors” for item i. The Indian buffet process (IBP) of Griffiths and Ghahramani (2005) is a non- parametric model for Z with infinitely many features, finitely many of which are possessed by any given item. The IBP has seen a number of applications, including models for protein complexes (Krause et al., 2006), gene expression data (Knowles and Ghahramani, 2007), causal graph structure (Wood et al., 2006), dyadic data (Meeds et al., 2007), similarity judgments (Navarro and Griffiths, 2007), network data (Miller et al., 2009), and multiple time-series (Fox et al., 2009). The IBP has some interesting theoretical properties (Thibaux and Jordan, 2007, Teh et al., 2007), and has inspired the development of a number of other models. The IBP can be derived as a certain infinite-dimensional limit of a family of finite feature models. Here, we show that by instead placing a prior on the number of features — that is, using a variable-dimension feature model — one obtains a distribution with many of the same attractive properties as the IBP: an exchangeable distribution on equivalence classes of binary matrices, representation via a simple buffet process, and approximate inference via the same Gibbs sampling algorithms. In fact, the IBP (as well as the two-parameter generalization of Ghahramani et al. (2007)) can be viewed as a limiting case, corresponding to a certain corner of the parameter space. For further background on the IBP, including motivation for the use of such priors in latent feature models, and several applications, we refer to Ghahramani et al. (2007) and references therein (especially Roweis and Ghahramani (1999)). Variable-dimension feature models have previously been developed by Lopes and 181 West (2004), Dunson (2006), and Ghosh and Dunson (2009), in the context of factor analysis, however, in such models the factors are continuous rather than binary- valued. For inference, Lopes and West (2004) used reversible jump Markov chain Monte Carlo, while Dunson (2006) and Ghosh and Dunson (2009) used a parameter- expanded sampler. Sparse variants of these models, in which many zeros appear in the entries of the factor loading matrix (rather than in the factors themselves), have been applied to gene expression profiling (West, 2003, Carvalho et al., 2008). 6.3.1 A distribution on binary matrices Let p(k) be a p.m.f. on {0, 1, 2, . . . }, and let a, b > 0. Let K ∼ p(k) iid π1 , . . . , πk ∼ Beta(a, b), given K = k (6.3.1) Zij ∼ Bernoulli(πj ) independently for i ∈ {1, . . . , n}, j ∈ {1, . . . , k}, given K = k and π1:k . We refer to this as a mixture of finite feature models (MFFM). When K is fixed, say K = k, this is precisely the finite feature model described by Griffiths and Ghahramani (2005), forming the basis for the IBP, which is then obtained by setting b = 1, a = α/k, and taking k → ∞. 6.3.2 A simple urn process In studying the properties of this model, it is useful to first consider a single column, iid leading to a simple urn process. Let π ∼ Beta(a, b), and U1 , . . . , Un ∼ Bernoulli(π) 182 given π. Then B(a + s, b + n − s) a(s) b(n−s) Z p(u1:n ) = p(u1:n |π)p(π)dπ = = (6.3.2) B(a, b) (a + b)(n) ui , and x(m) = P where B(a, b) = Γ(a)Γ(b)/Γ(a + b) is the beta function, s = i x(x + 1) · · · (x + m − 1) with the convention that x(0) = 1. Since this quantity appears so frequently, it is useful to introduce the notation a(s) b(n−s) wn (s) = . (6.3.3) (a + b)(n) From Equation 6.3.2, it is easy to see that U1 , . . . , Un can also be generated by taking the first n draws in the following urn process:  a  U1 ∼ Bernoulli a+b  a + Pm−1 U  i=1 i Um | Um−1 , . . . , U1 ∼ Bernoulli for m = 2, 3, . . . . a+b+m−1 From this perspective, a and b can be interpreted as “pseudo-counts”. Note that if  M = min i ∈ {1, 2, . . . } : Ui = 1 , then for any m ∈ {1, 2, . . . }, ab(m−1) P(M = m) = = wm (1), and (6.3.4) (a + b)(m) b(m) P(M > m) = = wm (0), (6.3.5) (a + b)(m) by Equation 6.3.2. 183 6.3.3 Equivalence classes of matrices 6.3.3.1 Removing columns of all zeros For any binary matrix z ∈ {0, 1}n×k , define z ∗ be the matrix obtained from z by removing any columns consisting of all zeros (and preserving the order of the re- maining columns). Define sj (z) = ni=1 zij . Let Z be as in Equation 6.3.1. Then, P as we will see, for z ∗ ∈ {0, 1}n×t , t ∈ {0, 1, 2, . . . }, such that s1 (z ∗ ), . . . , st (z ∗ ) > 0, t Y ∗ ∗ P(Z = z ) = vn (t) wn (sj (z ∗ )) (6.3.6) j=1 where wn (s) is as in Equation 6.3.3 above, and ∞   X k vn (t) = wn (0)k−t p(k). (6.3.7) k=0 t k  By convention, an empty product equals 1 and t = 0 when k < t. For comparison, under the IBP the corresponding distribution is t ∗ ∗ αt Y (sj (z ∗ ) − 1)! (n − sj (z ∗ ))! PIBP (Z = z ) = exp(−αHn ) (6.3.8) t! j=1 n! Pn where Hn = k=1 1/k is the nth harmonic number. To derive Equation 6.3.6, first note that by Equation 6.3.2, for z ∈ {0, 1}n×k , k Y t Y P(Z = z | k) = wn (sj (z)) = wn (0) k−t wn (sj (z ∗ )), (6.3.9) j=1 j=1 184 where t is the number of nonzero columns in z. Then, X P(Z ∗ = z ∗ | k) = P(Z = z0 | k) I(z0∗ = z ∗ ) z0 ∈{0,1}n×k n o t Y ∗ ∗ = # z0 ∈ {0, 1} n×k : z0 = z wn (0)k−t wn (sj (z ∗ )) j=1   t k Y = wn (0)k−t wn (sj (z ∗ )). (6.3.10) t j=1 P∞ Writing P(Z ∗ = z ∗ ) = k=0 P(Z ∗ = z ∗ | k)p(k) yields Equation 6.3.6. 6.3.3.2 Staircase form Let mj (z) = min{i : zij = 1 or i > n}. For z ∈ {0, 1}n×k , define z 0 = [B1 B2 · · · Bn ] where Bi is the submatrix of z consisting of the columns j such that mj (z) = i. In other words, z 0 is the matrix obtained by permuting the columns of z ∗ so that m1 (z 0 ) ≤ · · · ≤ mt (z 0 ), where t is the number of nonzero columns in z, while preserving the order of the columns within each block {j : mj (z ∗ ) = i}. We refer to this as staircase form. Then, the distribution of Z 0 is as follows: for z 0 ∈ {0, 1}n×t , t ∈ {0, 1, 2, . . . }, such that m1 (z 0 ) ≤ · · · ≤ mt (z 0 ) ≤ n, t 0 0 t! Y P(Z = z ) = vn (t) wn (sj (z 0 )) (6.3.11) t1 ! · · · tn ! j=1 where ti = #{j : mj (z 0 ) = i}. This is easily derived from Equation 6.3.6, since there are t!/(t1 ! · · · tn !) matrices z ∗ giving rise to z 0 , each of which has probability vn (t) tj=1 wn (sj (z 0 )). Q 185 Likewise, PIBP (Z 0 = z 0 ) is given by Equation 6.3.8 times t!/(t1 ! · · · tn !); for refer- ence, this is Equation (5) of Griffiths and Ghahramani (2005). A case of particular interest arises when p(k) = Poisson(k|λ). In this case, we have the closed-form expression vn (t) = exp(λwn (0)) Poisson(t|λ), (6.3.12) by a straightforward calculation. For later reference, we note that in this case, n t λt  X Y P(Z 0 = z 0 ) = exp − λ wm (1) wn (sj (z 0 )) (6.3.13) t1 ! · · · tn ! m=1 j=1 Pn since 1 − wn (0) = m=1 wm (1) by Equations 6.3.4 and 6.3.5. 6.3.3.3 Left-ordered form Following Griffiths and Ghahramani (2005), we also consider the reduction to “left- ordered form”. Each column can be viewed as representing a number in binary, taking the first entry to be the most significant bit. For z ∈ {0, 1}n×k , define z l to be the matrix obtained by permuting the columns of z ∗ so that the associated binary numbers are nonincreasing. Such a matrix is said to be in left-ordered form. Similar to the reduction to staircase form, the distribution of Z l is t t! Y P(Z l = z l ) = Q2n −1 vn (t) wn (sj (z l )) (6.3.14) h=1 τh ! j=1 for z l in left-ordered form (with no empty columns), where τh is the number of P2n −1 columns of z l representing h in binary, and t = h=1 τh is the total number of 186 columns in z l . S∞ n×k Consider any two matrices z1 , z2 ∈ k=0 {0, 1} to be equivalent if they have the same left-ordered form, i.e., z1l = z2l . Note that for any permutation σ of the rows, we have z1l = z2l if and only if σ(z1 )l = σ(z2 )l ; consequently, any such σ induces a permutation on the set of equivalence classes. In the same way as the IBP, it is straightforward to show that the distribution on equivalence classes induced by the distribution on left-ordered matrices above (Equation 6.3.14) is exchangeable in the rows, in the sense that it is invariant under such permutations. 6.3.4 Buffet process When p(k) = Poisson(k|λ), the distribution of Z 0 can be described by a simple buffet process bearing a close resemblance to the Indian buffet process (IBP). Consider the following buffet process: For m = 1, . . . , n: Customer m . . . a + sj (1) tries each previously-tried dish j with probability , where a+b+m−1 sj is the number of previous customers trying dish j, and (2) tries a Poisson(λwm (1)) number of new dishes. (Customer 1 simply tries a Poisson(λw1 (1)) number of dishes, since there are no previously-tried dishes.) In the same way as the Indian buffet process, this gives rise to a binary matrix z 0 in the staircase form described above; specifically, labeling 0 dishes in the order they are tried, zmj = 1 if and only if customer m tries dish j. Note in particular that, in the notation of Section 6.3.3.2, tm is the number of new 187 Pn dishes tried by customer m, and t = m=1 tm is the total number of dishes tried. We claim that a random matrix generated by the buffet process above has the same distribution as Z 0 (Equation 6.3.11); see below for this argument. The only difference between the buffet process above and the IBP is that in the a + sj IBP, is replaced by sj /m, and λwm (1) is replaced by α/m. In a sense, a+b+m−1 the IBP can be viewed as a special case obtained by setting b = 1, λ = α/a, and taking the limit as a → 0. Further, the two-parameter generalization of the IBP described by Ghahramani et al. (2007) can be obtained similarly by setting b = β, λ = αβ/a, and taking a → 0. This two-parameter generalization was motivated by the desire to control the prior on the total number of dishes tried, separately from the number of dishes per customer. In addition to providing such control, the buffet process described above has a further degree of freedom through the parameter a. This has two effects: (1) the urn process for each column is controlled by both parameters a and b, instead of just b, and (2) the rate of decay of the mean number of new dishes per customer can be controlled by a, since ab(m−1) λaΓ(a + b) λ wm (1) = λ (m) ∼ (a + b) ma+1 Γ(b) as m → ∞ (by Stirling’s approximation for the gamma function); in contrast, in the two-parameter IBP the mean number of new dishes for customer m is αβ/(β +m−1). With this increased flexibility, a diverse variety of matrices can be generated by the MFFM. To illustrate, Figure 6.3 shows samples of Z 0 for various parameter settings. To control the total number of dishes tried, we reparametrize using n X  b(n)  c := λwm (1) = λ(1 − wn (0)) = λ 1 − m=1 (a + b)(n) 188 Figure 6.3: Sample matrices from the MFFM buffet process instead of λ, so that the total number of dishes tried is a Poisson(c) random variable. 6.3.4.1 Justification of the buffet process We show that the buffet process above generates Z 0 when p(k) = Poisson(k|λ). To see this, fix a z 0 obtained via the buffet process. Note that there is a single sequence of choices leading to z 0 , so the probability of getting z 0 is simply the product of the probabilities of these choices. Suppose mj is the first customer to try dish j, i.e., mj = mj (z 0 ). For each customer after mj , the choice of whether or not to try dish j follows the urn process of Section 6.3.2, and therefore, defining U1 , U2 , . . . as in 189 Section 6.3.2, these choices contribute the factor 0 0  P Un = znj , . . . , Umj +1 = zm j +1,j | Umj = 1, Ui = 0 ∀i < mj 0 P(U1:n = z1:n,j ) wn (sj (z 0 )) = = (6.3.15) P(Umj = 1, Ui = 0 ∀i < mj ) wmj (1) to the probability of z 0 , where the last step follows by applying Equation 6.3.2 to both numerator and denominator. Multiplying these for j = 1, . . . , t, where t is the number of columns in z 0 , accounts for the choices in step (1) of the buffet process. For the choices in step (2), let tm denote the number of new dishes tried by customer m, and note that t = nm=1 tm . The probability of choosing t1 , . . . , tn is P n n Y Y (λwm (1))tm Poisson(tm |λwm (1)) = exp(−λwm (1)) m=1 m=1 tm !  Xn  λt n Y = exp − λ wm (1) wm (1)tm . m=1 t1 ! · · · t ! n m=1 Multiplying this by Equation 6.3.15 for each j = 1, . . . , t, the wm (1) factors cancel and we obtain Equation 6.3.13, as desired. 6.3.4.2 Exchangeable buffet process As with the IBP, although customers are not exchangeable under the buffet process above, this can be achieved by a certain modification (in the same way as the IBP). This involves generating a matrix in left-ordered form with the distribution of Z l , and considering the associated equivalence class. Suppose that the matrix generated by the customers before customer m has τh columns representing h in binary, for h =  a + sh  1, . . . , 2m−1 − 1. For each h, customer m draws L ∼ Binomial τh , a+b+m−1 and tries the first L of the τh associated dishes, where sh is the number of ones in 190 the binary representation of h. Then, as before, he tries a Poisson(λwm (1)) number of new dishes. It can be shown that the probability of obtaining z l by this process is P(Z l = z l ) (as in Equations 6.3.14 and 6.3.12) when p(k) = Poisson(k|λ). In the distribution on equivalence classes as in Section 6.3.3.3, the rows are exchangeable, and therefore the customers in this process are exchangeable in the sense that following the process for any ordering of the customers will yield the same distribution on equivalence classes. 6.3.5 Inference The Gibbs sampling algorithm for the IBP described by Griffiths and Ghahramani (2005) is easily adapted to the MFFM. However, it is well-known that mixing can be very slow when using Gibbs sampling with the IBP, and the same issue is present with the MFFM. Other inference algorithms have also been proposed for the IBP (e.g., Teh et al. (2007), Doshi-Velez et al. (2008)), and it might be interesting to see if these approaches can also be adapted to MFFM. Chapter Seven Conclusion 192 In this thesis, we have seen that variable-dimension models can be an appealing alternative to the commonly-used infinite-dimensional nonparametric models. In closing, we will speculatively discuss some open questions and possibilities for future work. The inconsistency results of Chapters 2 and 3 — and more fundamentally, the properties of the partition distribution that lead to this inconsistency — clearly in- dicate that when the data more closely resembles a finite mixture than an infinite mixture from the assumed family, naively estimating heterogeneity using the number of clusters in an infinite mixture can be problematic. In these situations, estimat- ing heterogeneity using a variable-dimension model may be preferable to using an infinite-dimensional model, however, the effect of misspecification of the component distributions needs to be carefully considered; finding a principled way to handle misspecification would be interesting. Although we have shown that the DPM posterior on the number of clusters does not concentrate, we have not been able to determine the precise limiting behavior. As a theoretical question, this is challenging and might be interesting to know, although the practical relevance may be limited. The empirical studies of Chapter 5 indicate that for density estimation, the MFM and the DPM seem essentially indistinguishable. On the other hand, they can be quite different for clustering, and perhaps the most significant difference is in the posteriors on the number of components and clusters, although interestingly, using a DPM with random α seems to be more similar to the MFM than a DPM with fixed α. It would be interesting to see practical examples of applications where these differences matter. 193 Regarding inference algorithms, there are several questions that would be in- teresting to explore. The slow mixing time of the incremental samplers when n is large is troublesome, but it seems likely that using one of the existing split-merge algorithms would resolve this issue. Using a smaller value of γ may also help in this regard. With the non-conjugate sampler, are there practical examples in which either DPM or MFM mixing can be significantly improved by approximating the single-cluster posterior, rather than sampling from the prior, for the auxiliary vari- able distribution? Would it make sense to use stick-breaking representations for inference in variable-dimensional models? It is striking that so many of the elegant properties of the three nonparametric models considered here — the Dirichlet process, the hierarchical Dirichlet process, and the Indian buffet process — are also exhibited by their variable-dimension coun- terparts (Chapters 4 and 6). It would be interesting to see how far this can be taken. Do all nonparametric models have a variable-dimension counterpart that exhibits similar properties? Is there a general principle at play here? Is there a broader perspective from which all of these models can be viewed? Bibliography M. Aitkin. Likelihood and Bayesian analysis of mixtures. Statistical Modelling, 1(4): 287–304, 2001. C. E. Antoniak. Mixtures of Dirichlet processes with applications to Bayesian non- parametric problems. The Annals of Statistics, 2(6):1152–1174, November 1974. A. Azzalini and A. Capitanio. Statistical applications of the multivariate skew- normal distribution. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 61(3):579–602, 1999. A. Azzalini and A. Dalla Valle. The multivariate skew-normal distribution. Biometrika, 83(4):715–726, 1996. A. Barron, M. J. Schervish, and L. Wasserman. The consistency of posterior dis- tributions in nonparametric problems. The Annals of Statistics, 27(2):536–561, 1999. D. Barry and J. A. Hartigan. Product partition models for change point problems. The Annals of Statistics, pages 260–279, 1992. D. A. Berry and R. Christensen. Empirical Bayes estimation of a binomial parameter via mixtures of Dirichlet processes. The Annals of Statistics, pages 558–568, 1979. J. Besag and P. J. Green. Spatial statistics and Bayesian computation. Journal of the Royal Statistical Society. Series B (Methodological), pages 25–37, 1993. A. Bhattacharya and D. B. Dunson. Nonparametric Bayesian density estimation on manifolds with applications to planar shapes. Biometrika, 97(4):851–865, 2010. D. Blackwell and J. B. MacQueen. Ferguson distributions via P´olya urn schemes. The Annals of Statistics, pages 353–355, 1973. D. M. Blei and M. I. Jordan. Variational inference for Dirichlet process mixtures. Bayesian Analysis, 1(1):121–143, 2006. L. Breiman. Probability (Classics in Applied Mathematics). Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1992. T. Broderick, M. I. Jordan, and J. Pitman. Beta processes, stick-breaking and power laws. Bayesian Analysis, 7(2):439–476, 2012. 194 195 C. A. Bush and S. N. MacEachern. A semiparametric Bayesian model for randomised block designs. Biometrika, 83(2):275–285, 1996. C. M. Carvalho, J. Chang, J. E. Lucas, J. R. Nevins, Q. Wang, and M. West. High- dimensional sparse factor modeling: applications in gene expression genomics. Journal of the American Statistical Association, 103(484), 2008. A. Cerquetti. Generalized Chinese restaurant construction of exchangeable Gibbs partitions and related results. arXiv:0805.3853, 2008. A. Cerquetti. Conditional α-diversity for exchangeable Gibbs partitions driven by the stable subordinator. arXiv:1105.0892, 2011. H. Chen, P. L. Morrell, V. E. T. M. Ashworth, M. de la Cruz, and M. T. Clegg. Tracing the geographic origins of major avocado cultivars. Journal of Heredity, 100(1):56–65, 2009. J. Chen, P. Li, and Y. Fu. Inference on the order of a normal mixture. Journal of the American Statistical Association, 107(499):1096–1105, 2012. W.-C. Chen. On the weak form of Zipf’s law. Journal of Applied Probability, pages 611–622, 1980. Y. Chung and D. B. Dunson. Nonparametric Bayes conditional distribution modeling with variable selection. Journal of the American Statistical Association, 104(488), 2009. D. B. Dahl. An improved merge-split sampler for conjugate Dirichlet process mix- ture models. Technical Report, Department of Statistics, University of Wisconsin Madison, 2003. D. B. Dahl. Sequentially-allocated merge-split sampler for conjugate and nonconju- gate Dirichlet process mixture models. Journal of Computational and Graphical Statistics, 11, 2005. D. B. Dahl. Model-based clustering for expression data via a Dirichlet process mix- ture model. Bayesian inference for gene expression and proteomics, pages 201–218, 2006. D. B. Dahl. Modal clustering in a class of product partition models. Bayesian Analysis, 4(2):243–264, 2009. N. G. De Bruijn. Asymptotic Methods in Analysis. North-Holland Publishing Co., Amsterdam, 1970. P. Diaconis and D. Freedman. On the consistency of Bayes estimates. The Annals of Statistics, 14(1):1–26, 1986. P. Diaconis and D. Ylvisaker. Conjugate priors for exponential families. The Annals of Statistics, 7(2):269–281, 1979. J. L. Doob. Application of the theory of martingales. In Actes du Colloque Inter- national Le Calcul des Probabilit´es et ses applications (Lyon, 28 Juin – 3 Juillet, 1948), pages 23–27. Paris CNRS, 1949. 196 F. Doshi-Velez, K. T. Miller, J. Van Gael, and Y. W. Teh. Variational inference for the Indian buffet process. In AISTATS 2009, pages 137–144, 2008. M. J. Drinkwater, Q. A. Parker, D. Proust, E. Slezak, and H. Quintana. The large scale distribution of galaxies in the Shapley supercluster. Publications of the As- tronomical Society of Australia, 21(1):89–96, 2004. D. B. Dunson. Efficient Bayesian model averaging in factor analysis. Technical Report, ISDS, Duke University, 2006. D. B. Dunson and J.-H. Park. Kernel stick-breaking processes. Biometrika, 95(2): 307–323, 2008. R. Durrett. Probability: Theory and Examples, volume 2. Cambridge University Press, 1996. R. G. Edwards and A. D. Sokal. Generalization of the Fortuin-Kasteleyn-Swendsen- Wang representation and Monte Carlo algorithm. Physical Review D, 38:2009– 2012, 1988. M. D. Escobar. Estimating the means of several normal populations by nonparametric estimation of the distribution of the means. PhD thesis, Yale University, 1988. M. D. Escobar and M. West. Bayesian density estimation and inference using mix- tures. Journal of the American Statistical Association, 90(430):577–588, 1995. M. D. Escobar and M. West. Computing nonparametric hierarchical models. In D. Dey, P. M¨ uller, and D. Sinha, editors, Practical Nonparametric and Semipara- metric Bayesian Statistics, pages 1–22. Springer-Verlag, New York, 1998. S. Favaro, A. Lijoi, and I. Pruenster. On the stick-breaking representation of nor- malized inverse Gaussian priors. Biometrika, 99(3):663–674, 2012. P. Fearnhead. Particle filters for mixture models with an unknown number of com- ponents. Statistics and Computing, 14(1):11–21, 2004. T. S. Ferguson. A Bayesian analysis of some nonparametric problems. The Annals of Statistics, pages 209–230, 1973. T. S. Ferguson. Bayesian density estimation by mixtures of normal distributions. In M. H. Rizvi, J. Rustagi, and D. Siegmund, editors, Recent Advances in Statistics, pages 287–302. Academic Press, 1983. E. B. Fox, E. B. Sudderth, and A. S. Willsky. Hierarchical Dirichlet processes for tracking maneuvering targets. In 10th International Conference on Information Fusion, 2007, pages 1–8. IEEE, 2007. E. B. Fox, E. B. Sudderth, M. I. Jordan, and A. S. Willsky. Sharing features among dynamical systems with beta processes. In Advances in Neural Information Pro- cessing Systems, pages 549–557, 2009. D. A. Freedman. On the asymptotic behavior of Bayes’ estimates in the discrete case. The Annals of Mathematical Statistics, pages 1386–1403, 1963. 197 Z. Ghahramani, T. L. Griffiths, and P. Sollich. Bayesian nonparametric latent feature models. 2007. S. Ghosal. The Dirichlet process, related priors and posterior asymptotics. In N. L. Hjort, C. Holmes, P. M¨ uller, and S. G. Walker, editors, Bayesian Nonparametrics, pages 36–83. Cambridge University Press, 2010. S. Ghosal and A. Van der Vaart. Posterior convergence rates of Dirichlet mixtures at smooth densities. The Annals of Statistics, 35(2):697–723, 2007. S. Ghosal and A. W. Van der Vaart. Entropies and rates of convergence for maximum likelihood and Bayes estimation for mixtures of normal densities. The Annals of Statistics, pages 1233–1263, 2001. S. Ghosal, J. K. Ghosh, and R. V. Ramamoorthi. Posterior consistency of Dirichlet mixtures in density estimation. The Annals of Statistics, 27(1):143–158, 1999. J. Ghosh and D. B. Dunson. Default prior distributions and efficient posterior com- putation in Bayesian factor analysis. Journal of Computational and Graphical Statistics, 18(2):306–320, 2009. J. K. Ghosh and R. V. Ramamoorthi. Bayesian Nonparametrics. Springer–Verlag, New York, 2003. S. K. Ghosh and S. Ghosal. Semiparametric accelerated failure time models for censored data. Bayesian Statistics and its Applications, pages 213–229, 2006. A. Gnedin. A species sampling model with finitely many types. Elect. Comm. Probab., 15:79–88, 2010. A. Gnedin and J. Pitman. Exchangeable Gibbs partitions and Stirling triangles. Journal of Mathematical Sciences, 138(3):5674–5685, 2006. E. G. Gonzalez and R. Zardoya. Relative role of life-history traits and historical factors in shaping genetic population structure of sardines (Sardina pilchardus). BMC evolutionary biology, 7(1):197, 2007. R. L. Graham, D. E. Knuth, and O. Patashnik. Concrete Mathematics. Addison- Wesley, 1989. P. J. Green and S. Richardson. Modeling heterogeneity with and without the Dirich- let process. Scandinavian Journal of Statistics, 28(2):355–375, June 2001. J. E. Griffin and M. J. Steel. Order-based dependent Dirichlet processes. Journal of the American Statistical Association, 101(473):179–194, 2006. T. L. Griffiths and Z. Ghahramani. Infinite latent feature models and the Indian buffet process. In Advances in Neural Information Processing Systems, volume 18, pages 475–482, 2005. T. L. Griffiths, K. R. Canini, A. N. Sanborn, and D. J. Navarro. Unifying rational models of categorization via the hierarchical Dirichlet process. In Proceedings of the 29th Annual Cognitive Science Society, pages 323–328, 2007. 198 B. Hansen and J. Pitman. Prediction rules for exchangeable sequences related to species sampling. Statistics & probability letters, 46(3):251–256, 2000. J. Henna. On estimating of the number of constituents of a finite mixture of con- tinuous distributions. Annals of the Institute of Statistical Mathematics, 37(1): 235–240, 1985. J. Henna. Estimation of the number of components of finite mixtures of multivariate distributions. Annals of the Institute of Statistical Mathematics, 57(4):655–664, 2005. N. L. Hjort. Bayesian analysis for a generalised Dirichlet process prior. Technical Report, University of Oslo, 2000. M.-W. Ho, L. F. James, and J. W. Lau. Gibbs partitions (EPPF’s) derived from a stable subordinator are Fox H and Meijer G transforms. arXiv:0708.0619, 2007. W. Hoeffding. The strong law of large numbers for U-statistics. Institute of Statistics, Univ. of N. Carolina, Mimeograph Series, 302, 1961. M. D. Hoffman, D. M. Blei, and P. R. Cook. Content-based musical similarity computation using the hierarchical Dirichlet process. In ISMIR, pages 349–354, 2008. J. Hoffmann-Jørgensen. Probability with a View toward Statistics, volume 2. Chap- man & Hall, 1994. J. P. Huelsenbeck and P. Andolfatto. Inference of population structure under a Dirichlet process model. Genetics, 175(4):1787–1802, 2007. H. Ishwaran and L. F. James. Gibbs sampling methods for stick-breaking priors. Journal of the American Statistical Association, 96(453), 2001. H. Ishwaran and L. F. James. Generalized weighted Chinese restaurant processes for species sampling mixture models. Statistica Sinica, 13(4):1211–1236, 2003. H. Ishwaran and M. Zarepour. Markov chain Monte Carlo in approximate Dirichlet and beta two-parameter process hierarchical models. Biometrika, 87(2):371–390, 2000. H. Ishwaran, L. F. James, and J. Sun. Bayesian model selection in finite mixtures by marginal density decompositions. Journal of the American Statistical Association, 96(456), 2001. S. Jain and R. M. Neal. A split-merge Markov chain Monte Carlo procedure for the Dirichlet process mixture model. Journal of Computational and Graphical Statistics, 13(1), 2004. S. Jain and R. M. Neal. Splitting and merging components of a nonconjugate Dirich- let process mixture model. Bayesian Analysis, 2(3):445–472, 2007. L. F. James. Large sample asymptotics for the two-parameter Poisson–Dirichlet process. In Pushing the Limits of Contemporary Statistics: Contributions in Honor of Jayanta K. Ghosh, pages 187–199. Institute of Mathematical Statistics, 2008. 199 L. F. James, C. E. Priebe, and D. J. Marchette. Consistent estimation of mixture complexity. The Annals of Statistics, pages 1281–1296, 2001. G. H. Jang, J. Lee, and S. Lee. Posterior consistency of species sampling priors. Statistica Sinica, 20(2):581, 2010. A. Jasra, C. Holmes, and D. Stephens. Markov chain Monte Carlo methods and the label switching problem in Bayesian mixture modeling. Statistical Science, pages 50–67, 2005. M. Kalli, J. E. Griffin, and S. G. Walker. Slice sampling mixture models. Statistics and computing, 21(1):93–105, 2011. C. Keribin. Consistent estimation of the order of mixture models. Sankhya Ser. A, 62(1):49–66, 2000. S. Khazaei, J. Rousseau, and F. Balabdaoui. Nonparametric Bayesian estimation of densities under monotonicity constraint. (preprint), 2012. S. Kim, M. G. Tadesse, and M. Vannucci. Variable selection in clustering via Dirichlet process mixture models. Biometrika, 93(4):877–893, 2006. J. J. Kivinen, E. B. Sudderth, and M. I. Jordan. Learning multiscale representations of natural scenes using Dirichlet processes. In Computer Vision, 2007. ICCV 2007. IEEE 11th International Conference on, pages 1–8. IEEE, 2007. A. W. Knapp. Basic Real Analysis. Birkh¨auser, 2005. D. Knowles and Z. Ghahramani. Infinite sparse factor analysis and infinite inde- pendent components analysis. In Independent Component Analysis and Signal Separation, pages 381–388. Springer, 2007. S. G. Krantz. Function Theory of Several Complex Variables. AMS Chelsea Pub- lishing, Providence, 1992. R. Krause, D. L. Wild, W. Chu, and Z. Ghahramani. Identifying protein complexes in high-throughput protein interaction screens using an infinite latent feature model. In Pacific Symposium on Biocomputing, volume 11, pages 231–242, 2006. W. Kruijer. Convergence rates in nonparametric Bayesian density estimation. PhD thesis, Department of Mathematics, Vrije Universiteit Amsterdam, 2008. W. Kruijer, J. Rousseau, and A. Van der Vaart. Adaptive Bayesian density estima- tion with location-scale mixtures. Electronic Journal of Statistics, 4:1225–1257, 2010. N. Lartillot and H. Philippe. A Bayesian mixture model for across-site heterogeneities in the amino-acid replacement process. Molecular Biology and Evolution, 21(6): 1095–1109, 2004. A. D. Leach´e and M. K. Fujita. Bayesian species delimitation in West African forest geckos (Hemidactylus fasciatus). Proceedings of the Royal Society B: Biological Sciences, 277(1697):3071–3077, 2010. 200 J. K. Lee, V. Danˇc´ık, and M. S. Waterman. Estimation for restriction sites observed by optical mapping using reversible-jump Markov chain Monte Carlo. Journal of Computational Biology, 5(3):505–515, 1998. B. G. Leroux. Consistent estimation of a mixing distribution. The Annals of Statis- tics, 20(3):1350–1360, 1992. P. Liang, S. Petrov, M. I. Jordan, and D. Klein. The infinite PCFG using hierarchical Dirichlet processes. In EMNLP-CoNLL, pages 688–697, 2007. A. Lijoi and I. Pr¨ unster. Models beyond the Dirichlet process. Bayesian Nonpara- metrics, 28:80, 2010. A. Lijoi, R. H. Mena, and I. Pr¨ unster. Hierarchical mixture modeling with normalized inverse-Gaussian priors. Journal of the American Statistical Association, 100(472): 1278–1291, 2005a. A. Lijoi, I. Pr¨ unster, and S. G. Walker. On consistency of nonparametric normal mixtures for Bayesian density estimation. Journal of the American Statistical Association, 100(472):1292–1296, 2005b. A. Lijoi, R. H. Mena, and I. Pr¨unster. Bayesian nonparametric estimation of the probability of discovering new species. Biometrika, 94(4):769–786, 2007. A. Lijoi, I. Pr¨ unster, and S. G. Walker. Bayesian nonparametric estimators derived from conditional Gibbs structures. The Annals of Applied Probability, 18(4):1519– 1547, 2008. J. S. Liu. The collapsed Gibbs sampler in Bayesian computations with applications to a gene regulation problem. Journal of the American Statistical Association, 89 (427):958–966, 1994. A. Y. Lo. On a class of Bayesian nonparametric estimates: I. Density estimates. The Annals of Statistics, 12(1):351–357, 1984. H. F. Lopes and M. West. Bayesian model assessment in factor analysis. Statistica Sinica, 14(1):41–68, 2004. E. D. Lorenzen, P. Arctander, and H. R. Siegismund. Regional genetic structuring and evolutionary history of the impala Aepyceros melampus. Journal of Heredity, 97(2):119–132, 2006. S. N. MacEachern. Estimating normal means with a conjugate style Dirichlet process prior. Communications in Statistics-Simulation and Computation, 23(3):727–741, 1994. S. N. MacEachern. Computational methods for mixture of Dirichlet process models. In Practical nonparametric and semiparametric Bayesian statistics, pages 23–43. Springer, 1998. S. N. MacEachern. Dependent dirichlet processes. Unpublished manuscript, Depart- ment of Statistics, The Ohio State University, 2000. 201 S. N. MacEachern and P. M¨ uller. Estimating mixture of Dirichlet process models. Journal of Computational and Graphical Statistics, 7(2):223–238, 1998. S. N. MacEachern, M. Clyde, and J. S. Liu. Sequential importance sampling for nonparametric Bayes models: The next generation. Canadian Journal of Statistics, 27(2):251–267, 1999. S. N. MacEachern, A. Kottas, and A. Gelfand. Spatial nonparametric Bayesian models. In Proceedings of the 2001 Joint Statistical Meetings, volume 3, 2001. J. Marin, M. Rodriguez-Bernal, and M. Wiper. Using Weibull mixture distributions to model heterogeneous survival data. Communications in StatisticsSimulation and Computation, 34(3):673–684, 2005. P. McCullagh and J. Yang. How many clusters? Bayesian Analysis, 3(1):101–120, 2008. M. Medvedovic and S. Sivaganesan. Bayesian infinite mixture model based clustering of gene expression profiles. Bioinformatics, 18(9):1194–1206, 2002. E. Meeds, Z. Ghahramani, R. M. Neal, and S. T. Roweis. Modeling dyadic data with binary latent factors. Advances in Neural Information Processing Systems, 19:977, 2007. S. P. Meyn and R. L. Tweedie. Markov Chains and Stochastic Stability. Cambridge University Press, 2009. J. W. Miller and M. T. Harrison. A simple example of Dirichlet process mixture inconsistency for the number of components. In Advances in Neural Information Processing Systems, Vol. 26, 2013a. J. W. Miller and M. T. Harrison. Inconsistency of Pitman–Yor process mixtures for the number of components. arXiv:1309.0024, 2013b. K. T. Miller, T. L. Griffiths, and M. I. Jordan. Nonparametric latent feature mod- els for link prediction. In Advances in Neural Information Processing Systems, volume 9, pages 1276–1284, 2009. P. M¨ uller and F. Quintana. Random partition models with regression on covariates. Journal of Statistical Planning and Inference, 140(10):2801–2808, 2010. P. M¨uller, F. Quintana, and G. Rosner. A method for combining inference across related nonparametric Bayesian models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 66(3):735–749, 2004. P. M¨ uller, F. Quintana, and G. L. Rosner. A product partition model with regression on covariates. Journal of Computational and Graphical Statistics, 20(1), 2011. L. Murino, C. Angelini, I. Bifulco, I. De Feis, G. Raiconi, and R. Tagliaferri. Mul- tiple clustering solutions analysis through least-squares consensus algorithms. In Computational Intelligence Methods for Bioinformatics and Biostatistics, pages 215–227. Springer, 2010. 202 D. J. Navarro and T. L. Griffiths. A nonparametric Bayesian method for inferring features from similarity judgments. Advances in Neural Information Processing Systems, 19:1033, 2007. R. M. Neal. Bayesian mixture modeling. In Maximum Entropy and Bayesian Meth- ods, pages 197–211. Springer, 1992. R. M. Neal. Markov chain sampling methods for Dirichlet process mixture models. Journal of Computational and Graphical Statistics, 9(2):249–265, 2000. D. P. Nguyen, L. M. Frank, and E. N. Brown. An application of reversible-jump Markov chain Monte Carlo to spike classification of multi-unit extracellular record- ings. Network (Bristol, England), 14(1):61–82, 2003. X. L. Nguyen. Convergence of latent mixing measures in finite and infinite mixture models. The Annals of Statistics, 41(1):370–400, 2013. A. Nobile. Bayesian Analysis of Finite Mixture Distributions. PhD thesis, Depart- ment of Statistics, Carnegie Mellon University, Pittsburgh, PA, 1994. A. Nobile. Bayesian finite mixtures: a note on prior specification and posterior computation. Technical Report, Department of Statistics, University of Glasgow, 2005. A. Nobile and A. T. Fearnside. Bayesian finite mixtures with an unknown number of components: The allocation sampler. Statistics and Computing, 17(2):147–162, 2007. A. Onogi, M. Nurimoto, and M. Morita. Characterization of a Bayesian genetic clustering algorithm based on a Dirichlet process prior and comparison among Bayesian clustering methods. BMC Bioinformatics, 12(1):263, 2011. E. Otranto and G. M. Gallo. A nonparametric Bayesian approach to detect the number of regimes in Markov switching models. Econometric Reviews, 21(4):477– 496, 2002. J. W. Paisley, A. K. Zaas, C. W. Woods, G. S. Ginsburg, and L. Carin. A stick- breaking construction of the beta process. In Proceedings of the 27th International Conference on Machine Learning, pages 847–854, 2010. O. Papaspiliopoulos and G. O. Roberts. Retrospective Markov chain Monte Carlo methods for Dirichlet process hierarchical models. Biometrika, 95(1):169–186, 2008. J.-H. Park and D. B. Dunson. Bayesian generalized product partition model. Sta- tistica Sinica, 20:1203–1226, 2010. D. Pati, D. B. Dunson, and S. T. Tokdar. Posterior consistency in conditional distribution estimation. Journal of Multivariate Analysis, 116:456–472, 2013. J. Pella and M. Masuda. The Gibbs and split-merge sampler for population mix- ture analysis from genetic data with incomplete baselines. Canadian Journal of Fisheries and Aquatic Sciences, 63(3):576–596, 2006. 203 M. Perman, J. Pitman, and M. Yor. Size-biased sampling of Poisson point processes and excursions. Probability Theory and Related Fields, 92(1):21–39, 1992. D. B. Phillips and A. F. M. Smith. Bayesian model comparison via jump diffusions. In Markov chain Monte Carlo in Practice, pages 215–239. Springer, 1996. J. Pitman. Exchangeable and partially exchangeable random partitions. Probability Theory and Related Fields, 102(2):145–158, 1995. J. Pitman. Some developments of the Blackwell-MacQueen urn scheme. Lecture Notes-Monograph Series, pages 245–267, 1996. J. Pitman. Combinatorial Stochastic Processes. Springer–Verlag, Berlin, 2006. J. Pitman and M. Yor. The two-parameter Poisson–Dirichlet distribution derived from a stable subordinator. The Annals of Probability, 25(2):855–900, 1997. J. K. Pritchard, M. Stephens, and P. Donnelly. Inference of population structure using multilocus genotype data. Genetics, 155(2):945–959, 2000. F. A. Quintana. A predictive view of Bayesian clustering. Journal of Statistical Planning and Inference, 136(8):2407–2429, 2006. F. A. Quintana and P. L. Iglesias. Bayesian clustering and product partition models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 65(2): 557–574, 2003. L. Ren, D. B. Dunson, and L. Carin. The dynamic hierarchical Dirichlet process. In Proceedings of the 25th International Conference on Machine Learning, pages 824–831. ACM, 2008. C. M. Richards, G. M. Volk, A. A. Reilley, A. D. Henk, D. R. Lockwood, P. A. Reeves, and P. L. Forsline. Genetic diversity and population structure in Malus sieversii, a wild progenitor species of domesticated apple. Tree Genetics & Genomes, 5(2): 339–347, 2009. S. Richardson and P. J. Green. On Bayesian analysis of mixtures with an unknown number of components. Journal of the Royal Statistical Society, Series B, 59(4): 731–792, 1997. C. P. Robert and G. Casella. Monte Carlo Statistical Methods, volume 319. Springer, 2004. R. T. Rockafellar. Convex Analysis. Princeton University Press, 1970. A. Rodriguez and D. B. Dunson. Nonparametric Bayesian models through probit stick-breaking processes. Bayesian Analysis, 6(1), 2011. A. Rodriguez, D. B. Dunson, and A. E. Gelfand. The nested Dirichlet process. Journal of the American Statistical Association, 103(483), 2008. K. Roeder. Density estimation with confidence sets exemplified by superclusters and voids in the galaxies. Journal of the American Statistical Association, 85(411): 617–624, 1990. 204 K. Roeder. A graphical technique for determining the number of components in a mixture of normals. Journal of the American Statistical Association, 89(426): 487–495, 1994. J. Rousseau and K. Mengersen. Asymptotic behaviour of the posterior distribution in overfitted mixture models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(5):689–710, 2011. S. Roweis and Z. Ghahramani. A unifying review of linear Gaussian models. Neural Computation, 11(2):305–345, 1999. T. Sapatinas. Identifiability of mixtures of power-series distributions and related characterizations. Annals of the Institute of Statistical Mathematics, 47(3):447– 459, 1995. C. Scricciolo. Adaptive Bayesian density estimation using Pitman–Yor or normalized inverse-Gaussian process kernel mixtures. arXiv:1210.8094, 2012. J. Sethuraman. A constructive definition of Dirichlet priors. Statistica Sinica, 4: 639–650, 1994. J. Sethuraman and R. C. Tiwari. Convergence of Dirichlet measures and the inter- pretation of their parameter. Technical Report, Department of Statistics, Florida State University, 1981. M. Stephens. Bayesian analysis of mixture models with an unknown number of components—An alternative to reversible jump methods. The Annals of Statistics, 28(1):40–74, 2000. Y. Tang and S. Ghosal. Posterior consistency of Dirichlet mixtures for estimating a transition density. Journal of Statistical Planning and Inference, 137(6):1711– 1726, June 2007. Y. W. Teh, M. I. Jordan, M. J. Beal, and D. M. Blei. Sharing clusters among related groups: hierarchical Dirichlet processes. In Advances in Neural Information Processing Systems, pages 1385–1392, 2004. Y. W. Teh, M. I. Jordan, M. J. Beal, and D. M. Blei. Hierarchical Dirichlet processes. Journal of the American Statistical Association, 101(476):1566–1581, 2006. Y. W. Teh, D. G¨or¨ur, and Z. Ghahramani. Stick-breaking construction for the Indian buffet process. In International Conference on Artificial Intelligence and Statistics, pages 556–563, 2007. H. Teicher. Identifiability of finite mixtures. The Annals of Mathematical Statistics, 34(4):1265–1269, 1963. R. Thibaux and M. I. Jordan. Hierarchical beta processes and the Indian buffet process. In International Conference on Artificial Intelligence and Statistics, pages 564–571, 2007. L. Tierney. Markov chains for exploring posterior distributions. The Annals of Statistics, pages 1701–1728, 1994. 205 S. T. Tokdar. Posterior consistency of Dirichlet location-scale mixture of normals in density estimation and regression. Sankhy¯a: The Indian Journal of Statistics, pages 90–110, 2006. G. Tomlinson and M. Escobar. Analysis of densities. Technical Report, University of Toronto, 1999. C. Vogl, F. Sanchez-Cabo, G. Stocker, S. Hubbard, O. Wolkenhauer, and Z. Tra- janoski. A fully Bayesian model to cluster gene-expression profiles. Bioinformatics, 21(suppl 2):ii130–ii136, 2005. S. G. Walker, A. Lijoi, and I. Pr¨ unster. On rates of convergence for posterior distri- butions in infinite-dimensional models. The Annals of Statistics, 35(2):738–746, 2007. M. West. Hyperparameter estimation in Dirichlet process mixture models. ISDS Discussion Paper #92-A03, Duke University, 1992. M. West. Bayesian factor regression models in the “large p, small n” paradigm. In J. M. Bernardo, M. J. Bayarri, J. O. Berger, A. P. Dawid, D. Heckerman, A. F. M. Smith, and M. West, editors, Bayesian Statistics 7. Oxford University Press, 2003. M. West, P. M¨uller, and M. D. Escobar. Hierarchical priors and mixture models, with application in regression and density estimation. In P. Freeman and A. F. Smith, editors, Aspects of Uncertainty: A Tribute to D.V. Lindley, pages 363–386. Wiley, 1994. M.-J. Woo and T. Sriram. Robust estimation of mixture complexity for count data. Computational Statistics and Data Analysis, 51(9):4379–4392, 2007. M.-J. Woo and T. N. Sriram. Robust estimation of mixture complexity. Journal of the American Statistical Association, 101(476), 2006. F. Wood, T. L. Griffiths, and Z. Ghahramani. A non-parametric Bayesian method for inferring hidden causes. In Proceedings of the 22nd Conference on Uncertainty in Artificial Intelligence (UAI), 2006. Y. Wu and S. Ghosal. The L1 -consistency of Dirichlet mixtures in multivariate Bayesian density estimation. Journal of Multivariate Analysis, 101(10):2411–2419, November 2010. E. P. Xing, K. A. Sohn, M. I. Jordan, and Y. W. Teh. Bayesian multi-population haplotype inference via a hierarchical Dirichlet process mixture. In Proceedings of the 23rd International Conference on Machine Learning, pages 1049–1056, 2006. S. J. Yakowitz and J. D. Spragins. On the identifiability of finite mixtures. The Annals of Mathematical Statistics, 39(1):209–214, 1968.