Science, Faculty of
Statistics, Department of
DSpace
UBCV
Li, Xinglong
Doctor of Philosophy - PhD
University of British Columbia
Bayesian nonparametric (BNP) models provide a flexible and powerful framework for statistical modeling by allowing the number of features or subgroups within a population to grow with the data volume. However, posterior inference in BNP models is challenging due to the infinite parameters involved, and the lack of general and efficient inference procedures impedes their practical application. Exact posterior inference methods either analytically marginalize out the infinitely many parameters, or introduce auxiliary variables to adaptively adjust the model size during inference. The former approach relies on conjugacy relationships between priors and likelihoods and suffers from high computational costs. Similarly, the latter approach is also computationally demanding, as it requires numerical integration during sampling for nonconjugate models. An alternative common practice in fitting BNP models involves approximating the nonparametric model with a parametric one, and subsequently applying a standard inference algorithm. While this is practical, parametric truncation can lead to significant unknown posterior approximation errors, particularly for BNP models with heavy tails that support the power-law behavior of the population. Previous work on truncated inference in BNP models has determined the truncation level via analysis of the forward generative model, which does not accurately reflect the error of approximation of the target posterior distribution. This thesis aims to develop approximate inference algorithms that can be directly used for posterior inference for general BNP models. We propose truncated inference methods and provide estimates of the posterior truncation error. Rather than setting the truncation level based on prior approximation error, we establish a desired posterior truncation error level, allowing the algorithm to adapt the truncation level until the desired truncation error is reached. The proposed algorithms are general in that they can be applied to a wide range of BNP models with completely random measure (CRM) priors. We have applied these algorithms to edge-exchangeable network models, where feature assignment variables are observed, and to latent feature models with latent feature assignment variables.
Approximate posterior inference forBayesian nonparametrics with totalvariation error boundsbyXinglong LiB.Sc., China University of Petroleum (East China), 2015M.Sc., Yunnan University, 2018A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinTHE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES(Statistics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)April 2025© Xinglong Li, 2025The following individuals certify that they have read, and recommend tothe Faculty of Graduate and Postdoctoral Studies for acceptance, the thesisentitled:Approximate posterior inference for Bayesian nonparametrics withtotal variation error boundssubmitted by Xinglong Li in partial fulfillment of the requirements for thedegree of Doctor of Philosophy in Statistics.Examining Committee:Trevor Campbell, Associate Professor, Department of Statistics, Universityof British ColumbiaSupervisorAlexandre Bouchard-Coˆte´, Professor, Department of Statistics, University ofBritish ColumbiaSupervisory Committee MemberBen Bloem-Reddy, Assistant Professor, Department of Statistics, Universityof British ColumbiaSupervisory Committee MemberPaul Gustafson, Professor, Department of Statistics, University of BritishColumbiaUniversity ExaminerFrank Wood, Professor, Department of Computer Science, University ofBritish ColumbiaUniversity ExaminerSinead Williamson, Associate Professor, Department of Statistics and DataSciences , The University of Texas at AustinExternal ExamineriiAbstractBayesian nonparametric (BNP) models provide a flexible and powerful frame-work for statistical modeling by allowing the number of features or subgroupswithin a population to grow with the data volume. However, posteriorinference in BNP models is challenging due to the infinite parameters in-volved, and the lack of general and efficient inference procedures impedestheir practical application. Exact posterior inference methods either analyti-cally marginalize out the infinitely many parameters, or introduce auxiliaryvariables to adaptively adjust the model size during inference. The formerapproach relies on conjugacy relationships between priors and likelihoodsand suffers from high computational costs. Similarly, the latter approach isalso computationally demanding, as it requires numerical integration duringsampling for nonconjugate models.An alternative common practice in fitting BNP models involves approxi-mating the nonparametric model with a parametric one, and subsequentlyapplying a standard inference algorithm. While this is practical, parametrictruncation can lead to significant unknown posterior approximation errors,particularly for BNP models with heavy tails that support the power-lawbehavior of the population. Previous work on truncated inference in BNPmodels has determined the truncation level via analysis of the forward gen-erative model, which does not accurately reflect the error of approximationof the target posterior distribution. This thesis aims to develop approximateinference algorithms that can be directly used for posterior inference forgeneral BNP models. We propose truncated inference methods and provideestimates of the posterior truncation error. Rather than setting the trunca-tion level based on prior approximation error, we establish a desired posteriortruncation error level, allowing the algorithm to adapt the truncation leveluntil the desired truncation error is reached. The proposed algorithms aregeneral in that they can be applied to a wide range of BNP models withcompletely random measure (CRM) priors. We have applied these algorithmsto edge-exchangeable network models, where feature assignment variablesare observed, and to latent feature models with latent feature assignmentvariables.iiLay SummaryBayesian nonparametric (BNP) models provide a flexible and powerful frame-work for statistical modeling by allowing the complexity of the model togrow with the data volume. However, a BNP model contains infinite numberof parameters, making it generally intractable in real applications. Thisthesis aims to develop generally applicable algorithms for BNP models bysubstituting the original model with a tractable approximate model havingonly finite parameters. The approximation error incurred is mathematicallyquantified and numerically estimated in applications. Within this framework,this thesis: (a) proposes one BNP model for network data; (b) develops thealgorithm for network data analysis using the proposed model; (c) derives atheoretically sound algorithm for learning underlying features of data usingBNP models.iiiPrefaceThe thesis is an original intellectual work of the author, Xinglong Li, withthe guidance and supervision of Prof. Trevor Campbell. The researchtopics and the new methods proposed were discussed with the supervisor.Throughout the preparation of this thesis, the supervisor made suggestions onimprovement of structures, and technical details of the thesis. The supervisoralso helped revising the first draft. Mathematical proofs of theoretical resultsin this thesis, as well as numerical experiments, were mainly completed bythe author.Chapter 3 and Chapter 4 are mainly based on a published paper: Li andCampbell. (2021) “Truncated simulation and inference in edge-exchangeablenetworks” in Electronic Journal of Statistics. The supervisors identified theresearch problem and suggested the idea of algorithms. Under guidance ofthe supervisors, the author derived theoretical results, conducted numericalexperiments, analyzed results, and drafted the manuscript.Chapter 5 is motivated as an generalization and extension of the work inChapter 4. The supervisor suggested the idea. The author conducted theresearch work and analyzed the results.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiNomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xivDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.1 Poisson processes . . . . . . . . . . . . . . . . . . . . . . . . 92.1.1 Properties of Poisson processes . . . . . . . . . . . . . 92.2 Completely random measures . . . . . . . . . . . . . . . . . . 132.2.1 Normalized CRMs . . . . . . . . . . . . . . . . . . . . 162.3 Bayesian nonparametric models based on CRMs . . . . . . . 182.3.1 BNP models based on CRMs . . . . . . . . . . . . . . 182.3.2 BNP models based on NCRMs . . . . . . . . . . . . . 232.3.3 Exchangeable partitions probability functions inducedby an NCRM . . . . . . . . . . . . . . . . . . . . . . 242.3.4 Power law behavior . . . . . . . . . . . . . . . . . . . 262.3.5 Hierarchical BNP models . . . . . . . . . . . . . . . . 272.4 Series representations of CRM . . . . . . . . . . . . . . . . . 292.5 Equivalence between inverse Le´vy and rejection representa-tions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32v2.6 Inference for BNP models . . . . . . . . . . . . . . . . . . . . 332.6.1 Exact posterior inference . . . . . . . . . . . . . . . . 342.6.2 Truncated inference . . . . . . . . . . . . . . . . . . . 372.7 Truncation errors . . . . . . . . . . . . . . . . . . . . . . . . 393 Self-product CRM and edge-exchangeable networks . . . . 413.1 BNP on transformed CRMs with observed feature assignments 423.2 Self-product CRM and edge-exchangeable network models . 433.2.1 Edge-exchangeable network model . . . . . . . . . . . 443.3 Truncation of edge-exchangeable network models . . . . . . . 473.4 Prior truncation error . . . . . . . . . . . . . . . . . . . . . . 483.5 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 523.5.1 Beta-Bernoulli process network . . . . . . . . . . . . . 533.5.2 Gamma-independent Poisson network . . . . . . . . . 543.5.3 Normalized-gamma process network . . . . . . . . . . 564 Truncated inference for BNP models with observed labels 584.1 Posterior truncation error . . . . . . . . . . . . . . . . . . . . 594.2 Adaptive truncated Metropolis–Hastings . . . . . . . . . . . 634.3 Posterior inference experiments . . . . . . . . . . . . . . . . . 644.3.1 Synthetic data . . . . . . . . . . . . . . . . . . . . . . 644.3.2 Real network data . . . . . . . . . . . . . . . . . . . . 665 Truncated inference for BNP models with latent labels . . 695.1 Truncated BNP models with latent feature assignment . . . 695.2 Approximate inference via truncation . . . . . . . . . . . . . 725.3 Posterior truncation error . . . . . . . . . . . . . . . . . . . . 735.3.1 Approximation error between arbitrary distributions . 735.3.2 Approximation error for latent feature BNP models . 765.4 Adaptively truncated inference . . . . . . . . . . . . . . . . . 795.5 Adjusted Gibbs sampling algorithms . . . . . . . . . . . . . . 796 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90AppendicesA Proofs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93A.1 Equivalence between rejection and inverse Le´vy representa-tions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93A.2 Truncation error bounds for self-product measures . . . . . . 96viA.2.1 Proof of Theorem 3.4.5 . . . . . . . . . . . . . . . . . 100A.3 Error bounds for edge-exchangeable networks . . . . . . . . . 103A.3.1 Rejection representation . . . . . . . . . . . . . . . . 103A.3.2 Beta-independent Bernoulli process network . . . . . 105A.3.3 Gamma-independent Poisson network . . . . . . . . . 111A.4 Errors of truncated BNP model with latent features . . . . . 115Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120viiList of Figures2.1 The marking and the thinning of a Poisson process. In subplot(a), the points on the horizontal axis represent the originalPoisson process, and the points in the two dimensional spacerepresent the marked Poisson process. In subplot (b), somepoints in the original Poisson process are removed duringthinning and these solid black dots left on the horizontal axisrepresent the thinned Poisson process. . . . . . . . . . . . . . 122.2 Marking and thinning of Poisson process . . . . . . . . . . . . 202.3 The tail behavior of CRP and Pitman-Yor process. The x-axisis log10(N). Each line represent the trace of 1 of 20 simulations 272.4 Graphical model representation of BNP models and hierarchi-cal BNP models. . . . . . . . . . . . . . . . . . . . . . . . . . 282.5 The generative process of a truncated BNP model. . . . . . . 383.1 The generative process of edge-exchangeable network models 463.2 An illustration of truncation of self-product CRMs with d = 2,and the edge-exchangeable network built on it. . . . . . . . . 473.3 Beta-independent Bernoulli network . . . . . . . . . . . . . . 533.4 Gamma-independent Poisson graph . . . . . . . . . . . . . . . 553.5 Normalized-gamma process graph . . . . . . . . . . . . . . . . 574.1 An illustration of the posterior truncation error . . . . . . . . 624.2 Posterior of α, γ, λ for the dense simulated network. Orangehistograms depict the first adaptation iteration, and bluehistograms depict the second and final iteration. . . . . . . . 654.3 Posterior of α, γ, λ for the sparse simulated network. Orangehistograms depict the first adaptation iteration, and bluehistograms depict the second and final iteration. . . . . . . . 65viii4.4 Visualization of the truncation expansion procedure. Blackcircles denote the truncation level and error in each iterationof the adaptation. The adaptation stops when the truncationerror falls below the desired threshold (here log10(0.01) = −2).Grey dashed lines and circles visualize the predicted truncationerror using rates generated from the sequential representation.The vertical dotted line shows that the algorithm selects avalue of K that attempts to match the desired log10 errorthreshold of −2 using the predictions. . . . . . . . . . . . . . 664.5 Posterior of α, γ, λ in the Facebook-like network . . . . . . . 674.6 (4.6a,4.6b): The characteristics of the observed Facebook-likenetwork (black) and 100 samples from the posterior predictivedistribution (blue). (4.6c): The truncation adaptation process,with truncation level and error bound in each iteration (blackcircles) and predictive truncation errors (gray dashed lines). . 685.1 An illustration of the posterior truncation error . . . . . . . . 745.2 Sample means of log10(θk) for samples generated from theregular Gibbs sampler. Each curve represents a trace oflog10(θk) for one truncation level K. Samples of differenttruncation levels are generated independently from the samemodel. Subplots (a) and (b) are samples from the dense modelwith α = 0, and the subplots (c) and (d) are samples from thesparse model with α = 0.3. Subplots (b) and (d) are magnifiedparts of the subplots (a) and (c) respectively, for the purposeto show the pattern of samples of θk’s clearly. It can be seenfrom the plots that sample means of log10(θk) for samplesgenerated from the standard Gibbs sampler, demonstratingthe problematic flat regions that appear to due the inabilityof the sampler to reorder the columns. . . . . . . . . . . . . . 82ix5.3 Sample means of log10(θk) for samples generated from theGibbs sampler with column swaps. Each curve representsa trace of log10(θk) for one truncation level K. Samples ofdifferent truncation levels are generated independently fromthe same model. The black dotted lines are curves of thetrue values of log10(θk) used to generate the synthetic data.Subplots (a) and (b) are samples from the dense model withα = 0, and the subplots (c) and (d) are samples from the sparsemodel with α = 0.3. Subplots (b) and (d) are magnified partsof the subplots (a) and (c) respectively. It can be seen thatthese flat regions in the original Gibbs sampler is alleviateda lot by these two additional column swap steps. However,the estimates of the θk’s are still not close enough to the truevalues, implying that there might be multimodal issue due tothe high dimension of the truncated model. . . . . . . . . . . 845.4 Sample means of log10(θk) for samples generated from theGibbs sampler with column swaps and warm-up stage witha fixed annealing scheme. Each curve represents a trace oflog10(θk) for one truncation level K. Samples of differenttruncation levels are generated independently from the samemodel. The black dotted lines are curves of the true valuesof log10(θk) used to generate the synthetic data. Subplots (a)and (b) are samples from the dense model with α = 0, andthe subplots (c) and (d) are samples from the sparse modelwith α = 0.3. Subplots (b) and (d) are magnified parts ofthe subplots (a) and (c) respectively. It can be seen that theannealed stage generates significant improvement than theGibbs algorithm with feature swap steps. . . . . . . . . . . . 865.5 Sample mean of log10(θK) varying with the truncation levelK, for both the dense model with α = 0 and the sparse modelwith α = 0.3. Samples are generated from the Gibbs samplerwith column swaps and warm-up stage with a fixed annealingscheme. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 885.6 TV error bound of the truncated model varying with thetruncation level K. Samples are generated from the Gibbssampler with column swaps and warm-up stage with a fixedannealing scheme. The subplot (b) is a magnified part ofsubplot (a) where the error bound is smaller than 1. Theshaded regions are 95% Bootstrap confidence intervals. . . . . 88x5.7 TV error bound of the truncated model varying with thetruncation level K. Samples are generated from the Gibbssampler with column swaps and warm-up stage with a fixedquadratic annealing scheme. The shaded regions (which isalmost invisible) is the 95% Bootstrap confidence intervals.When K = 1000, log10(DTV(Πˆµ,Πµ)) ≈ 0.1. . . . . . . . . . 89xiNomenclatureG Discrete distribution of feature assignment on the space N ∪ {0}.g Probability mass function of the feature assignment distributionG.H Diffusive distribution of features on the feature space.h Probability density distribution of the feature distribution H.Z Random variable of feature assignment.BP Beta process.Λ Poisson process.λ Intensity measure of Poisson process.ΓP Gamma process.LP Likelihood process.P Probability of an event.R+ Non-negative real coordinate space.S State space.S The Borel σ-field of S.µ Completely random measure.µ˜ Normalized completely random measure.µK Truncated completely random measure with truncation level K.µ˜K Normalized truncated random measure with truncation level K.ν Intensity (Le´vy) measure of the completely random measure.xiiΨ Feature space of a BNP model.ξ Transformation of a completely random measure.Θ Population rate of features.Πµ Posterior distribution of the BNP model.πµ Unnormalized density function of Πµ.Πˆµ Posterior distribution of the truncated BNP model.πˆµ Unnormalized density function of Πˆµ.xiiiAcknowledgementsA Ph.D. in statistics is not short of randomness. I would not have reachedthis point without the support and love of many people. First, I express mydeepest gratitude to my supervisor, Trevor Campbell, for patiently guidingme in exploring the field of Bayesian statistics. None of this work wouldhave been possible without his detailed guidance. I also want to thankProfessor Jiahua Chen, who has been kindly supporting and caring for mesince before I started my Ph.D. journey. I am grateful to Professor JamesZidek. Working with him in the last year of my Ph.D. has been a fortunateexperience. His understanding, trust, and full support in all aspects trulywarmed my heart. His humility and lasting curiosity about progresses inscience has greatly inspired me. I thank Dr. Kevin Murphy for offeringme the opportunity to contribute to his renowned book. I also extend mygratitude to my supervisory committee members, Alexandre Bouchard-Coˆte´and Ben Bloem-Reddy, for their valuable time reviewing my thesis.Teaching was a significant part of my Ph.D. journey. I thank the instruc-tors in our department, especially Ben Bur, for their support in my teachingendeavors. Additionally, I am grateful to all past and current staff membersin our department, particularly Ali Hauschildt, Peggy Ng, Andrea Sollberger,Kevin Lin, Mairea´d Roche, and Ah Young Kim, for their help.I also want to thank my friends, especially Yichen Zhang, Lin Zhang,David Xu, Chloe You, Naitong Chen, Archer Zhang, Julian Ho, and HaodiLiang, for their encouragement and companionship, particularly during thechallenging periods of my journey. Lastly, I want to thank my ex-girlfriend,Jia Sheng. Her love has been a source of strength for me in confronting theuncertainties of life.xivDedication献给我的家人。To my family.xvChapter 1IntroductionStatistical inference is the process of using data to understand and quantifyuncertainty in its underlying patterns and structures. This goal is usuallyachieved by selecting an appropriate model class, and then updating one’sknowledge of a set of parameters within the selected model class using thedata. For example, in biology, after collecting a set of samples from anunknown population of individuals belonging to some species, scientists areinterested in estimating the population frequency of the observed species,and predicting the amount of observed and unobserved species that will beseen in future samples. These problems are referred to as species samplingproblems and can be formulated as statistical problems related to estimatingprobabilities of subgroups of a population (Fisher et al., 1943; Efron andThisted, 1976). As another example from the field of information retrieval,researchers are often curious about the topics to which each document in acorpus belongs. It is common to view words in a document as arising froma latent number of clusters where a topic can be modeled as a multinomialdistribution of words (Blei et al., 2003a). Depending on the problem, adocument might be classified according to its main topic, or to all topicscovered. The first case is referred to as clustering analysis in statistics, and thesecond case is referred to as factor analysis. Additionally, statistical modelshave also been developed to allow one observation to have different levels ofbelonging to multiple groups. One scenario where such models can be appliedto is in social network data analysis, where one participant can send multiplemessages to many groups. Given a set of observations, people may not onlywant to cluster observations into subgroups but also desire to characterizefeatures underlying each subgroup as well. For example, in genetics, scientistsare often interested in clustering subjects into subgroups according to theirgene expression sequencing data. In addition, another common problemis to characterize latent factors underlying the DNA sequences of differentsubpopulations. The latter problem is usually formulated as statisticalproblems involving learning parameters of mixture models or analyzinglatent factors.Model selection is a key step in all of the above example scenarios. The1model must be of proper complexity to avoid over-fitting or under-fitting thedataset, but also for the discovery of meaningful structure underlying thedata. When fitting a mixture model in a clustering analysis, for example,one must determine an appropriate number of mixture components beforefitting the model and estimating parameters. The same problem exists inlatent factor analysis, where the amount of latent factors should be selectedcarefully.However, model selection is a challenging and computationally demandingtask in general. Depending on the type of problem and statistical frameworkbeing used, it can take a number of different forms. In the case of a predictivemodelling exercise in the context of machine learning, model selection usuallyinvolves fitting the given dataset with models of different complexities andchoosing the model with the highest estimated prediction quality via crossvalidation (Arlot and Celisse, 2010). When working with Bayesian models,one often selects the model with the largest data-marginal probability (Gel-man et al., 2013). In any case, model comparison in traditional parametricmodels is usually constrained to a limited number of candidates chosen apriori, and the model pool must be reset and the model comparison proceduremust be rerun when data changes. However, from information retrieval torecommender systems, from bioinformatics to financial market analysis, theamount of data available has exploded in recent years. As the amount ofdata in a particular dataset grows, the data typically reveal more nuancedlatent structure. In a clustering scenario, the data will exhibit rarer and rarerclusters; in a latent feature assignment scenario, the data will exhibit rarerand rarer features; and so on. Given the challenges of model selection andthe ever-expanding volume of available data in many scenarios, it is desirableto have statistical models that can adaptively adjust model complexity asthe dataset grows.Bayesian nonparametric (BNP) models offer one general framework toaddress this challenge. BNP models have an infinite-dimensional latentparameter space, which can be used to characterize arbitrarily complexpopulation structures. But for any given finite dataset, BNP models assertthat only finitely many of the parameters are involved in the data generatingprocess. This subset of parameters is assumed to grow as the size of the dataincreases.One cornerstone of BNP models is the Dirichlet process (DP) and relatedextensions (Ferguson, 1973). Draws from a DP are random discrete proba-bility measures that contain a countably infinite number of atoms. DPs aretypically used as nonparametric priors for the latent mixture parameters andmixing weights in mixture models used for density approximation or for clus-2tering analysis. The resulting Dirichlet process mixture models (Antoniak,1974) generalize parametric mixture models to contain a countably infinitenumber of components; these models therefore do not require the number ofmixture components to be set a priori, but instead assert that the numberof expressed components grows with the number of observations. The DPand its extensions that allow for power-law behavior (e.g., the Pitman-Yorprocess (Perman et al., 1992; Pitman and Yor, 1997)) of the population havebeen successfully applied to species sampling problems (Pitman, 1996; Lijoiet al., 2007a) and clustering analysis in other fields (Pitman, 2006).BNP models have also been broadly applied in latent factor analysis.In the scenario of feature allocation, where each data point can belongto multiple groups with each belonging indicated by a binary indicator,the Indian buffet process (IBP) (Griffiths and Ghahramani, 2011) and itsheavy-tailed extension (Teh and Gorur, 2009; Broderick et al., 2012) offer anonparametric generalization of latent factor analysis models. While eachobservation can belong to at most a finite number of groups and express afinite number of features, the IBP model allows the total number of groupsto grow with the size of the dataset.In other scenarios, we want individual observations to have differentlevels of belonging to multiple groups or express multiple features to differentextents. BNP models dealing with such problems include the beta-negativebinomial process (Broderick et al., 2014; Zhou et al., 2012), gamma-Poissonprocess (Titsias, 2007), and more general combinatorial trait allocationmodels (Campbell et al., 2018).Another important family of BNP models is hierarchical BNP models,which involve composing multiple BNP priors together to expose layeredlatent structure in data. For example, in a clustering scenario, as the datasetexpands it is often possible to detect more nuanced subclusters within eachcluster in the larger population. The hierarchical Dirichlet process (Teh et al.,2006; Teh and Jordan, 2010), hierarchical beta process (Thibaux and Jordan,2007), and other hierarchical BNP models (Lijoi et al., 2014; Argiento et al.,2019) have successful applications in hierarchical clustering and latent factoranalysis.In statistical modeling, one common basic assumption about the datasetis that observations are exchangeable, meaning that the joint probabilityof the dataset does not depend on the order of arrival of each observation.De Finetti’s theorem (De Finetti, 1929) states, roughly, that exchangeableobservations are conditionally i.i.d. distributed given some random object.The distribution of this random object is usually referred to as the de Finettimixing distribution. It turns out that many popular BNP models can be3formulated as marginal distributions of observations where the de Finettimixing distributions are completely random measures (CRM) (Kingman,1967), or normalized completely random measures (NCRM). For example, ithas been shown that the Beta process (Hjort, 1990) is the de Finetti mixingdistribution underlying the IBP. Given the beta process, feature allocationsof each observation are i.i.d. distributed according to a Bernoulli process(Thibaux and Jordan, 2007). The DP itself can be obtained by normalizingthe gamma process (Ferguson, 1973), and at the same time, it is the de Finettimixing distribution of the Chinese restaurant process (CRP) (Blackwell andMacQueen, 1973), a popular model for infinitely-exchangeable partitionsof a countable set. Given the DP, the observations are conditionally i.i.d.following a multinomial distribution with infinite classes with class weightsspecified by the DP.CRMs provide a general framework within which one can build andanalyze various BNP models; specific models with different properties can beobtained by changing various parameters of the underlying CRM/NCRM, orby choosing different likelihood processes conditioned on the CRM/NCRM.Past work has used this generality of CRMs extensively to study the generalfeatures of BNP models (Regazzini et al., 2003; Thibaux, 2008; Camerlenghiet al., 2019; Campbell et al., 2019; Lee et al., 2023), and to develop posteriorinference methods that are not constrained to specific models (James et al.,2009; Griffin and Walker, 2011; Lijoi et al., 2014; Broderick et al., 2018; Zhuet al., 2020). In this work, we take a similar approach, obtaining results andmethods that generalize to a broad class of BNP models based on CRMs.One major challenge with BNP models that this thesis addresses iscomputational posterior inference. While flexible and expressive, the posteriorinference for BNP models is notoriously hard, for the simple reason that itis impossible to process infinitely many parameters in practice. The lack ofeasily applicable inference algorithms is one of the main factors that impedethe application of BNP models in practice. Since by statistical modeling weaim to infer the structures and patterns underlying a given dataset, and BNPmodels only invoke a finite subset of infinite parameters given finitely manyobservations, posterior inference for BNP models only concerns the finitesubset of parameters expressed by the given dataset. To exclude the infinitetail of inactive parameters from consideration during inference, existingmethods either marginalize them out exactly, or remove them via exact orapproximate truncation.In one family of inference methods, sometimes referred to as collapsedinference methods (Neal, 2000; Ishwaran and James, 2001; Griffiths andGhahramani, 2011), the CRM is marginalized out, and inference is performed4over the posterior distribution of combinatorial variables such as cluster labels(Neal, 2000) or feature/trait assignments (Griffiths and Ghahramani, 2011).While these methods provide exact inferential results, they require that thelikelihood process for cluster/feature assignments be conjugate to the CRMor NCRM, thus limiting the choice of available models. Another crucialproblem with collapsed inference methods is that the feature assignments ofindividual observations are no longer independent after marginalizing outthe de Finetti mixing distribution. Since all parameters are dependent in thecollapsed posterior, the update of the state of any random variable involvesthe whole dataset, making efficient computation (e.g., via parallelization)difficult to design in general. Prohibitive computational cost often renderscollapsed inference for BNP nonviable for large datasets. For example, thecollapsed Gibbs sampler for the IBP has a running time cubic in the numberof observations (Griffiths and Ghahramani, 2011). Efforts to accelerate thecomputational efficiency of collapsed inference algorithms are usually specificto particular (usually conjugate) models (Doshi-Velez and Ghahramani, 2009;Doshi-Velez et al., 2009).Another approach to exact inference is to include inference for the CRMas well. To avoid having to instantiate or update parameters in an infinite-dimensional parameter space, this approach usually introduces auxiliaryvariables that adjust a finite subset of “active” parameters to be updatedin each iteration during sampling (Teh et al., 2007; Walker, 2007; Griffinand Walker, 2011; Zhu et al., 2020). In this case, the conditionally i.i.d.structure of the observations given the latent CRM is maintained, and so thecomputational cost of this approach is usually linear in the size of the dataset.However, this approach often suffers from slow mixing, and in full generalityinvolves expensive numerical integration during sampling (Zhu et al., 2020).There are also exact inference methods that combine both collapsed andconditional sampling and perform posterior inference in a partially collapsedmanner Dubey et al. (2020).In practice, probably the most straightforward approach to posteriorinference is to truncate the BNP model into a parametric one that has onlyfinite parameters and use this parametric model as a surrogate to the trueposterior in performing inference. The truncated model is purely parametric,and inference methods for parametric Bayesian models can be easily appliedto the truncated model. This approach requires a series representationof the CRM (Ferguson and Klass, 1972; Bondesson, 1982; Rosinski, 1990;Rosin´ski, 2001), i.e., a generative description of the CRM involving generatingone atom at a time, which can be terminated in finitely many steps toproduce a computationally feasible parametric approximation. The first5series representations were referred to as “stick-breaking representations,”and were used in designing inference methods for Dirichlet and beta process-related models (Sethuraman, 1994; Ishwaran and James, 2001; Teh et al.,2007), but general series representations have been developed that apply tothe broad class of CRM-based models (Rosin´ski, 2001). Throughout thisthesis, we will consider truncation of general CRM series representations.As an approximate distribution, the truncated CRM induces truncationerrors, and it is essential to quantify and control the truncation error forthe truncated model to be useful. Too small a trunction will yield too muchapproximation error, while too large a trunction will yield a computationallyexpensive inference procedure. Therefore, one key problem when adoptingthe truncated inference is to quantify the truncation error and then choosethe proper truncation level that balances the approximation quality and thecomputation efficiency.Previous work on truncated CRMs only quantify the error in the forwarddata-generative process (Sethuraman, 1994; Muliere and Tardella, 1998;Ishwaran and James, 2001; Ishwaran and Zarepour, 2002; Zarepour andAl Labadi, 2012; Campbell et al., 2019), by examining error in the CRMitself (Ishwaran and Zarepour, 2000; Doshi et al., 2009; Argiento et al., 2016)Truncation methods include sorting the weights of the CRM from posteriorsamples and keeping only large ones (Caron and Fox, 2017; Zhang andDassios, 2023), or setting the truncation level by monitoring the effectivesample sizes of truncated models at different truncation levels (Griffin, 2016).However, for the goal of performing posterior inference for the BNP model,the real important quantity is the approximation error of the truncatedposterior distribution from the true posterior distribution of the BNP model.The truncation level of the truncated model should be determined based onthe posterior truncation error, rather than the approximation quality of theprior model, which does not contain any information about the dataset.The contributions of this thesis are as follows. We provide the first analysisof the a priori/forward generative truncation error in CRM self-products,which have been used previously in edge-exchangeable network models withweighted or unweighted edges (Li and Campbell, 2021). We propose atruncated simulation algorithm for these edge-exchangeable network modelsbased on the prior truncation error analysis. We also provide the first analysisof posterior truncation error in BNP models as a function of the truncationlevel, dataset size, and BNP model properties. We begin by analyzing modelswith observed feature variables, such as the aforementioned self-product edge-exchangeable network models. Using the bounds on posterior truncationerror, we develop an adaptive truncated sampling procedure. Instead of6setting the truncation level in advance or based on prior approximation error,we set a desired level of posterior truncation error. The sampling algorithmstarts from a small truncated model, estimates the truncation error, andadaptively increases the truncation level until the desired error bounds aremet. The truncated posterior distribution in the final step is guaranteedto have small approximation errors to the target posterior distribution ofthe BNP model. We then extend this analysis to BNP models where thefeature assignment variables are latent. Examples include DP mixture modelsand latent factor/trait models. With this result, we propose an adaptivelytruncated sampling algorithm for general BNP models with latent featureassignments built on CRMs, as we did earlier for BNP models with observedfeature assignments.While various metrics can be used to measure the distance (or discrepancy)between two distributions, we quantify posterior approximation error usingthe total variation (TV) distance between the target posterior distributionand its truncated approximation. Theoretically, a small TV distance ensuresthat the probabilities of any event and the expectations of any boundedfunctions are similar under the target and the approximate distribution. Thisstrong theoretical guarantee implies that a truncated model with a low TVapproximation error remains reliable in practice for posterior inference anddecision making. For instance, a small TV approximation error bounds thedifference in quantiles between the true and the truncated model, and limitsthe discrepancy in probability mass assigned to any credible intervals. InBayesian decision theory, it also ensures that posterior expected losses canbe accurately estimated using the truncated model.The thesis is organized as follows. We begin with the background inChapter 2, where we introduce basic concepts and properties of Poissonprocesses, CRMs, and BNP models. We also conduct a brief review of BNPinference algorithms at the end of the chapter. The work in Chapter 3and Chapter 4 are published in Li and Campbell (2021). In Chapter 3 weintroduce a family of edge-exchangeable network models based on seriesrepresentations of BNP models. We propose a truncated simulation methodand provide the analysis of truncation errors of BNP models based on self-product CRMs. In Chapter 4 we analyze the posterior truncation error whichcan be estimated with the truncated model and the prior truncation error.We also propose a truncated inference algorithm for the edge-exchangeablenetwork models, which adaptively increases the truncation level until thedesired approximation error level is reached. Chapter 5 is the unpublishedwork, where we will analyze the approximation error of truncated BNPmodels with latent variables. Chapter 6 is the conclusion and future work.7Chapter 2BackgroundA basic assumption underpinning many statistical models is that the dataare exchangeable, meaning that the joint probability distribution of the datais not affected by the order of the observations. With this basic assumptionof exchangeability, the de Finetti theorem states that the observations areconditionally i.i.d. (independent and identically distributed) given somerandom quantity. We will begin our discussion of BNP models with thisfundamental theorem.Theorem (de Finetti theorem). 1 Let (Xi)i≥1 be an (S,S)-valued stochasticprocess, where (S,S) is Polish. We say (Xi)i≥1 is exchangeable ifP(X1 ∈ S1, . . . , Xn ∈ Sn) = P(Xσ(1) ∈ S1, . . . Xσ(n) ∈ Sn)for any permutation σ of 1, . . . , n, all S1, . . . , Sn ∈ S and any n ≥ 1. The pro-cess is exchangeable if and only if there exists a random measure µ, measurablein the tail σ-field of (Xi)i≥1, such that given µ, X1, X2, . . . are independentand distributed according to µ. Furthermore, let µn(S) =∑ni=1 1S(Xi)/n, ifthe sequence is exchangeable, then the distribution of the random measure µis unique, and µn(S) converges to µ(S) almost surely for each S ∈ S.We refer to the distribution of the underlying random measure µ to bethe de Finetti mixing distribution. It turns out that the de Finetti mixingdistributions underlying many common BNP models are the distributions ofunnormalized and normalized completely random measures (CRMs). There-fore, we focus on BNP models based on CRMs and will introduce CRMsand their properties that affect the properties of BNPs. Since CRMs areclosely related to Poisson point processes, we begin with a brief introductionto Poisson processes and their properties, which will be useful in derivingtheoretical results regarding CRMs and BNP models throughout the thesis.We will then introduce some popular BNP models and their features. Atthe end of this chapter, we will give a brief overview of popular posteriorinference methods for BNP models.1Theorem 1.49 in Schervish (2012).82.1. POISSON PROCESSES2.1 Poisson processesThroughout this thesis, we assume that the state space S is a metric spaceand S is the Borel σ-field B(S). In particular, the singleton {x} ∈ S for allx ∈ S. Let λ be a σ-finite measure on S. A Poisson process (PP) Λ withintensity measure λ is a random measure on S taking the formΛ =κ∑k=1δXk ,where X1, X2, . . . are random elements of S, κ ∈ N ∪∞, and Λ satisfies thefollowing two properties:1. For anym ∈ N disjoint measurable subsets A1, A2, . . . , Am of S, the ran-dom variables Λ(A1),Λ(A2), . . . ,Λ(An) are independent (the completeindependence property of the Poisson process);2. The random counts Λ(A) follows a Poisson distribution Poiss(λ(A))with rate λ(A), where 0 ≤ λ(A) ≤ ∞.We denote a Poisson process with intensity measure λ as Λ ∼ PP(λ).2.1.1 Properties of Poisson processesPoisson processes have many desirable properties that make them popularstatistical models in applications, and these properties are also useful instudying the theoretical aspects of completely random measures and BNPmodels. First, it is important to note that the properties of a Poisson processare preserved under many basic operations. For instance, a mapping fromthe original state space to another space induces a Poisson process on thenew state space.Theorem (Mapping theorem (Kingman, 1992)). Let Λ =∑κk=1 δXk bea Poisson process on S with intensity measure λ, and let T : S → V bemeasurable. Let T (Λ) be a mapping of the PP, Λ, which is defined asT (Λ) :=κ∑k=1δT (Xk).Then T (Λ) is a Poisson process on V with intensity measure T (λ), whereT (λ)(·) := λ(T−1(·)).92.1. POISSON PROCESSESThe mapping theorem indicates that we can obtain a Poisson processby transforming samples from another Poisson process, and the intensitymeasure of the resulting Poisson process corresponds to the transformedrate measure of the previous PP. This result is useful for generating Poissonprocesses with complex intensity measures from simpler ones. As we willsee in Section 2.4, the sequence of rates of a completely random measureis a PP, and its samples can be generated by transforming samples from ahomogeneous Poisson process on R+.The superposition theorem and the restriction theorem states that we stillhave a PP by adding to or cutting on a PP. In particular, the superpositiontheorem states that the sum of any PPs is still a PP with intensity given bythe sum of individual intensities.Theorem (Superposition theorem (Kingman, 1992)). Let Λi, i ∈ N, be asequence of independent Poisson processes on S with intensity measure λi.ThenΛ =∞∑i=1Λiis a Poisson process with intensity measure λ = λ1 + λ2 + . . ..It is often the case that Λi are independent Poisson processes on disjointsubsets of S. As we will discuss in Section 2.2, completely random measurescan be modeled as Poisson processes on S × R+. The superposition theo-rem implies that independent CRMs on disjoint parts of S sum to form acompletely random measure on a larger set.The inverse conclusion also holds. In particular, the restriction theoremstates that a Poisson process on a subset of the state space is itself a Poissonprocess with the intensity measure being the restriction of the originalintensity measure to the subset. Furthermore, restrictions of a Poissonprocess on disjoint subsets are independent Poisson processes.Theorem (Restriction theorem (Last and Penrose, 2017)). Let Λ be aPoisson process on S with σ-finite intensity measure λ and let C1, C2, . . . ∈ Sbe pairwise disjoint. Then ΛC1 ,ΛC2 , . . . are independent Poisson processeswith intensity measures λC1 , λC2 , . . . , respectively, where λCi(A) := λ(Ci ∩A)for any A ∈ S.Again, with the fact that completely random measures are PPs on S×R+,the restriction theorem for Poisson processes implies that the restriction of acompletely random measure to a subspace of S remains a CRM. Furthermore,102.1. POISSON PROCESSESthis theorem also implies that the restriction of a CRM by limiting its weightsto a neighborhood of 0 in R+ is still a CRM. This implies that the tail of aCRM retains the properties of a CRM.One key property that connects the Poisson process with completelyrandom measures is the marking property, which states that a Poissonprocess on a higher dimensional space can be constructed by augmentingthe atoms of a PP on a lower dimensional space with independently drawn“marks.” Conversely, a PP on a lower dimensional space can be obtainedby slicing a PP on a higher dimensional space. In particular, we will see inSection 2.2 that a CRM can be constructed as a marked PP.Specifically, let Λ =∑κn=1 δxk be a Poisson process on S, and K be aprobability kernel from S to V. Let V1, V2, . . . be random elements in V andassume that the conditional distribution of (Vn)n≤κ given (Xn)n≤κ is that ofindependent random variables with distributions K(Xn, ·) for n ≤ κ. Thenthe point processξ :=κ∑n=1δ(Xn,Vn)is called a K−marking of Λ. If there is a probability measure Q on V suchthat K(x, ·) = Q for all x ∈ S, then ξ is called an independent Q-marking ofΛ.Theorem (Marking theorem). 2 Let ξ be a K-marking of a Poisson processΛ with intensity measure λ. Then ξ is a Poisson process on S × V withintensity measure λ×K.The inverse operation of superposition is thinning, which removes atomsof a PP randomly. Let q : S→ [0, 1] be measurable. Consider the probabilitykernel K from S→ {0, 1} defined byKq(x, ·) := (1− q(x))δ0 + q(x)δ1, x ∈ SIf ξ is a Kq-marking of a proper point process Λ, then ξ(· × {1}) is called aq-thinning of Λ. The thinning theorem states that both the thinned atomsand remaining ones of the Poisson process are still Poisson processes.Theorem (Thinning theorem (Last and Penrose, 2017)). Let q : S→ [0, 1]be measurable and let Λq be a q-thinning of a Poisson process Λ. ThenΛq := ξ(· × {1}) and Λ− Λq are independent Poisson processes.2Theorem 5.6 in Last and Penrose (2017).112.1. POISSON PROCESSES(a) Marked Poisson process (b) Thinned Poisson processFigure 2.1: The marking and the thinning of a Poisson process. In subplot(a), the points on the horizontal axis represent the original Poisson process,and the points in the two dimensional space represent the marked Poissonprocess. In subplot (b), some points in the original Poisson process areremoved during thinning and these solid black dots left on the horizontalaxis represent the thinned Poisson process.The marking property allows the construction of a Poisson process on ahigher dimensional space, while the thinning property states that a randomsubset of a PP could also be a PP, and thus can be analyzed using allproperties of PPs. Fig. 2.1a demonstrates the process of marking a PPby accompanying each atom of the PP with independent random variables,resulting in a new PP on the expanded space. Fig. 2.1b illustrates the processof thinning a PP by independently retaining each point mass according to acertain probabilistic distribution 3.An important problem in the statistical inference of a Poisson processis to estimate the expectation of the summation f , a real-valued functionon the state space S, with respect to Λ. Campbell’s theorem states that theexpectation of the sum of a function over a Poisson process equals the integralof that function over the intensity measure of the PP. This equivalence allowsfor the calculation of the expectation, variance, and other statistics of therandom sum.Theorem (Campbell’s Theorem). 4 For a Poisson process Λ ∼ PP(λ) and ameasurable function f : S→ R, the random sum ∑ f(x) is absolutely conver-gent with probability one if and only if the integral∫S min(|f(x)|, 1)λ(dx) <∞.3Chapter 5 in Last and Penrose (2017)4Chapter 3 in Kingman (1992)122.2. COMPLETELY RANDOM MEASURESProvided that this integral is finite, thenE[∫Sf(x)Λ(dx)]=∫Sf(x)λ(dx),andE[exp(∫Sf(x)Λ(dx))]= exp[∫S(ef(x) − 1)λ(dx)].Campbell’s theorem is frequently used in analyzing properties of BNPmodels based on CRMs. As we will explore in Section 2.3, CRMs serve asthe de Finetti mixing distributions for many BNP models. Since CRMs areessentially Poisson processes, the statistical properties of expectations ofBNP samples over the completely random measure are often analyzed usingCampbell’s theorem. For instance, Campbell’s theorem is a widely usedtheoretical tool in studying the prior truncation errors of truncated BNPmodels (Doshi et al., 2009; Campbell et al., 2019).Another useful result for the Poisson process is the Mecke equation, whichprovides a method for computing the expectation of sums with respect to aPoisson process, where the integrand can depend on both the point processand specific points in the state space. In our work, the Mecke equation willbe useful when we analyze the truncation error of truncated BNP networkmodels based on self-products of CRMs in Chapter 3.Theorem (Mecke equation). 5 Let Λ be a Poisson process on S with intensitymeasure λ and let m ∈ N. Denote N as the space of counting measures on S.Then for any f : N× Sm → [0,∞],E[∫f(x1, . . . , xm,Λ\δx1\ . . . \δxm)Λ(m)(d(x1, . . . , xm))]=∫E[f(x1, . . . , xm,Λ)]λ(x1)dx1 . . . λ(xm)dxm.2.2 Completely random measuresCompletely random measures and random probability measures obtained bynormalizing completely random measures serve as de Finetti mixing distribu-tions for many popular BNP models, under the fundamental assumption thatobservations are exchangeable. In this section, we introduce the concept of a5Theorem 4.5 in Last and Penrose (2017).132.2. COMPLETELY RANDOM MEASUREScompletely random measure and its properties that are crucial in analyzingBNP models.A completely random measure (CRM) µ on a measurable space (S,S) isa random measure such that for any collection of K ∈ N disjoint measurablesets A1, . . . , AK ⊂ S, the values µ(A1), . . . , µ(AK) are mutually independentrandom variables (Kingman, 1967).When the state space is R, the completely random measure is known asa random measure with independent increments (Kingman, 1975), or Le´vyprocess. The concept of a completely random measure extends the idea of arandom measure with independent increments to arbitrary state spaces. Ithas been shown (Kingman, 1967) that any completely random measure µcan be decomposed into three components:µ = µdet + µfix + µord = µdet +J∑j=1Vjδψj +κ∑k=1θkδψk ,where µdet is a deterministic measure; µfix is called a fixed measure with{ψ1, . . . , ψJ} being fixed components of S. The collections may be countablyinfinite and Vj are mutually independent random variables on R+. The termµord is referred to as the ordinary component, and∑κk=1 δ(θk,ψk) is a Poissonprocess on R+ × S. The principle of complete randomness within CRMsstipulates that there should be no random diffusive component in any CRM,ensuring the independence of increments across disjoint sets.For the purpose of statistical inference, the deterministic component µdetis not of primary interest, and µfix is easy to handle with standard methods.Throughout this thesis, we focus on CRMs that contain only the ordinarycomponent, which we use as BNP priors. We denote it asµ =κ∑kθkδψk ∼ CRM(ν,H).We refer to the ψk as features and the θk as rates of features. Even thoughwe consider only the ordinary component of the CRM when selecting priorsfor BNP models, the posterior distribution of µ given samples drawn from µ,is a CRM that includes both a fixed component and an ordinary component(Broderick et al., 2018). In the posterior distribution, these samples serve asthe features in the fixed component.Typically µ is constructed as a Poisson process on R+ ×Ψ with intensitymeasure ν(dθ) × H(dψ), where ν is commonly referred to as the Le´vymeasure on R+ and H is a probability distribution function on S. Specifically,142.2. COMPLETELY RANDOM MEASURES∑κk=1 δθk is a PP on R+ with intensity measure ν and (ψ1, . . . , ψκ)i.i.d.∼ Hare marks at (θ1, . . . , θκ). Treating µ as constructed from a marked Poissonprocess allows us to impose assumptions on the intensity measure ν on R+and H on S to make sure that µ exhibits desired properties as a CRM. Inparticular, we assume that ν is a σ-finite measure satisfying the followingconditions:ν(R+) =∞ and∫R+min(1, θ)ν(dθ) <∞. (2.1)The assumption ν(R+) = ∞ guarantees that κ = ∞ almost surely in thePoisson process∑κk=1 δψk . Together with the infinite mass assumption, thesecond assumption in Eq. (2.1) implies that there are at most finitely manypoints of θk falling in [1,∞) and that µ(S) < ∞ almost surely. We alsoassume that H(dψ) is a diffusive distribution so that there are no overlappedmass points in µ. The assumption detailed in Eq. (2.1) is also employed inRegazzini et al. (2003) for constructing a random measure with independentincrements on R+, which can be normalized to construct a random probabilitymeasure.Two commonly used CRMs as priors for BNP models in statisticalinference and machine learning are the beta process and the gamma process.The Le´vy measures of both of them satisfy the conditions in Eq. (2.1) withinthe range of their hyperparameters.Beta processThe beta process BP(γ, λ, α) (Hjort, 1990; Thibaux and Jordan, 2007) withthe discount parameter α ∈ [0, 1), the concentration parameter λ > −α, andthe mass parameter γ > 0 is a CRM with rate measureνB(dθ) = γΓ(λ+ 1)Γ(1− α)Γ(λ+ α)θ−1−α(1− θ)λ+α−1dθ, θ ∈ (0, 1). (2.2)The beta process generates random measures µ with rates θk ∈ [0, 1]. As afunction of θ, the Le´vy measure νB(dθ) is a degenerate Beta distribution,which justifies its name. This process is conjugate to the Bernoulli likelihoodprocess (Thibaux and Jordan, 2007), and the joint beta-Bernoulli model iswidely used in BNP models for feature allocations (Broderick et al., 2012).Additionally, another likelihood that is conjugate to the beta process is thenegative-binomial distribution (Broderick et al., 2014). The beta-negative-binomial model is a preferred choice for trait allocations, accommodatingdifferent levels of belonging of each observation to different groups (Broderick152.2. COMPLETELY RANDOM MEASURESet al., 2014). We will discuss more about the beta-Bernoulli process inSection 2.3.Gamma processThe gamma process ΓP(γ, λ, α) (Brix, 1999; Titsias, 2007) with the discountparameter α ∈ [0, 1), scale parameter λ > 0 and mass parameter γ > 0 isthe completely random measure with the rate measureνG(dθ) = γλ1−αΓ(1− α)θ−α−1e−λθdθ, θ > 0. (2.3)The gamma process generates random measures µ with rates θk ∈ [0,∞).In analogy to the beta process, the name gamma process comes from thefact that its Le´vy measure is proportional to the probability mass functionof a degenerate gamma distribution. Note that in some other referencesthe CRM with Le´vy measure in Eq. (2.3) is referred to as the generalizedgamma process (Hougaard, 1986; Brix, 1999), while the term gamma processspecifically refers to the special case where α = 0. Some other special casesof the gamma process are well studied in the literature. In particular, theCRM with a Le´vy measure proportional to Eq. (2.3) but with λ = 0 andσ ∈ (0, 1) is referred to as the stable process (Perman, 1990; Perman et al.,1992). Additionally, the special case where α = 1 in Eq. (2.3) corresponds tothe inverse Gaussian process (Lijoi et al., 2005).The gamma process features a Le´vy measure that is conjugate to thePoisson distribution (Titsias, 2007). This conjugate relationship makes thegamma-Poisson process a popular choice of BNP models for trait allocations.Additionally, it is known that the Dirichlet process is essentially a normalizedversion of the gamma process (Ferguson, 1973), which underscores the utilityof the gamma process in clustering analysis.2.2.1 Normalized CRMsDiscrete random probability measures are widely used in Bayesian modelingfor clustering analysis and for density approximation using mixture models.One method to obtain a random probability measure is by normalizing acompletely random measure. Let µ be a completely random measure on Sthat satisfies the conditions in Eq. (2.1). Given that µ(S) < ∞, one canconstruct a normalized completely random measure (NCRM) by normalizing162.2. COMPLETELY RANDOM MEASURESthe weights:µ =κ∑k=1θkδψk ∼ CRM(ν,H)TS = µ(S)µ˜ =µTS=κ∑k=1pkδψk , (2.4)where pk := θk/TS . We denote the normalized random measure as µ˜ ∼NCRM(ν,H). Note that when S = R, this definition coincides with the con-cept of a normalized random measure with independent increments (Regazziniet al., 2003; Lijoi et al., 2007b; James et al., 2009).Many random probability measures can be represented as normalizedcompletely random measures. For example, the Dirichlet process DP(α,H)with base measure H and concentration parameter α is a normalized gammaprocess µ ∼ ΓP(γ, 1, 0) (Antoniak, 1974).One generalization of NCRMs is called the Poisson-Kingman processintroduced by Pitman (2003). This generalization allows for the total massTS of the CRM to come from a distribution different from the one associatedwith the CRM. Specifically, we denote by CRM(ν,H | t) the random measureobtained by conditioning the CRM on having a total mass µ(S) = t. Let ρtbe a distribution on R+; a Poisson-Kingman process is a random probabilitymeasure constructed as followsTS ∼ ρt(dt)µ | TS ∼ CRM(ν,H | TS)µ˜ =µTS.We note that not all random probability measures can be obtained by nor-malizing a completely random measure; some are defined via other processes.For example, the popular Pitman-Yor process, which is a generalization ofthe Dirichlet process, cannot be obtained by normalizing any CRM, but is aspecial case of the Poisson-Kingman process (Pitman, 2003). In this thesis,we restrict our discussion to BNP models based on CRMs and NCRMs.Nonetheless, we will demonstrate that BNP models based on CRMs aresufficiently rich to model various patterns in observations.172.3. BAYESIAN NONPARAMETRIC MODELS BASED ON CRMS2.3 Bayesian nonparametric models based onCRMsA key application of CRMs and NCRMs is in the construction of Bayesiannonparametric (BNP) models. Many BNP models—and in particular, thoseappearing in this thesis—can be constructed as models with a countablyinfinite-dimensional space, but for which only a finite subset of these parame-ters are invoked to generate a finite dataset, and that finite subset grows withthe number of observations. In these models, the (normalized) completelyrandom measure appears as a de Finetti mixing distribution. For instance,the Dirichlet process, which is a normalized gamma process, serves as thede Finetti mixing distribution underlying the Chinese restaurant process.Another example is the Beta process, which acts as the de Finetti prior forthe Indian buffet process. The connection between BNP models and CRMs,as established by the de Finetti theorem, enables the construction of otherBNP models from CRMs (Lijoi et al., 2010).In this section, we will demonstrate how BNP models for feature allocationand trait allocation can be constructed from CRMs, and how models forclustering analysis or density approximation using mixture models are builtupon NCRMs.2.3.1 BNP models based on CRMsWe consider two types of factor analysis model: the feature allocation modeland the trait allocation model. In the feature allocation model, we assumethat each observation can be assigned to multiple groups (or express multiplefeatures), with each group assignment indicated by a binary variable. Incontrast, the trait allocation model allows each observation to exhibit differentlevels of belonging to various groups, which is typically indicated by annonnegative integer-valued variable.We define Bayesian Nonparametric (BNP) models for factor analysis byoutlining the generative procedures for feature allocations or trait allocations.This generative definition begins with a completely random measure, whichserves as the de Finetti mixing distribution in the BNP model. Given a CRM,each observation is then assigned to features independently and identically,following a stochastic process that determines the presence and weights ofthese features. In particular, we letµ =∞∑k=1θkδψk182.3. BAYESIAN NONPARAMETRIC MODELS BASED ON CRMSbe a sample from a completely random measure CRM(ν,H). Based onassumptions in Eq. (2.1) regarding ν, we know that the discrete measure µhas infinite atoms and that µ(S) =∑∞k=1 θk <∞ almost surely. Note thatµ is not a probability measure and, therefore, cannot directly be used togenerate samples. To generate feature assignments conditioned on µ, we letthe feature assignment for each observation be a completely random measurewith only a fixed component located at {ψ1, ψ2, ...} corresponding to atomsfrom µ. In particular, letXn =∞∑k=1Zn,kδψk , Zn,kindep∼ G(θk), n = 1, 2, . . . , k = 1, 2, . . . .We say Xn follows a likelihood process with the base measure µ and theassignment distribution G, and denote it as Xni.i.d.∼ LP(G,µ). We referto feature assignment variables Zn,k as labels of observations. Fig. 2.2aillustrates the process of feature assignments of a BNP model conditionedon a CRM.Most BNP models choose the feature assignment distribution G from adiscrete distribution on N. By construction,∑κnk=1 δ{θk,Zn,k} forms a markedPoisson process on S× N with an intensity measure ν(dθ)×G(dz | θ). Con-versely, if the probability mass function of G is g(z | θ), then ∑Zn,k=z δ{θk,z}constitutes a thinned Poisson process with intensity measure ν(dθ)g(z|θ).If all these thinned Poisson processes corresponding to Z > 0 have finitelymany atoms in total, then each observation can belong to at most a finitenumber of groups, or express at most a finite number of features. This isensured by assuming that∫R+(1− g(0 | θ))ν(dθ) <∞. (2.5)This requirement aligns with Assumption A2 in Broderick et al. (2018), andwe will see that it is commonly satisfied by prevalent BNP models. Therefore,we will continue to make this assumption in our subsequent discussions.One example of BNP models of this type is the topic model (Blei et al.,2003b; Teh et al., 2006). In these models, each ψk represents a topic in agiven corpus. The variable Zn,k either represents the presence (or absence)of the topic ψk in the document Xn if Zn,k ∈ {0, 1}, or it represents thenumber of words from topic ψk in Xn if Zn,k ∈ N.To summarize, BNP models that use a CRM µ and a distribution G for192.3. BAYESIAN NONPARAMETRIC MODELS BASED ON CRMS(a) Feature assignment of BNP modelsconditioned on a CRM.(b) Feature assignments of BNP modelsconditioned on an NCRM.Figure 2.2: Marking and thinning of Poisson processexplicit feature assignments can be formulated as follows:µ =∞∑k=1θkδψk ∼ CRM(ν,H)Xn | µ =∞∑k=1Zn,kδψki.i.d.∼ LP(G,µ), n = 1, . . . , N.However, in some applications, features and feature assignments arelatent. The latent features ψk can be continuous, as in the latent factoranalysis model (Roweis and Ghahramani, 1999) and the probabilistic principalcomponent analysis model (Tipping and Bishop, 1999), or discrete, as in thecooperative vector quantization model (Zemel and Hinton, 1993; Ghahramani,1994).In this thesis, we consider a straightforward case where we assume thateach observation follows a distribution parameterized by a linear combinationof latent features assigned to this observation. For example, in the classicalelimination by aspects model (Tversky, 1972) in psychology for choice be-havior, people’s choices among multiple options are influenced by multiplefeatures of each option. For example, when choosing which smartphone tobuy, factors such as brand, CPU, operating system, screen size, and cost mayinfluence decisions. In applications where features and feature assignmentsare latent, we assume that each observation Yn follows a distribution Fparameterized by∑k zn,kψk. The BNP model in this case can be formulated202.3. BAYESIAN NONPARAMETRIC MODELS BASED ON CRMSas follows:µ =∞∑k=1θkδψk ∼ CRM(ν,H) (2.6)Xn | µ =∞∑k=1Zn,kδψki.i.d.∼ LP(G,µ), n = 1, . . . , NYn | Xn indep∼ F(κn∑kzn,kψk).Conjugate likelihood processFor feature assignment, we can choose any discrete probability distribution Gon N that is parameterized by θk. However, a commonly preferred choice isto select G whose probability mass function is conjugate to the Le´vy measureof the CRM. Similar to parametric modeling, choosing a likelihood processconjugate to the CRM enhances BNP models by yielding desirable propertiesand elegant posterior distributions that facilitate posterior inference.For example, if the intensity measure ν and the assignment distributionG are both from an exponential family and are conjugate to each other, itis shown (Broderick et al., 2018) that given observations X1, . . . , XN , theposterior distribution of µ comprises a CRM with one fixed location compo-nent and one ordinary component. The locations of this fixed component aredetermined by the active feature assignments in X1, . . . XN , and the intensitymeasure of the posterior ordinary component remains within the same expo-nential family as the prior CRM. This configuration includes hyper-parameterupdate rules analogous to those found in parametric exponential distributionfamilies. Additionally, when likelihood processes are conjugate to the Le´vymeasure of the CRM, it becomes feasible to derive marginal distributions ofobservations, further enhancing the tractability of these models.We will discuss two commonly used BNP models based on conjugateCRMs and their corresponding assignment distributions. These includethe beta-Bernoulli process, which is popular for feature allocation, and thegamma-Poisson process, widely used for trait allocation.Beta-Bernoulli processOne distribution conjugate to the beta distribution is the Bernoulli distribu-tion. The beta-Bernoulli process is the BNP model constructed on the betaprocess as the intensity measure in Eq. (2.2), µ ∼ BP(γ, λ, α), and uses the212.3. BAYESIAN NONPARAMETRIC MODELS BASED ON CRMSBernoulli distribution Bern(θ) as the feature assignment distribution in thelikelihood process.Suppose that µ ∼ BP(γ, λ, α) and that Xn | µ i.i.d.∼ LP(Bern, µ) for n =1, . . . , N , then the marginal distribution of X1:N follows the two-parametersIndian buffet process (Teh and Gorur, 2009; Broderick et al., 2013, 2018).For n = 1, . . . , N , the conditional marginal distributionXn | X1, . . . , Xn−1is determined as follows:1. For active features ψk’s in X1, . . . , Xn−1,Zn,k | X1:n−1 ∼ Bern(−α+∑n−1i=1 zi,kλ+ n− 1)2. Xn has Mn new features not expressed yet by X1:n−1, whereρn = γΓ(λ+ 1)Γ(λ+ α+ n− 1)Γ(λ+ n)Γ(λ+ α)Mn ∼ Poiss(ρn)ψn,ki.i.d.∼ H(dψ), k = 1, . . . ,Mn.In the special case α = 0 and λ = 1, the marginal distribution of X1:N isthe classical Indian buffet process (Griffiths and Ghahramani, 2011). Theconstruction process illustrates that the each observation has a finite numberof active features, yet the total number of active features expressed by allobservations X1:N grows with N without limit.Gamma-Poisson processIt is desirable for models to allow single observations to belong to multiplegroups to different extents or express multiple features at varying levels. Thiscan be achieved using discrete feature assignment distribution G on N inthe likelihood process. Commonly used distributions include the Poissondistribution, which is conjugate to the gamma process, and the negativebinomial distribution, which is conjugate to the beta process.In particular, the gamma-Poisson process utilizes a gamma process asthe prior. Features are then assigned to observations according to a Poissondistribution, which allows for varying counts of feature representation. Thegamma-Poisson process is formulated as follows:µ ∼ CRM(νG, H)Xni.i.d.∼ LP(Poiss, µ), n = 1, . . . , N.222.3. BAYESIAN NONPARAMETRIC MODELS BASED ON CRMS2.3.2 BNP models based on NCRMsIn clustering analysis, each observation is assigned to precisely one group, withthe primary goals often being to analyze the proportions of different groupswithin the population or to predict the labels of unseen observations. BNPmodels generalize traditional parametric clustering models by utilizing aninfinite-dimensional parameter space, thereby allowing the number of groupsto grow with the size of the dataset. BNP models for clustering analysis anddensity estimation can be constructed using normalized completely randommeasures (NCRMs).In particular, a BNP model for an exchangeable sequence of observationsX1, . . . , XN can be constructed generatively as followsµ˜ =∞∑k=1pkδψk ∼ NCRM(ν,H)Xn | µ˜ i.i.d.∼ µ˜, n = 1, . . . , N.In contrast to BNP models designed for feature allocations, there is no needto introduce a separate likelihood process for clustering analysis models.We can directly sample each observation Xn from µ˜, as it is a probabilitymeasure that already specifies a distribution on the state space S. Thisprocess is demonstrated in Fig. 2.2b.In density approximation or clustering analysis using mixture models,the features Xn are latent, and the generative process for the observations isµ˜ =∞∑k=1pkδψk ∼ NCRM(ν,H)Xn | µ˜ i.i.d.∼ µ˜, n = 1, . . . , NYn | Xn indep∼ F (Xn).Dirichlet process and Chinese restaurant processPerhaps the most popular BNP model in clustering analysis and densityestimation is the Dirichlet process and its extensions. It has been shown thatthe Dirichlet process is a normalized Gamma process (Ferguson, 1973) withan intensity measure νG specified in Eq. (2.3) with a discount parameterα = 0. Since the multinomial distribution is conjugate to the gamma process,Ferguson showed that the posterior distribution of a DP remains a DP withan updated parameter and an updated base measure. Additionally, themarginal distribution of X1, . . . , Xn has a closed form.232.3. BAYESIAN NONPARAMETRIC MODELS BASED ON CRMSIn particular, the initial sample X1 ∼ H. Conditioned on X1, . . . , Xn−1,the distribution for Xn either equals an existing feature with probabilityproportional to its frequency in X1, . . . , Xn−1, or it is sampled from H:Xn | X1:n−1 ∼ αn− 1H +n−1∑i=1δxn .This process of generating X1:N from the marginal distribution is referredto as the Po´lya urn scheme. If we disregard the feature values ψk and onlyconsider the partition of 1, . . . , N induced by Z1:N , we obtain the well-knownChinese restaurant process. The total number of clusters grows with theamount of observations. Additionally, the order 1, . . . , N of observationsthat are clustered defines a partition of the sequence 1, . . . , N , with thedistribution over these partitions being equivalent to the one describedabove.2.3.3 Exchangeable partitions probability functions inducedby an NCRMIn species sampling models and certain clustering analysis applications,the focus is primarily on the partition structure of the observations ratherthan the features themselves. For example, in species sampling models, theobjective is to characterize the population intensity of different species andpredict the number of unseen species based on previous samples. In thesecontexts, there is no need to model the features of each species explicitly,and the order in which species appear does not matter.If we label each observation by the order that it is sampled, differentgroups form a partition of the counting sequence [N ] = 1, . . . , N . Wedenote ΠN as the a partition of size n(ΠN ), meaning that ΠN contains ngroups, and let nj represent the number of elements in the jth group of thepartition. Considering that observations are assumed to be exchangeable,the distribution of a random partition should not depend on the order ofthe observations. In other words, a partition should also be exchangeable,meaning that its distribution is invariant under permutations of samples,depending only on the number of sets it contains and sizes of these sets.Specifically, for any permutation σ of [N ], we denote σ(ΠN ) as the partitionachieved by permuting the members of clusters in π. A partition ΠN is saidto be exchangeable ifP(ΠN = πN ) = P(ΠN = σ(πN )).242.3. BAYESIAN NONPARAMETRIC MODELS BASED ON CRMSEquivalently, the distribution of ΠN is such that for each particular partition{A1, . . . , AJ} of {1, . . . , N}, with |Aj | = nj for 1 ≤ j ≤ J , where nj ≥ 1 and∑Jj=1 nj = n:P (ΠN = {A1, . . . , AJ}) = f[n](n1, . . . , nJ),for some symmetric function f[N ] of a sequence of positive integers.In applications, it is a natural requirement that each successive partitionin a sequence of partitions of observational data should encompass all previouspartitions. This requirement is known as projectivity. Specifically, a sequenceof random partitions Π1,Π2 . . . of [1], [2] . . . is said to be projective ifP(J[m](Πn) = πm) = P(Πm = πm), m < n,where J[m](ΠN ) is the partition of [m] obtained by removing all elements ofclusters in ΠN not in [m]. It can be seen that a projective sequence of parti-tions can be constructed by iteratively adding numbers into it. Specifically,f[n](n1, . . . , nk) = f[n+1](n1, . . . , nk, 1) +k∑j=1f[n+1](n1, . . . , nj + 1, . . . , nk).(2.7)A function f that takes a finite sequence of positive integers as arguments iscalled an exchangeable partition probability function (EPPF) if f is symmetricin its arguments and f satisfies the addition rule in Eq. (2.7). It has beenshown by Kingman (1978) that given a projective sequence of exchangeablepartitions Π1,Π2, . . ., there is a unique exchangeable random partition Π∞of N, such that Πn ∼ J[n](Π∞).The discrete nature of DP and other NCRM makes them well suited forplacing priors on EPPF (James et al., 2009; Pitman, 2003, 2006). The EPPFinduced by the NCRM is given in Corollary 6 by Pitman (2003) (see also(James et al., 2009, in Proposition 3)):f(n1, . . . , nk) =(−1)n−kΓ(n)∫ ∞0λn−1e−ψ(λ)k∏k=1ψni(λ)dλwhere ψ(λ) =∫∞0 (1 − e−λx)ν(x)dx is the Laplace exponent and ψm(λ) =(−1)m−1 ∫∞0 xme−λxν(x)dx.252.3. BAYESIAN NONPARAMETRIC MODELS BASED ON CRMS2.3.4 Power law behaviorIt has been observed that observations from many fields demonstrate apower law phenomenon, meaning that the number of subgroups (clusters,features, etc.) in the dataset is asymptotically proportional to the number ofobservations raised to some positive power. For example, word frequenciesin natural languages typically exhibit a power law pattern (Zipf, 2013).Similarly, in the field of image processing, the empirical distribution of labelfrequencies in a database also shows power law behavior (Sudderth andJordan, 2008). In ecology, the frequency of species is well known to followa power law distribution (Gnedin et al., 2007). It is therefore important inBNP models to ensure that the prior process has the capacity to capturethis type of behaviour.However, many traditional BNP models are known to generate manyfewer subgroups in the prior process. For example, samples generated fromthe Dirichlet process or its marginal distribution, exhibit a number of clustersK that grows at a logarithmic pace with sample size, K ∝ log(N). It canbe seen from Section 2.3.4 that there is a linear relationship between K andlog10(N). This logarithmic relationship between the number of features andsample sizes also holds in the Indian buffet process (Griffiths and Ghahramani,2011).This limitation has motivated the development of more general models.One generalization of the DP that supports power-law behavior is the Pitman-Yor process (Perman et al., 1992; Pitman and Yor, 1997). Fig. 2.3 (b) showstwenty simulations from the Pitman-Yor process, where the linear relationshipbetween the logarithm-logarithm transformation clearly indicates a power-lawbehavior relationship between the number of clusters and sample sizes. Thepower-law behavior has made the Pitman-Yor process a popular alternativeto the DP in applications (Goldwater et al., 2005; Teh, 2006; Wood et al.,2009). While the Pitman-Yor process cannot be represented as an NCRM, ithas been shown (Lijoi et al., 2007b) that the expected number of clusters ofsamples from a normalized gamma distribution with a discount parameterα > 0 also grows as a power of the dataset size, similar to the Pitman-Yorprocess. Recall that the Dirichlet process is the normalized gamma processwith a discount parameter α = 0. Therefore, the normalized gamma processsupports both power-law and logarithm-law behaviors by changing the valueof the discount parameter.In feature allocation models, the beta-Bernoulli model with a discountparameter α > 0 also supports power-law behavior (Teh and Gorur, 2009;Hougaard, 1986; Brix, 1999; Broderick et al., 2014). Since the Indian Buffet262.3. BAYESIAN NONPARAMETRIC MODELS BASED ON CRMS(a) Chinese Restaurant Process (b) Pitman-Yor processFigure 2.3: The tail behavior of CRP and Pitman-Yor process. The x-axis islog10(N). Each line represent the trace of 1 of 20 simulationsProcess (IBP) is the marginal distribution of the beta-Bernoulli model witha discount parameter α = 0, it has been shown that the family of beta-Bernoulli processes supports both power-law and logarithm-law phenomenain datasets, and different behaviors are achieved by changing the value of amodel parameter (Broderick et al., 2012).2.3.5 Hierarchical BNP modelsAs the size of a population grows, it is common to see small subgroups withinlarger clusters. It is desirable for models to separate observations into groupswhile simultaneously maintaining connections between these larger clusters.In Bayesian statistics, such a model structure is naturally achieved throughhierarchical modeling. It is assumed that the population is subdivided intosub-populations, and within each sub-population, there are clusters thatcapture the distribution structure of that sub-population.One example is in topic models, where words from each document are clus-tered according to the topics to which they belong, allowing each documentto cover multiple topics. Even though words from different documents do notfollow exactly the same distribution, different documents from one corpusmay share topics. Therefore, in this case, it is desirable to share commoninformation among the distributions of topics when analyzing clusters oftopics in each document. This can be achieved using hierarchical modeling.Hierarchical modeling also provides a solution when the assumptionof exchangeability of observations is violated. For example, in a socialnetwork, the distribution of edges within a social group is exchangeable, but272.3. BAYESIAN NONPARAMETRIC MODELS BASED ON CRMS(a) Representation of BNP models. (b) Representation of HBNP models.Figure 2.4: Graphical model representation of BNP models and hierarchicalBNP models.edges or connections from different groups may not be. According to deFinetti’s theorem, only observations from one cluster are conditionally i.i.d.distributed. Hierarchical modeling addresses this type of data by allowingdifferent clusters to have separate mixing distributions, which are sampledfrom another overarching distribution.The BNP solution to such problems is addressed with hierarchical BNPmodels. As illustrated in Section 2.3.5, the BNP models we have discussedso far have feature assignments generated from µ in an i.i.d. manner. Thehierarchical BNP model is constructed by treating µ as a base measurefor generating a completely random measure µj ∼ CRM(ν, µ), as shown inSection 2.3.5. Since the completely random measure µ is discrete almostsurely, µj share common features ψk ∈ S.Hierarchical clustering analysisThe hierarchical NCRMs are used in clustering analysis when subpopulationsin different groups are clustered following different rules yet the clusters areshared across the entire population. Generally, a hierarchical normalizedCRM can be constructed by first generating a completely random measureµ˜0 and then using it as the base measure to generate an NCRM for each282.4. SERIES REPRESENTATIONS OF CRMsubgroup:µ˜0 ∼ NCRM(ν0, H)µ˜j |µ˜0 i.i.d.∼ NCRM(νj , µ˜0) j = 1, · · · , d.For example, in hierarchical Dirichlet processes (Teh et al., 2006), µ˜0 is a DPwith a diffuse base measure H and µ˜j are DPs with base measure µ˜0. Thisconstruction ensures that one marginally has E[µ˜j |µ˜0] = µ˜0. Conditionedon the NCRMs (µ˜1, · · · , µ˜d), the observations Xj,i within each group areindependent and identically distributed, but only independent across groups.The DP in the hierarchical clustering model could be replaced with othernormalized completely random measures (NCRMs) to support the power lawbehavior of the population (Teh and Jordan, 2010).Just as with the BNP models discussed earlier, we need to marginalize outthe infinite parameters in the hierarchical BNP models to perform posteriorinference. However, closed forms of the predictive distribution only exist forspecific choices of the CRMs used to define the hierarchical structures (Tehet al., 2006).Hierarchical feature allocationThe dependent discrete random measures for feature allocations can beconstructed hierarchically:µ0 ∼ CRM(ν0, H)µj |µ0 ∼ CRM(νi, µ0), j = 1, · · · , d.and each group has its own observation process Xj,i based on its CRM µj :Xj,ii.i.d.∼ LP(g, µj), i = 1, · · · , nj ,where nj is the number of observations from the group ℓ, and LP(G,µj) isthe likelihood process that generate labels.2.4 Series representations of CRMTo simulate a realization µ ∼ CRM(ν,H), the rates θk and features ψk may begenerated sequentially using a series representation (Rosinski, 1990; Rosin´ski,2001) of the CRM. The truncated simulation and inference for CRMs thatappear throughout this thesis are also based heavily on series representations292.4. SERIES REPRESENTATIONS OF CRMof CRMs. In this section we introduce different series representations ofCRMs.As is seen in Section 2.2,∑∞k=1 δ(θk,ψk) is a H-marking Poisson processwith intensity measure ν×H, and the Poisson process Θ =∑∞k=1 δθk ∼ PP(ν).The marking process representation provides a generative process for samplingfrom a CRM. In particular, to sample µ ∼ CRM(ν,H), we can sample Θfirst and then mark Θ by sampling from the distribution H.Specifically, we start by simulating the jumps of a unit-rate homogeneousPoisson process (Γk)∞k=1 on R+ in increasing order. For a given distributiong on R+ and non-negative measurable function τ : R+ × R+ → R+, we setµ =∞∑k=1θkδψk , θk = τ(Uk,Γk), Uki.i.d.∼ g, ψk i.i.d.∼ H. (2.8)The measure µ constructed in this manner is a CRM with a rate measure νthat depends on the choice of τ and g. Different choices of τ and g lead todifferent series representations of one CRM (Campbell et al., 2019).The inverse Le´vy representationThe inverse Le´vy representation (Ferguson and Klass, 1972) takes the formθk = ν←(Γk), ν←(x) = inf {y : ν ([y,∞)) ≤ x} . (2.9)The θk obtained using the inverse Le´vy representation are atoms from thedesired Poisson process with the intensity measure ν, which follows from themapping theorem for the Poisson process in Section 2.1.Evaluation of ν←(x) is usually intractable, making it hard to generate θkin this manner in practice. However, in Section 2.5, we will show that thecompletely random measures simulated using the rejection representationhave the same distribution with those simulated using the inverse Le´vyrepresentation of the same CRM.Inverse Le´vy representation of Beta process The analytical form ofthe inverse Le´vy representation for the beta process can be obtained undercertain choices of the hyperparameters. For example, let νB be the ratemeasure of BP(γ, 1− α, α), then if α > 0,ν[x,∞) = γ′α−1(x−α − 1), ν←(Γk) = (1 + αΓk/γ′)−1/α, (2.10)and if α = 0,ν[x,∞) = −γ′ log(x), ν←(Γk) = e−Γk/γ′ , (2.11)302.4. SERIES REPRESENTATIONS OF CRMwhere γ′ := γ Γ(λ+1)Γ(1−α)Γ(λ+α) . It can be seen from the inverse Le´vy representa-tion of the beta process that θk decreases exponentially when α = 0 while θkdecreases polynomially if α > 0, which demonstrates the power law behavior.Rejection representationWhen the inverse Le´vy representation is not tractable, we can generate aseries of rates from CRM(ν,H) with the rejection representation (Rosin´ski,2001), which takes the formθk = Tk1(dνdνˆ(Tk) ≥ Uk), Tk = νˆ←(Γk), Uk i.i.d.∼ Unif[0, 1], (2.12)where νˆ is a measure on R+ chosen such that dνdνˆ ≤ 1 uniformly and νˆ←(·) iseasy to calculate in closed-form.Note that although the rejection representation is often more tractablethan the inverse Le´vy representation to simulate, there may be some zerorates θk = 0 generated with the rejection representation; this may or maynot be an issue depending on the application.Rejection representation of Beta process To simulate the beta process,we use the rejection representation with a proposal rate measure νˆ given byνˆB(dθ) = γ′1 [θ ≤ 1] θ−1−αdθ, γ′ = γ Γ(λ+ 1)Γ(1− α)Γ(λ+ α) .With this choice of νˆB, we can evaluate νˆ←B anddνBdνˆBin closed form. Specifi-cally, when α = 0,νˆ←B (u) = e−u/γ′ ,dνBdνˆB= (1− θ)λ−1; (2.13)and when α ∈ (0, 1),νˆ←B (u) =(1 +αuγ′)−1/α,dνBdνˆB= (1− θ)λ+α−1. (2.14)Rejection representation of Gamma process For the simulation ofthe gamma process, when α = 0, we use the rejection representation with aproposal rateνˆG(dθ) = γλθ−1(1 + λθ)−1dθ,312.5. EQUIVALENCE BETWEEN INVERSE LE´VY AND REJECTION REPRESENTATIONSfor which we can obtain the closed form of νˆ←G anddνGdνˆG,νˆ←G (u) =1λ(e(γλ)−1u − 1) , dνGdνˆG (1 + λθ)e−λθ. (2.15)When α ∈ (0, 1), we use the rejection representation with the proposal rateνˆG(dθ) = γλ1−αΓ(1− α)θ−1−αdθ,for whichνˆ←G (u)(γ′u−1)1/d, γ′ = γλ1−ααΓ(1− α) ,dνGdνˆG= e−λθ. (2.16)While there are many other sequential representations of CRMs (Campbellet al., 2019), the inverse Le´vy and rejection representations in Eqs. (2.9)and (2.12) are broadly applicable and will play a key role in our theoreticalanalysis.2.5 Equivalence between inverse Le´vy andrejection representationsEvaluation of the inverse Le´vy representation θk = ν←(Γk) is usually in-tractable. In this section we show that the nonzero rates generated fromthe rejection representation (Rosin´ski, 2001) of a Poisson process have thesame distribution as the well-known but typically intractable inverse Le´vyor Ferguson-Klass representation (Ferguson and Klass, 1972). This providesa novel method for simulating the inverse Le´vy representation, which hasa wide variety of uses in applications of Poisson processes. This theoreticalresult also facilitates the study of truncation errors by allowing us to alwaysassume the accessibility of the inverse Le´vy representation.Note that this is a novel contribution of this thesis, instead of backgroundmaterial. We are including it here as its introduction along with the inverseLe´vy and rejection representations is most natural.Theorem 2.5.1. Let µ ∼ CRM(ν,H) be a completely random measure on S.Let θ1, . . . , θK be the first K rates from the inverse Le´vy representation ofµ, as described in Eq. (2.9), and let θ′1, . . . , θ′K be the first K nonzero ratesfrom any rejection representation of µ as outlined in Eq. (2.12). Then(θ1, . . . , θK)d= (θ′1, . . . , θ′K).322.6. INFERENCE FOR BNP MODELSTheorem 2.5.1 implies that we can always use the inverse Le´vy represen-tation when simulating a CRM and in analyzing properties of the truncatedmodel, even when the inverse Le´vy representation does not have a closed form.In particular, we can always use the rejection representation to simulate therates θk, but instead of simulating for iterations k = 1, . . . ,K, we simulateuntil we obtain K nonzero rates. This is no longer a sample of a truncatedrejection representation; but Theorem 2.5.1 states that the first K nonzerorates have the same distribution as simulating K iterations of the inverseLe´vy representation.Even though it is not easy to know in advance which rejection represen-tation rates will be 0, Theorem 2.5.1 is a decent analogue that is worth usingsince the number of zeros is typically finite in the whole process (Campbellet al., 2019) and the assumption of the inverse Le´vy representation simplifiesanalysis quite a bit.2.6 Inference for BNP modelsWhile Bayesian nonparametric models can express structures of unboundedcardinality underlying the data, the fact that BNP models have an infinitenumber of parameters makes posterior inference notoriously difficult. Thelack of efficient inference methods for general BNP models has impeded theapplication of BNP models in practice.However, posterior inference is still in theory viable, since only a finitesubset of parameters of the BNP model are invoked by a given dataset almostsurely. For example, a dataset with N observations can have at most Nclusters. In latent factor analysis, only a finite number of features can beexpressed by a dataset under the assumption specified in Eq. (2.5), whichstates that each observation can express only a finite number of features orbelong to a finite number of groups. Ideally, posterior inference is performedonly over the features that reflect the structure and patterns underlyingthe observed dataset. The infinite inactive parameters that come with theCRM prior and are not related to the dataset must be somehow avoided incomputational posterior inference methods.To avoid processing an infinite number of parameters in the posteriordistribution, previous methods that perform inference for the exact poste-rior either conduct marginal inference—marginalizing out infinite inactiveparameters in the posterior—or introduce auxiliary variables to decide theset of active features in each iteration during posterior sampling.A different approach is to perform inference using approximate distri-332.6. INFERENCE FOR BNP MODELSbutions based on truncated CRMs. The BNP models based on truncatedCRMs are parametric, to which posterior inference algorithms can be applied.However, truncating the BNP model induces approximation errors. To ensurethat the truncated posterior is a good approximation of the true posterior,the approximation error between the truncated posterior distribution andthe true posterior distribution must be quantified and controlled.2.6.1 Exact posterior inferenceCollapsed inferenceAn important family of inference methods for BNP models is the collapsedinference, which integrate out the infinite inactive features of the BNPmodels and perform inference on only the marginal posterior distributionof parameters invoked by the data. For the purpose of posterior inferenceof (latent) feature assignment, the predictive distribution p(Xn|X\Xn) issufficient. The closed form of the marginal distribution usually requires theconjugate relationship between the CRM and the likelihood process.When feature assignment labels Z[N ] can be observed, sampling from theposterior marginal distribution is straightforward with a conjugate CRM andlikelihood process. For example, the Chinese restaurant process describes themarginal predictive distribution of a new sample Zn conditioned on previoussamples Z1:n−1. Under the assumption of exchangeability, any observationcan be treated as the last sample ZN and therefore the CRP scheme providesa Gibbs sampling rule for the posterior distribution. In BNP models forfeature assignment, the Indian buffet process is the marginal distributionof the beta-Bernoulli process. With the assumption of exchangeability, thechoice of dishes for the latest customer provides a Gibbs updating rule forsampling from the posterior marginal distribution. The CRP and IBP shareanalogous posterior updating rules. In particular, in each iteration of Gibbssampling, an observation belongs to an existing group with a probabilityproportional to the size of that group, and to some new groups with aprobability determined by the parameters of the CRM.Broderick et al. (2018) extended the marginal distribution of IBP to BNPmodels with a CRM from an exponential distribution family CRM(ν,H) anda conjugate likelihood process LP(µ,G). For completely random measureswith exponential family distributions and accompanied by the conjugatelikelihood, the posterior distribution of the CRM is in the same familyas the CRM, with analogous parameter update schemes to those in theparametric exponential family distributions. In particular, for n = 2, . . . , N ,342.6. INFERENCE FOR BNP MODELSlet k = 1, . . .Kn−1 be the labels of active features expressed by existingobservations X[n−1]. Then Zn,k has the following conditional distributiong(znk | z[n−1],k) =∫R+ g(znk | θ)∏n−1m=1 g(zmk | θ)ν(dθ)∫R+∏n−1m=1 g(zmk | θ)ν(dθ).Additionally, the observation Xn adds K+n,z new features for each possiblevalues z > 0, where K+n,z ∼ Poiss(ρn,z) withρn,z =∫R+g(z | θ)n−1∏m=1g(0 | θ)ν(dθ).It can be seen that the marginal distribution and the posterior inferenceof feature assignment relies on the integral∫R+ g(z | θ)ν(dθ) having a closedform, which only exists when ν and g are conjugate. Therefore, the collapsedinference for feature assignment is only possible in limited cases, and spe-cial inference algorithms need to be designed for each specific distribution.Furthermore, when Z can take infinite values, it is practically challenging tointroduce new features via K+n,z ∼ Poiss(ρn,z) as feature assignment wouldneed to be introduced for each possible value of Z.It is common for feature assignments Z[N ] to be latent, as in DP mixturemodels and latent feature allocation models based on the IBP. Posteriorinference for BNP models with conjugate CRM priors requires updating thefeature assignments Z[N ], and values of features of active groups when needed. Compared to models with explicit feature assignment, the probability thatone observation belongs to a group is proportional to the size of that groupmultiplied by the likelihood of the observations in that group. (Neal, 2000;Griffiths and Ghahramani, 2011).The marginal distribution is only available for BNP models with a conju-gate CRM and likelihood process. Another problem with this approach isthat Z, . . . , ZN are no longer i.i.d. in the marginal distribution. Therefore, allthe rest of the feature assignments, Z1, . . . , Zn−1, Zn+1, . . . ZN , are involvedin updating the posterior distribution of Zn in each iteration, leading tosignificant computation costs. For example, the computational complexity ofeach iteration of Gibbs sampling for an IBP latent feature model is O(N3)(Griffiths and Ghahramani, 2011), which prohibits its application to largedatasets in practice.352.6. INFERENCE FOR BNP MODELSConditional inference based on slice samplingThe prohibitive computation cost of collapsed inference algorithms and theirrestrictions to conjugate distribution are the main motivations for developingconditional inference algorithms. Conditional inference refers to methodsthat perform posterior inference for both Θ and Z[N ]. One approach to avoiddealing with an infinite number of θk’s is to use slice sampling, which addsauxiliary variables whose values set conditional ranges for the number ofactive features for each observation (Walker, 2007; Zhu et al., 2020).For example, in posterior inference for BNP models based on CRMswith latent feature allocation in 2.6, Zhu et al. (2020) introduces auxiliaryvariables U1, . . . , UN in such a way thatUnindep∼ Unif[0, η(κn)], n = 1, . . . , N,where κn := max{k ∈ N : Zn,k > 0} ∪ {0}, and η(·) is any function such thatη : N→ R+, η is monotone decreasing and limn→∞ η(n) = 0. Conditioned onUn, Xn has only active features in the range 1, . . . , κn such that η(κn) > Un.By construction, the auxiliary variables U1, . . . , UN do not affect either theprior distribution of the likelihood process, or the likelihood of observations.Therefore, the posterior distribution ofX1, . . . , XN obtained by marginalizingout these auxiliary variables is the true posterior distribution.One advantage of using conditional inference is that feature assignmentsare conditionally i.i.d., making the sampling algorithm efficient and suitablefor parallel algorithms. However, this method requires evaluating the proba-bility P(Z[N ],k>K = 0 | θK) in each iteration, whereK is the global truncationlevel in each iteration. For general CRMs, this probability typically cannotbe computed exactly.Partially collapsed inferenceIt is also possible to combine the collapsed inference and the conditionalinference to take advantage of the quick mixing of the collapsed inferencealgorithms and the computational efficiency of conditional inference. This isreferred to as partially collapsed inference and is proposed by Dubey et al.(2020) as an alternative to collapsed inference and conditional inference.The partially collapsed inference algorithm applies to BNP models whoseposterior is a CRM, including an ordinary component as well as a fixedlocation component. The posterior CRM can therefore be decomposed into asummation of a parametric part and a nonparametric part. The parametricpart corresponds to the fixed location component of the posterior CRM,362.6. INFERENCE FOR BNP MODELSwhile the nonparametric part corresponds to the ordinary component of theposterior CRM. Collapsed inference is performed on the nonparametric part,and the conditional inference is performed on the parametric part. A majoradvantage of this methodology is that computation is not wasted on themany empty features that often appear as part of a conditional samplingmethod. However, the collapsed inference for the ordinary component of theposterior CRM still depends on the conjugate relationship between the CRMand the likelihood process.2.6.2 Truncated inferenceA different approach to the posterior inference for BNP models that avoidsdealing with an infinite number of parameters is to truncate the CRMs toa finite sequence that contains only finitely many parameters. Posteriorinference is then performed using the truncated BNP model instead ofthe original model. The idea behind truncated inference is that when thetruncation level is large enough, the truncated posterior distribution is a goodsurrogate for the original posterior since only finite features are invoked by agiven dataset. After truncation, general Bayesian inference methods can beemployed for posterior inference because the truncated model is parametricand contains only finite parameters.The truncated BNP model is obtained by truncating the CRM and build-ing the likelihood process on the truncated CRM, as we do in constructingBNP models on CRMs. In particular, the truncation of a BNP model beginswith a series representation of the CRM µ ∼ CRM(ν,H).µ =∞∑k=1θkδψk , θk = τ(Vk,Γk),where (Γk)∞k=1 are ordered jumps of the unit rate homogeneous Poissonprocess on R+, Vk is some random variable, and τ is a deterministic functionof V and Γ. The choice of the distribution of the random variable Vk andthe function τ depends on the series representation of the CRM, as we seein Section 2.4.The truncation of a CRM µ with truncation level K is obtained bydiscarding all atoms θkδψk for k > K in the series representation of µ. Wedenote the truncated CRM as µK , and the discarded tail as µK+. Based ona series representation, truncation splits µ into two parts as follows:µ = µK + µK+, µK =K∑k=1θkδψk , µK+ =∞∑k=K+1θkδψk .372.6. INFERENCE FOR BNP MODELSFigure 2.5: The generative process of a truncated BNP model.In practice, we simulate µK using a series representation of µ and stop atthe Kth step. The tail part µK+ will be used in analysing the approximationerror between the truncated model and the true model.Fig. 2.5 illustrates the generative definition of the truncated BNP model.Given the truncated CRM µK , we construct BNP models by simulatingfeature assignments from µK in an i.i.d. manner. Specifically, for BNPmodels with explicit feature assignments, the generative process isµ =∞∑k=1θkδψk ∼ CRM(ν,H)µK =K∑k=1θkδψkXn | µK =K∑k=1Zn,kδψki.i.d.∼ LP(G,µK), i = 1, . . . , N.If features and feature assignments are latent, there is an additional layerin the generative process:Yn | Xn indep∼ F(K∑k=1zn,kψk), n = 1, . . . , N.where we assume that the likelihood for observation Yn is parameterized bya linear combination of latent features.Similarly, we obtain a truncated normalized CRM by normalizing thetruncated CRM µK :µ˜K =K∑k=1θkµK(S)δψk =K∑k=1θkTKδψk ,382.7. TRUNCATION ERRORSwhere TK =∑Kk=1 θK . Conditioned on µ˜K , the feature assignments ofobservations can be simulated directly from the discrete finite distributionµ˜K :Xn | µ˜K i.i.d.∼ µ˜K , n = 1, . . . , N.Similarly, if the features and feature assignments are latent, we add anotherlayer to the generative process for the observation model:Yn | Xn indep∼ F (Xn).The truncated BNP model is purely parametric, and therefore generalinference methods for parametric BNP models can be applied for posteriorinference. Additionally, truncated inference only requires a series repre-sentation of the CRM and does not rely on a conjugate likelihood process.Therefore, truncated inference applies to general BNP models and enjoysthe computational efficiency of conditional inference methods. However, thetruncated posterior distribution is not precisely the desired BNP posterior;quantifying the amount of approximation error due to truncation is a primarygoal of this thesis.2.7 Truncation errorsFixed truncation is typically easy to apply to any CRM-based model; onesimply replaces the underlying CRM µ with its truncation µK . However,substituting µ with µK incurs approximation error, since features in the tailhas been removed together with µK , as can be seen from Fig. 2.5. Thereforethe truncated BNP model is only useful when the truncation error can bequantified and is controlled.Previous work on truncated inference for BNP models determine trunca-tion levels either based on the statistics of the truncated model (Ishwaranand James, 2001) or based on the prior truncation error. In particular,researchers have quantified the truncation error in terms of the L1 (i.e.,total variation) distance between marginal data distributions based on theoriginal CRM and the truncated CRM (Sethuraman, 1994; Ishwaran andJames, 2001; Ishwaran and Zarepour, 2002; Blei and Jordan, 2006; Doshiet al., 2009; Paisley et al., 2012; Roy, 2014; Roychowdhury and Kulis, 2015;Campbell et al., 2019). The truncation error is typically a function of theprior process µ, the likelihood process, the level of truncation K, and thenumber of observations N .392.7. TRUNCATION ERRORSMost analyses of truncation error begin by bounding the probability thatdata drawn from the full (nonparametric) model will use a feature that isnot available to the truncated model (Campbell et al., 2019). In particular,for the full and truncated likelihood processXn|µ i.i.d.∼ LP(G,µ) XK,n|µK i.i.d.∼ LP(G,µK),Yn|Xn indep∼ F (·|Xn) Wn|XK,n indep∼ F (·|XK,n),where µ is a CRM and µK is its truncation, if pN,∞ and pN,K are the marginaldensities of Y1:N and W1:N , respectively, then the truncation error can bebounded by Campbell et al. (2019)12∥pN,∞ − pN,K∥ ≤ 1− P(supp(X1:N ) ⊂ supp(µK)). (2.17)This basic inequality is used to derive specific bounds for particular seriesrepresentations and likelihood processes. We will show that the strategyof bounding the truncation error by the probability that there are featuresdiscarded by truncation works in general in other types BNP models. Forexample, the inequality in Eq. (2.17) will be used in Chapter 4 in derivingthe prior truncation error for edge-exchangeable network models.There are two major limitations to this approach. First and foremost,in practice the error in the marginal data distribution is not of interest;in the Bayesian inferential setting, we really care about the error incurredby truncation in the posterior distribution. Second, truncation error hasonly been studied in relatively simple BNP models; for example, past workdoes not apply to models using a wide variety of transformations of CRMs,e.g. the edge-exchangeable network model based on the self-product of CRM.The goal of this work is therefore to develop a truncated inference procedurethat is applicable to general BNP models and provide tractable bounds forthe posterior truncation error.40Chapter 3Self-product CRM andedge-exchangeable networksBayesian nonparametric models have been widely adopted in clusteringanalysis, probability density approximation, feature allocation, and traitallocation. As discussed in Chapter 2, popular BNP models have beenconstructed using completely random measures or normalized completelyrandom measures as priors. Approximate inference methods by truncationfor BNP models based on (normalized) CRMs have been widely adoptedas a generally applicable and computationally efficient alternative to exactposterior inference methods. The truncation error of the truncated models istypically analyzed by quantifying the probability that no features are lost inthe tail region of the CRM, and this probability is usually quantified usingCampbell’s theorem.Previous BNP models based on (normalized) CRMs are commonly con-structed directly on the completely random measures. In this chapter, weextend our analysis to the family of BNP models based on transformations ofCRMs. In particular, we take as an example the self-product of a completelyrandom measure and the family of edge-exchangeable network models basedon such priors (Williamson, 2016; Cai et al., 2016), which is a new field ofapplication for BNP models proposed recently.We will begin this chapter by introducing BNP models based on trans-formations of CRMs, where we also assume in this chapter that featureassignment variables of the BNP models are observed. BNP models withlatent feature assignments will be discussed in Chapter 5. In particular, wewill focus on a specific transformation, the self-product, of a CRM and theedge-exchangeable network models based on it. We propose a truncatedsimulation algorithm for edge-exchangeable network models and analyze theprior truncation error. With the equivalence relationship between the inverseLe´vy representation and the nonzero rates of the rejection representationdescribed in Section 2.5, we analyze properties of truncated models assumingthe inverse Le´vy representation in this chapter.413.1. BNP ON TRANSFORMED CRMS WITH OBSERVED FEATURE ASSIGNMENTS3.1 BNP on transformed CRMs with observedfeature assignmentsBayesian nonparametric models are usually constructed using completelyrandom measures as priors. In this work, we consider a general family ofBNP models in which the prior is a discrete random measure ξ on somemeasurable space (S,S) that is a transformed CRM µ ∼ CRM(ν,H),ξ =∑kϑkδϕk ,where ϕk ∈ S are features that each can be expressed in the data, and ϑk ∈ R+are rates that determine the population frequency of the appearance of featureϕk.BNP models typically set ξ = µ (e.g. in feature allocation models) orξ = µ˜ (e.g. in mixture models), where recall that µ˜ is the normalizationof µ given in Eq. (2.4). In particular, we denote Ψ as the state space offeatures ψk, previous BNP models usually set S = Ψ, ϕk = ψk, and set ϑkto either θk or θk/µ(Ψ). In this chapter we consider more general cases of ξand BNP models based on it, for example, the edge-exchangeable networkmodels where ξ is the product measure ξ = µ × µ on S = Ψ2 (Cai et al.,2016; Williamson, 2016).In all of the models we consider, each observation n ∈ N has a variableXn indicating features assigned to nth observation, which itself is a discretemeasure on the same features as the underlying prior measure,Xn =∑kZnkδϕk ,where Znk ∈ N ∪ {0} is a nonnegative integer count corresponding to howmany “occurrences” of feature ϕk is expressed in observation n.There are two common likelihood processes, i.e., ways that the featureassignment Xn is generated from the prior measure ξ. In models whereeach observation is assigned to exactly one feature, we generate Xn from acategorical distribution. In particular, as long as ξ(S) <∞, we may generateXn via a categorical likelihood process Xn ∼ Categorical(ξ) defined byXn = δϕIn , In ∼ Categorical((ϑkξ(S))k∈N). (3.1)In models where each observation n can be assigned to multiple features,we generate Xn via an independent likelihood process Xn ∼ LP(G, ξ) defined423.2. SELF-PRODUCT CRM AND EDGE-EXCHANGEABLE NETWORK MODELSbyXn =∑kZn,kδϕk , Zn,kindep∼ G(·|ϑk), (3.2)where Zn,k = j denotes that there were j copies of feature ψk assigned tonth observation, and for all ϑ ∈ R+, G(·|ϑ) is a probability distributionon N ∪ {0}. These two are akin to the categorical likelihood process andindependent likelihood processes from standard CRM-based BNP models,but now the likelihood process is based on the transformed latent discretemeasure ξ.3.2 Self-product CRM and edge-exchangeablenetwork modelsIn this section, we consider the self-product of CRMs and an edge-exchangeablenetwork model built on it. We propose truncated simulation of the edge-exchangeable network in Section 3.3 and analyze the prior truncation errorlater in Section 3.4.Network data is increasingly available in many fields, such as socialnetworks and biological networks. Analyzing the distribution of connectionsand predicting future connections in these networks is of great interest. Amain tool for analyzing network data is the random graph, which treatsindividuals and connections in the network as vertices and edges of thegraph. Random graphs are traditionally described as exchangeable randomarrays with the property that the graph is vertex exchangeable, meaningthe distribution of its edges is invariant to any finite permutation of itsvertices. While vertex exchangeability is a natural requirement of the graph,it follows from the Aldous-Hoover theorem (Aldous, 1981; Hoover, 1979)that any vertex exchangeable graph is either empty or dense, implying thatthe number of edges grows quadratically with the number of vertices. Sincemost network data in real life is sparse, most Bayesian models of networkdata are inherently misspecified (Orbanz and Roy, 2014).To construct sparse graphs for network data, recent works have startedusing completely random measures to build priors for graph models andintroduced concepts of exchangeability other than vertex exchangeability. Apioneering work in this direction (Caron and Fox, 2017) models the graphbased on a discrete random measure on R2, where vertices of the graphare represented by locations of a CRM in R and the probability of an edge(edge rate) between a pair of vertices is given by the product of rates of433.2. SELF-PRODUCT CRM AND EDGE-EXCHANGEABLE NETWORK MODELSthe CRM at corresponding locations. Graphs constructed this way haveKallenberg exchangeability, which refers to the invariance of the distributionof the random measure to the orders of small increments in the domain. Thisgraph model achieves sparsity and exchangeability simultaneously and theproperties of this type of random graphs are thoroughly studied in recentworks (Veitch and Roy, 2015; Borgs et al., 2018; Caron and Fox, 2017). Basedon the CRM framework introduced by Caron and Fox (2017), more flexiblegraph models have been proposed that can model the clustering structure ofthe network data (Todeschini et al., 2020; Herlau et al., 2016; Naik et al.,2022).Other works (Cai et al., 2016; Williamson, 2016; Crane and Dempsey,2018; Campbell et al., 2018; Janson, 2018) model the distribution overgraph edges directly based on self-products of CRMs. In particular, theprior distribution of edge rates is determined by the self-product of a CRM,and edges are added to the graph sequentially by a likelihood process. Byconstruction, edges in the graph are i.i.d. samples based on edge rates, andhence the distribution of edges is invariant to the order in which they join thegraph. Graphs constructed this way are therefore called edge-exchangeablegraphs. One example of such graphs is based on the self-product of a Dirichletprocess, and given the normalized edge rates, one directed edge is generatedfrom the multinomial distribution and added to the graph in each step(Williamson, 2016). Another work (Cai et al., 2016) employs the self-productof discrete CRMs on R+, which is not normalized. Given the edge rates,multiple edges that are generated independently from a likelihood processare added to the graph in each step. By adjusting hyperparameters of theCRM, the edge-exchangeable graph can be either dense or sparse.3.2.1 Edge-exchangeable network modelWe use the edge-exchangeable network model as an example of BNP modelsbased on transformations of CRMs. In particular, we are interested inprobabilistic network models that utilize self-product CRMs as priors (Caiet al., 2016; Williamson, 2016). In these models, each label ψk in µ representsa vertex in the network, and each rate θk in µ is a “sociability” for vertex ψk.We generate (hyper)edges in the network with frequency that corresponds tothe product of the sociabilities for all vertices participating in the (hyper)edge.More precisely, ξ is a random measure on a d-dimensional space S = Ψd443.2. SELF-PRODUCT CRM AND EDGE-EXCHANGEABLE NETWORK MODELSthat has the formξ = µ(d) =∑i∈Nd̸=ϑiδϕi , ϑi =d∏j=1θij , ϕi = (ψi1 , ψi2 , ..., ψid), (3.3)here i is a d-dimensional multi-index, and Nd̸= is the set of such indices withall distinct components. Note that ξ = µ(d) is not itself a CRM, as it doesnot satisfy the independence condition.Each tuple ϕi = (ψi1 , . . . , ψid) corresponds to a potential (hyper)edge inthe network. We construct a growing sequence of networks by sequentiallygenerating edges from µ(d) in rounds n = 1, . . . , N . In each round n, edgescan be added to the graph using the independent likelihood process asdescribed in Eq. (3.2). Fig. 3.1 illustrates the generative process of an edge-exchangeable network based on the self-product of an completely randommeasure. We denote the probability of 0 under g to be g0(ϑ) := g(0|ϑ) forconvenience. By the Slivnyak-Mecke theorem (Last and Penrose, 2017), if gsatisfies ∫Rd+1− g0( d∏j=1θj) d∏j=1ν(dθj) <∞, (3.4)then finitely many edges are added to the graph in each round. Thisassumption is analogous to Eq. (2.5) in Section 2.3, and we will assume thiscondition holds in the following.Alternatively, under the finite mass assumption about the CRM:∫R+min(1, θ)ν(dθ) <∞,it holds that Ω := µ(d)(Ψd) <∞. In this case we may add only a single edgeper round n via a categorical likelihood process as outlined in Eq. (3.1). Thisconstruction has appeared in Williamson (2016), where µ follows a Dirichletprocess.Using either the independent likelihood process or the categorical likeli-hood process, the edges in the network after N rounds areN∑n=1Xn =∑i∈Nd̸=Ziδϕi , where Zi =N∑n=1Zn,i,453.2. SELF-PRODUCT CRM AND EDGE-EXCHANGEABLE NETWORK MODELSFigure 3.1: The generative process of edge-exchangeable network modelswhere Zi ∈ N ∪ {0} represents the count of edge i after N rounds. Fig. 3.1illustrates the generative process of edge-exchangeable network models.There are three points to note about this formulation. First, since theatom locations ϕi are not usually used, we can represent the network usingonly its array of edge countsEN = (Zi)i∈Nd̸= ∈ Nd, (3.5)where Nd denotes the set of integer arrays indexed by Nd̸= with finitely manynonzero entries. Note that Nd is a countable set. Second, by construction,the distribution of EN is invariant to reordering of the arrival of edges, andthus the network is edge-exchangeable (Cai et al., 2016; Crane and Dempsey,2018; Campbell et al., 2018; Janson, 2018). Finally, note that the networkEN , as formulated in Eq. (3.5), is generally a directed multigraph with noself-loops (due to the restriction to indices i ∈ Nd̸= rather than Nd).Similar to the species sampling model, a common goal in modelingnetwork data is to predict the edges that will appear in the future, whichcan be done by sampling from the posterior predictive distribution of ZN+1given Z1, . . . , ZN . Again, this requires us to marginalize out the latent edgefrequencies {ϑi}. However, this marginalization is not typically possible toobtain in closed-form; in fact, unlike the species sampling model, it is noteasy even when the underlying rate measure ν comes from an exponentialfamily. Additionally, one cannot simply first sample from the posteriordistribution of ϑ given Z1, . . . , ZN and then sample from ZN+1 given ϑ inpractice, because ϑ has infinitely many parameters.463.3. TRUNCATION OF EDGE-EXCHANGEABLE NETWORK MODELS(a) Truncation of µ(2) (b) Truncation of networkFigure 3.2: An illustration of truncation of self-product CRMs with d = 2,and the edge-exchangeable network built on it.3.3 Truncation of edge-exchangeable networkmodelsIn this section, we propose a truncated simulation algorithm for edge-exchangeable network models. We then analyze the prior truncation error,which is not only useful for generative simulation, but will also appear lateron as a key component of the truncated posterior inference algorithm andtruncation error analysis in subsequent chapters.To construct a truncated self-product measure, we first split the underly-ing CRM µ into a truncation µK and a tail part µK+µ = µK + µK+ =K∑k=1θkδψk +∞∑k=K+1θkδψkand construct the self-product ξK = µ(d)K from the truncation µK as inEq. (3.3). Fig. 3.2a provides an illustration of the truncation of µ(2), fromwhich it can be seen that that µ(2) can be decomposed into a sum of fourparts,µ(2) = (µK + µK+)(2) = µ(2)K +(µK × µK+ + µK+ × µK + µ2K+).Thus, while we only discard µK+ in truncating µ to µK , we discard three473.4. PRIOR TRUNCATION ERRORparts in truncating µ(2)K to µ(2)K ; and in general, we discard 2d − 1 parts ofµ(d) when truncating it to µ(d)K .Once the measure is truncated, a truncated network—based on µ(d)K —can be constructed in the same manner as the original network using theindependent likelihood process Eq. (3.1) or categorical likelihood processEq. (3.2), as illustrated in Fig. 3.2b.Note that if the rejection representation is used in simulating the trun-cated CRM µ, it is possible that µK = 0 and so is µ(2)K . This causes anadditional challenge for the normalization of the truncated networks inEq. (3.1). Fortunately, the inverse Le´vy representation provides an avenue toaddress this issue. The rates θk are all guaranteed to be nonzero—meaning aslong as K ≥ d, the categorical likelihood process is well-defined—and are de-creasing, which enables our theoretical analysis. Additionally, Theorem 2.5.1guarantees that we can always assume that the inverse Le´vy representationis accessible, even though we are actually using the rejection representationin practice.As in the full network model, we do not typically use the (hyper)edgefeatures, and hence we can represent the truncated network at round Nusing only its array of edge countsEN,K = (Zi,K)i∈Nd̸= ∈ Nd,where Nd denotes the set of integer arrays indexed by Nd̸= with finitelymany nonzero entries, and Zi,K = 0 for any index i ∈ Nd̸= such that somecomponent ij > K. We keep EN and EN,K in the same space in order tocompare their distributions in analyzing the truncation error.3.4 Prior truncation errorIn this section, we use the truncated network model as an example andanalyze the truncation error of truncated BNP models with observed featureassignments. However, the method can be easily extended to other modelswith observed feature assignments, such as the species sampling models(Fisher et al., 1943).We formulate the approximation error incurred by truncation in forwardsimulation as the total variation distance between the marginal networkdistributions, i.e., of EN and EN,K . Similar marginal truncation errorbounds for various CRM models have been developed in previous work(Campbell et al., 2019). However, existing results on CRMs do not apply483.4. PRIOR TRUNCATION ERRORto the network model here, which is constructed on the self-product CRMand leads to different truncation structure as discussed in Section 3.3. Inthis section, we provide truncation error bounds for a general d-dimensionalnetwork model.The first step in the analysis of truncation error is to show that the totalvariation distance between the marginal distributions is bounded above bythe probability that there are edges in the full network EN involving verticesbeyond the truncation level K. To this end, we denote the maximum vertexindex of EN to beIN := maxi∈Nd̸=(maxj∈[d]ij)s.t. Zi > 0,and note that by definition, IN ≤ K if and only if all edges in EN fall in thetruncated region. This result is illustrated in Lemma 3.4.1.Lemma 3.4.1. Let µ =∑∞k=1 θkδψk ∼ CRM(ν,H) be a completely randommeasure, and µK =∑Kk=1 θkδψk be its truncation to K atoms. Let µ(d) bethe self-product of µ, and µ(d)K be the self-product of µK . Let PN and PN,Kbe the marginal distributions of edge sets EN , EN,K ∈ Nd under either theindependent or categorical likelihood process. ThenDTV (PN , PN,K) ≤ 1− P (IN ≤ K) .As mentioned above, the network EN is in general a directed multigraphwith no self loops. However, Lemma 3.4.1—and the downstream truncationerror bounds presented —also apply to any graph E′N = (Z′i)i∈Nd̸= that is afunction of the original graph E′N = f(EN ) such that truncation commuteswith the function, i.e., E′N,K = f(EN,K). For example, to obtain a truncationerror bound for the common setting of undirected binary graphs, we generatethe directed multigraph EN as above and construct the undirected binarygraph E′N viaZ ′i = 1Zi>0 · 1i1<i2<···<id , i ∈ Nd̸=. (3.6)Corollary 3.4.2 provides the precise statement of the result; note that thebound is identical to that from Lemma 3.4.1.Corollary 3.4.2. Let E′N = (Z′i)i∈Nd̸= ∈ Nd be the set of edges of a networkwith truncation E′N,K ∈ Nd, and denote their marginal distributions P ′N andP ′N,K . If there exists a measurable function f such thatE′N = f(EN ) and E′N,K = f(EN,K),493.4. PRIOR TRUNCATION ERRORthenDTV(P ′N , P′N,K) ≤ 1− P (IN ≤ K) .Edges added via independent likelihood processWe now specialize Lemma 3.4.1 to the setting where µ is a CRM generatedby a series representation of the form Eq. (2.8), and the network is generatedvia the independent likelihood process from Eq. (3.2). As a first step towardsa bound on the truncation error for general hypergraphs with d > 1 inTheorem 3.4.4, we present a simpler corollary in the case where d = 2,which is of direct interest in analyzing the truncation error of popular edge-exchangeable networks.Corollary 3.4.3. In the setting of Lemma 3.4.1, suppose µ is a CRMgenerated from the series representation Eq. (2.8), edges are generated fromthe independent likelihood process Eq. (3.2), and d = 2 ≤ K. ThenP (IN ≤ K) ≥ exp (−N ·BK) ,whereBK = BK,1 +BK,2 +BK,3BK,1 = E[∫R2+− log g0 (τ(U1, γ1 + ΓK)τ(U2, γ2 + ΓK)) dγ1dγ2]BK,2 = 2E[∫ ΓK0(K − 1)ΓK∫ ∞ΓK− log g0 (τ(U1, γ1)τ(U2, γ2)) dγ1dγ2]BK,3 = 2E[∫R+− log g0 (τ(U1,ΓK)τ(U2, γ + ΓK)) dγ].The proof of this result (and Theorem 3.4.4 below) in the appendixfollows by representing the tail of the CRM as a homogeneous (unit-rate)Poisson process on [ΓK ,∞). Though perhaps complicated at first glance,an intuitive interpretation of the truncation error terms BK,i is providedby Fig. 3.2b. The term BK,1 corresponds to the truncation error arisingfrom the upper right quadrant, where both vertices participating in the edgewere in the discarded tail region. The term BK,2 is the truncation errorarising from the bottom right and upper left quadrants, where one of thetwo vertices participating in the edge was in the truncation, and the otherwas in the tail. Finally, BK,3 represents the truncation error arising from503.4. PRIOR TRUNCATION ERRORedges in which one vertex was at the boundary of tail and truncation, andthe other was in the tail.Theorem 3.4.4 is the generalization of Corollary 3.4.3 from d = 2 to thegeneral setting of arbitrary hypergraphs with d > 1. The bound is analogousto that in Corollary 3.4.3—with BK expressed as a sum of terms, eachcorresponding to whether vertices were in the tail, boundary, or truncationregion—except that there are d > 1 vertices participating in each edge,resulting in more terms in the sum. Note that Theorem 3.4.4 also guaranteesthat the bound is not vacuous, and indeed converges to 0 as the truncationlevel K →∞ as expected.Theorem 3.4.4. In the setting of Lemma 3.4.1, suppose µ is a CRMgenerated from the series representation Eq. (2.8), edges are generated fromthe independent likelihood process Eq. (3.2), and 1 < d ≤ K. ThenP (IN ≤ K) ≥ exp (−N ·BK) ,whereBK = BK,1 +BK,2 +BK,3, (3.7)BK,1 = E[∫[ΓK ,∞)d− log g0(θ˜)dγ]BK,2 =d−1∑ℓ=1(dℓ)E[(K − 1) !(K − 1− ℓ) !Γ−ℓK∫[0,ΓK ]ℓ×[ΓK ,∞)d−ℓ− log g0(θ˜)dγ]BK,3 =d−1∑ℓ=1ℓ(dℓ)E[(K − 1) !(K − ℓ) ! Γ−(ℓ−1)K∫[0,ΓK ]ℓ×[ΓK ,∞)d−ℓ−δγℓ=ΓK log g0(θ˜)dγ],where δ· is the Dirac delta, dγ =∏dj=1 dγj, and θ˜ =∏dj=1 τ(Uj , γj). Fur-thermore, limK→∞BK = 0.The same geometric intuition from the d = 2-dimensional case appliesto the general hypergraph truncation error in Eq. (3.7). The term BK,1corresponds to the error arising from the edges whose vertices all belongto the tail region µK+. Each term in the summation in BK,2 correspondsto the error arising from edges that have ℓ out of d vertices belonging tothe truncation µK . Each term in the summation in BK,3 corresponds tothe error arising from the edges that have ℓ− 1 out of d vertices belongingto the truncation µK and have one vertex exactly on the boundary of thetruncation. Note that we obtain Corollary 3.4.3 by taking d = 2 in Eq. (3.7).513.5. EXAMPLESEdges added via categorical likelihood processWe may also specialize Lemma 3.4.1 to the setting where the networkis generated via the single-edge-per-step categorical likelihood process inEq. (3.1). However, truncation with the categorical likelihood process posesa few key challenges. From a practical angle, certain choices of seriesrepresentation for generating µ may be problematic. For instance, whenusing the rejection representation Eq. (2.12) in the typical case where νˆ ̸= ν,there is a nonzero probability thatK∑k=11 [θk > 0] < d,meaning there aren’t enough accepted vertices in the truncation to generatea single edge. In this case, the categorical likelihood process—which mustgenerate exactly one edge per step—is ill-defined. We therefore assumein Theorem 3.4.5 that the truncated CRM is generated via inverse Le´vyrepresentation.Theorem 3.4.5. In the setting of Lemma 3.4.1, suppose µ is a CRMgenerated from the inverse Le´vy representation Eq. (2.9), edges are generatedfrom the categorical likelihood process Eq. (3.1), and 1 < d ≤ K. ThenP (IN ≤ K) ≥ (1−BK)Nd ≥ 0,whereBK = E[∫ ∞−∞Q(ΓK , x)(∫ 10Q(ΓKu, x)du)K−d(ddxe∫∞0 Q(ΓK+γ, x)−1 dγ)dx],andQ(u, t) = exp(−ν←(u)e−t) and ΓK ∼ Gam(K, 1).Furthermore, limK→∞BK = 0.3.5 ExamplesWe now apply the results of Section 3.4 to binary undirected networks (d = 2)constructed via Eq. (3.6) from common edge-exchangeable networks (Caiet al., 2016). We derive the convergence rate of truncation error, and providesimulations of the expected number of edges and vertices under the infiniteand truncated network. In each simulation we use the rejection representationto simulate µ, and run N = 10, 000 steps of network construction.523.5. EXAMPLES(a) Dense graph with α = 0 (b) Sparse graph with α = 0.6Figure 3.3: Beta-independent Bernoulli network3.5.1 Beta-Bernoulli process networkDense networkWhen α = 0, the binary beta-Bernoulli graph is dense. We use the rejectionrepresentation specified in Eq. (2.13), where the inverse of the Le´vy measureof the proposal isνˆ←B (u) = e−u/γ′ ,dνdνˆB= (1− θ)λ−1.Therefore the rejection representation Eq. (2.12) of BP(γ, λ, 0) can be writtenasθk = Tk1(Uk ≤ (1− Tk)λ−1), Tk = e−Γk/γ′ .In the appendix, we show that there exists K0 ∈ N such that∀K ≥ K0, BK ≤ 12γ(1− e−1)λ−2(γ′1 + γ′)K.This implies that the truncation error of the dense binary beta-Bernoullinetwork converges to 0 geometrically in K. Fig. 3.3a corroborates this resultin simulation with λ = 2 and γ = 1; it can be seen that for the dense beta-Bernoulli graph, truncated graphs with relatively low truncation level—inthis case, K ≈ 50—approximate the true network model well.533.5. EXAMPLESSparse networkWhen α ∈ (0, 1), the binary beta-Bernoulli graph is sparse. We use therejection representation specified in Eq. (2.14), and the inverse of the Le´vymeasure of the proposal distribution isνˆ←B (u) =(1 +αuγ′)−1/α,dνdνˆB= (1− θ)λ+α−1.Therefore the rejection representation Eq. (2.12) of BP(γ, λ, α) can be writtenasθk = Tk1(Uk ≤ (1− Tk)λ+α−1), Tk = (1 + αΓk/γ′)−1/α.In the appendix, we show that there exists K0 ∈ N such that∀K ≥ K0, BK ≤ 6α(γ′α−1)1/α eγ′α−1 (K − 1)α−1α .This bound suggests that the truncation error for the sparse binary beta-Bernoulli network converges to 0 much more slowly than for the dense graph.Fig. 3.3b corroborates this result in simulation with λ = 2, γ = 1, andα = 0.6; it can be seen that for the sparse beta-Bernoulli graph, truncatedgraphs behave significantly differently from the true graph for moderatetruncation levels. In practice, one should select an appropriate (large) valueof K using our error bounds as guidance, and use sparse data structures toavoid undue memory burden.3.5.2 Gamma-independent Poisson networkDense networkWhen α = 0, the gamma-Poisson graph is dense, and we choose the proposalmeasure νˆG to be the one in 2.15 such thatνˆ←G (u) = 1/(λ(e(γλ)−1u − 1)),dνdνˆG= (1 + λθ)e−λθ.Therefore, the rejection representation in Eq. (2.12) has the formθk = Tk1(Uk ≤ (1 + λTk)e−λTk), Tk =1λ(e(γλ)−1Γk − 1) .In the appendix, we show that there exists K0 ∈ N such that∀K ≥ K0, BK ≤ 6γλ(γλ1 + γλ)K−1.543.5. EXAMPLES(a) Dense graph with α = 0 (b) Sparse graph with α = 0.6Figure 3.4: Gamma-independent Poisson graphAgain, for the dense network BK converges to 0 geometrically, indicating thattruncation may provide reasonable approximations to the original network.Fig. 3.4a corroborates this result in simulation with λ = 2 and γ = 1; forK ≈ 50, no vertices are discarded on average by truncation.Sparse networkWhen α ∈ (0, 1), the gamma-Poisson graph is sparse, and we choose theproposal measure νˆG to be the one specified in Eq. (2.16) such thatνˆ←G (u) = (γ′u−1)1/d, γ′ = γλ1−ααΓ(1− α) ,dνdνˆG= e−λθ.Therefore the rejection representation of the Gamma process in Eq. (2.12)has the formθk = Tk1(Uk ≤ e−λTk), Tk =(γ′Γ−1k)1/α.In the appendix, we show there exists K0 ∈ N such that∀K ≥ K0, BK ≤ 12γ2λ1−α(1− α)Γ(1− α)(3γ′K − 1) 1−αα.Again, for the sparse network BK converges to 0 slowly, suggesting thatthe truncation error for the sparse binary gamma-Poisson graph convergesmore slowly than for the dense graph. Fig. 3.4b corroborates this resultin simulation with λ = 2, γ = 1, and α = 0.6; for a moderate range of553.5. EXAMPLEStruncation values K ≤ 100, the truncated graph behaves very differentlyfrom the true graph. Therefore in practice, one should follow the aboveguidance for the beta-Bernoulli network: select a value of K using our errorbounds, and avoid intractable memory requirements by using sparse datastructures.3.5.3 Normalized-gamma process networkIn this example we generate the network via the categorical likelihood processand use the self-product of the normalized gamma process as the prior of edgerates. The rates of edges are simulated using rejection representation samewith that used in Section 3.5.2 and truncation at different levels are obtainedaccordingly. In each round only one edge sampled from the categoricaldistribution is added to the network.When the discount parameter α = 0, the normalized gamma processis the Dirichlet process. When α > 0, the normalized gamma process hassimilar behavior to the Pitman-Yor process, (Lijoi et al., 2007b) which couldnot be represented as a BNP model based on CRM.By comparing the pattern in 3.5 to the pattern in Fig. 3.4, we seethat the networks generated via categorical likelihood process demonstratesimilar structures to networks generated using independently likelihoodprocess, conditioned on the same CRM. Specifically, Fig. 3.5a indicatesthat the network constructed using the Dirichlet process is dense, implyingthat truncated models even with small truncation level provide descentapproximations to the true network. In contrast, Fig. 3.5b shows thatthat the truncated networks approximate the true network poorly when thenetwork is sparse, and this problem is not solved by increasing the truncationlevel. This simulation indicates that truncated models does not provide goodapproximation to BNP models showing power law behavior.563.5. EXAMPLES(a) Dense graph with α = 0 (b) Sparse graph with α = 0.6Figure 3.5: Normalized-gamma process graph57Chapter 4Truncated inference for BNPmodels with observed labelsPosterior inference for Bayesian nonparametric models is not trivial as BNPmodels have an infinite number of parameters. However, statistical inferenceis concerned only with parameters related to the underlying patterns andstructures on the data set to be fitted. Since only a finite subset of parametersof a BNP models are invoked by any given dataset, it is both necessary anddesirable to avoid the infinite parameters of the BNP model that are inactive.Collapsed inference methods marginalize out all inactive parameters andperform inference over only the invoked features. However, such methodsare restricted to conjugate BNP models and suffer from high computationalcost as feature assignments of different observations no longer hold the i.i.d.structure. Conditional inference via slice sampling adjusts the subset ofactive parameters in each iteration during sampling. However, numericalintegration is required for each feature added when expanding the truncationlevel, leading to high computational cost.Approximate inference methods based on truncation offer a generally ap-plicable and computationally efficient alternative to exact posterior inference.Since posterior inference is performed on the truncated model, it is onlyuseful when an appropriate truncation level can be selected such that theapproximation error between the truncated model and the true posterior canbe controlled. Previous work on truncated inference chooses the truncationlevel either based on the prior approximation error or based on statistics ofthe truncated model, which do not necessarily reflect the closeness betweenthe truncated posterior and the true posterior.In this chapter, we continue the study of BNP models with observedfeature assignments and turn to analyze the posterior approximation error.We make use of the fact that when feature assignments are observed, thejoint likelihood of the truncated model is an upper bound of that of thetrue BNP model. We use this fact to derive an upper bound of the totalvariation distance between the truncated posterior and the true posterior,which decreases with K as expected and increases with N . We then propose584.1. POSTERIOR TRUNCATION ERRORan adaptive truncated sampling algorithm that increases the truncation leveluntil the desired error level is achieve.4.1 Posterior truncation errorFor BNP models with observed feature assignments, the target of posteriorinference is the weights of atoms of the CRM µ, which represent popula-tion frequency of the corresponding clusters. For BNP models based ontransformations of CRMs, the same transformation of the posterior CRMsleads to the posterior of the quantifies of interest. For example, we build theedge-exchangeable network model on the self-product of the CRM µ, andin this case the posterior of µ determines the rates of edges in the networkmodel.We perform posterior inference via truncation, and to measure the approx-imation error, we derive computable bounds on the total variation distancebetween the truncated and full BNP posterior. We illustrate the developmentand application of our method for inference in edge-exchangeable networkmodels, where exact marginalization is not possible.In this chapter we continue using the edge-exchangeable network modelas an example in deriving the truncation error of BNP models based ontransformations of CRMs. We begin by examining the density of the posteriordistribution of the Kth rate θK from the inverse Le´vy representation for somefixed K ∈ N, the unordered collection of rates (θk)K−1k=1 such that θk ≥ θK ,and the parameters ϱ ∈ Rm of the CRM rate measure νϱ given the observedset of edges EN . We denote νϱ(θ) to be the density of νϱ(dθ) and p(ϱ) to bethe prior density of ϱ, both with respect to the Lebesgue measure. Giventhese definitions the posterior density can be expressed using Bayes’ rule asp(ϱ, θ1:K | z1:N ) ∝p(ϱ) · e−νϱ[θK ,∞) ·K∏k=11[θK ≤ θk]νϱ(θk)· N∏n=1K∏i∈[K]d̸=g(zni |d∏j=1θij ) · p(zK+1:N |θ1:K , ϱ), (4.1)where [K]d̸= is the subset of Nd̸= such that maxj∈[d] ij ≤ K, and zK+1:N isshorthand for the set of (zni)n∈[N ],i∈Nd̸= such that i ∈ [K]d̸=.Note that all the factors on the first row of Eq. (4.1) correspond to theprior of ϱ and (θk)Kk=1, and the first factor on the second row corresponds to594.1. POSTERIOR TRUNCATION ERRORthe likelihood of the portion of the network involving only vertices 1, · · · ,K;these are straightforward to evaluate, though νϱ[θK ,∞) will occasionallyrequire a 1-dimensional numerical integration. The last factor correspondsto the likelihood of labels involving related to the tail part of the CRM, andis not generally possible to evaluate exactly.To handle this term, suppose that K is large enough that zK+1:N = 0. Thenp(zK+1:N |θ1:K , ϱ) = P (IN ≤ K | Γ1:K , ϱ), i.e., the probability that no portionof the network involves vertices of index larger than K. Theorem 4.1.1provides upper and lower bounds on this probability. This theorem alsomakes the dependence of the bound on the rate measure parameters ϱnotationally explicit via parameterized series representation components τϱand gϱ from Eq. (2.8). Finally, though Theorem 4.1.1 applies to general seriesrepresentations, we require only the specific instantiation for the inverse Le´vyrepresentation in the present context.Theorem 4.1.1. The conditional probability P (IN ≤ K | U1:K ,Γ1:K , ϱ) sat-isfies1 ≥ P (IN ≤ K | U1:K ,Γ1:K , ϱ) ≥ exp (−N ·B(U1:K ,Γ1:K , ϱ)) ,whereB(U1:K ,Γ1:K , ϱ) =d−1∑ℓ=0(dℓ) ∑L⊆[K]|L|=ℓB(U1:K ,Γ1:K , ϱ,L)B(U1:K ,Γ1:K , ϱ,L) =∫[ΓK ,∞)d−ℓE− log g0∏j∈Lθjd−ℓ∏j=1τϱ(U′j , γj) | θ1:Kdγ,where dγ =∏d−ℓj=1 dγj and U′ji.i.d.∼ gϱ.Since U1:K is unused in the inverse Le´vy representation, in the presentcontext we use the notation B(Γ1:K , ϱ) for brevity. The bound in Theo-rem 4.1.1 implies that as long as K is set large enough such that bothZK+1:N = 0 and B(Γ1:K , ϱ) ≤ ϵ/N for some ϵ > 0 then1− ϵ ≤ p(zK+1:N |θ1:K , ϱ) ≤ 1.Therefore as K increases, this term should become ≈ 1 and have a negligibleeffect on the posterior density. We use this intuition to propose a truncatedMetropolis–Hastings algorithm that sets K > IN , ignores the p(zK+1:N |θ1:K , ϱ)term in the acceptance ratio, and fixes ZK+1:N to 0.604.1. POSTERIOR TRUNCATION ERRORDenote the truncated posterior distribution of ϱ, θ1:K to be Πˆµ, withdensity proportional toπˆµ(ϱ, θ1:K) =p(ϱ) · e−νϱ[θK ,∞) ·K∏k=11[θK ≤ θk]νϱ(θk)·N∏n=1K∏i∈[K]d̸=g(zni |d∏j=1θij ),and the exact posterior distribution of the truncated part (ϱ, θ1:K) to be Πµwith density proportional toπµ(ϱ, θ1:K) = πˆµ(ϱ, θ1:K) · p(zK+1:N |θ1:K , ϱ).Theorem 4.1.2 provides a rigorous analysis of the error involved in using thetruncated sampler.Theorem 4.1.2. Fix K > IN . If for some η ∈ [0, 1),Πˆµ {B(Γ1:K , ϱ) ≤ ϵ/N} ≥ 1− η,thenDTV(Πˆµ,Πµ)≤ 3(ϵ+ η)2− ϵη. (4.2)Crucially, as long as one can obtain samples from the truncated posteriordistribution, one can estimate the bound in Eq. (4.2) using sample estimatesof the tail probability Πˆµ {B(Γ1:K , ϱ) ≤ ϵ/N} ≥ 1 − η. Therefore, one cancompute the bound Eq. (4.2) without needing to evaluate p(zK+1:N |θ1:K , ϱ)exactly.Fig. 4.1 illustrates the idea behind Theorem 4.1.2. The blue curverepresents πµ (proportional to the true posterior density), while the orangecurve represents πˆµ (proportional to the approximate posterior density). Wesplit the posterior parameter space into two pieces: a central region whereB(Γ1:K , ϱ) ≤ ϵ/N , and the tail region corresponding to the complement. Inthe central region, we know that 1− ϵ ≤ p(zK+1:N |θ1:K , ϱ) ≤ 1, and so πµ ≈ πˆµwith error uniformly bounded by ϵ. In the tail region, we have that πµ ≤ πˆµ;therefore the mass of πˆµ (roughly η) dominates the mass of πµ in this secondregion. Hence we can use samples from Πˆµ to bound the tail mass of Πµ. Bycombining these bounds on the tail mass and the error within the centralregion, we obtain a bound on the total variation between DTV(Πˆ,Π).614.1. POSTERIOR TRUNCATION ERRORFigure 4.1: An illustration of the posterior truncation errorProof of Theorem 4.1.2Here we derive the total variation distance between Πˆµ and Πµ in moredetail. Given any fixed K ∈ N and ϵ > 0, we define the subset of state spaceA ={ZK+1:N = 0, B(Γ1:K , ϱ) ≤ ϵN+1}. By assumption,Πˆµ (A) ≥ 1− η.Further, by applying the bound from Theorem 4.1.1, we know that for statesin A,1− ϵ ≤ p(ZK+1:N = 0 | θ1:K , ϱ) ≤ 1.Suppose πµ is the right hand side (RHS) of 4.1 that is proportional tothe density of Πµ, and πˆµ is the RHS of Eq. (4.1) removing the termp(zK+1:N |θ1:K , ϱ), which is proportional to the density of Πˆµ. Then we haveπµ ≤ πˆµ, and within A, πµ ≥ (1− ϵ)πˆµ.We define the normalization constants Cµ, Cˆµ such that∫π/Cµ = 1 and∫πˆµ/Cˆµ = 1. Then the above bounds yield the upper and the lower boundfor Cµ:Cµ =∫πµ ≤∫πˆµ = CˆµCµ =∫Aπµ +∫Acπµ ≥∫A(1− ϵ)πˆµ = (1− ϵ)∫A πˆµ∫πˆµ∫πˆµ ≥ (1− ϵ)(1− η)Cˆµ.624.2. ADAPTIVE TRUNCATED METROPOLIS–HASTINGSTherefore, (1− ϵ)(1− η)Cˆµ ≤ Cµ ≤ Cˆµ, and hence the subset A also haslarge probability under the true posterior distribution:Πµ(A) =∫Aπµ/Cµ ≥ (1− ϵ)∫Aπˆµ/Cµ ≥ (1− ϵ)∫Aπˆµ/Cˆµ ≥ (1− ϵ)(1− η).The above results yield the total variation bound via12∫ ∣∣∣∣∣ πµCµ − πˆµCˆµ∣∣∣∣∣=12∫A∣∣∣∣∣ πµCµ − πˆµCˆµ∣∣∣∣∣+ 12∫Ac∣∣∣∣∣ πµCµ − πˆµCˆµ∣∣∣∣∣≤12∫A∣∣∣∣∣ πµCµ − πµCˆµ∣∣∣∣∣+ 12∫A∣∣∣∣∣ πµCˆµ − πˆµCˆµ∣∣∣∣∣+ 12 (Πµ(Ac) + Πˆµ(Ac))=12∣∣∣∣∣ 1Cµ − 1Cˆµ∣∣∣∣∣∫Aπµ +121Cˆµ∫A(πˆµ − πµ) + 12(Πµ(Ac) + Πˆµ(Ac))≤12(1− (1− ϵ)(1− η)) + 12ϵ+12(1− (1− η)(1− ϵ) + η)=3(ϵ+ η)2− ϵη.4.2 Adaptive truncated Metropolis–HastingsTheorem 4.1.2 suggests the following adaptive truncation procedure:1. Pick a value of K > IN and desired total variation error ξ ∈ (0, 1).2. Obtain samples from the truncated posterior.3. Minimize the bound in Eq. (4.2) over ϵ ∈ (0, 1), using the samples toestimate η = 1− Πˆµ {B(Γ1:K , ϱ) ≤ ϵ/N} for each value of ϵ.4. If the minimum bound exceeds ξ, increase K and return to step 2.5. Otherwise, return the samples.In practice, we can start by initializing K to IN + 1. In order to decidehow much to increase K by in each iteration, we can use the sequentialrepresentation to extrapolate the total variation bound Eq. (4.2) to largervalues of K without actually performing MCMC sampling. In particular, for634.3. POSTERIOR INFERENCE EXPERIMENTSeach posterior sample, we use its hyperparameters to generate additional ratesfrom the sequential representation. We then use these extended samples tocompute the total variation error guarantee as per step 3 above. We continueto generate additional rates (typically doubling the number each time) untilthe predicted total variation guarantee is below the desired threshold. Finally,we use linear interpolation between the last two predicted errors to find thenext truncation level K that matches the desired (log) error threshold. Thisfixes a new value of K; at this point, we return to step 2 above and iterate.4.3 Posterior inference experimentsIn this section, we examine the properties of the proposed adaptive truncatedinference algorithm for the beta-independent Bernoulli network model withdiscount α ∈ (0, 1), concentration λ > 1, mass γ > 0, unordered collectionof rates θ1:K−1, and Kth rate from the sequential representation θK . Inorder to simplify inference, we transform each of these parameters to anunconstrained version:α =exp(αu)1 + exp(αu), λ = 1 + exp(λu), γ = exp(γu),θK =exp(θu,K)1 + exp(θu,K), θk =θK + exp(θu,k)1 + exp(θu,k)k = 1, . . . ,K − 1.We use a Markov chain Monte Carlo algorithm that includes an exactGibbs sampling move for γ, and separate Gaussian random-walk Metropolis–Hastings moves for αu, λu, θu,K , all θu,k such that vertex k has degree 0(jointly), and each θu,k such that vertex k has nonzero degree (individually).4.3.1 Synthetic dataWe first apply the model to synthetic data simulated from the generativemodel. We simulate a sparse network with parameters λ = 2, γ = 1, α = 0.2,and N = 105, and a dense network with λ = 2, γ = 1, α = 0, and N = 107.In both settings we set the truncation level for data generation to 500, thedesired total variation bound from Eq. (4.2) to 0.01, and initialize the samplerwith α = 0.4, λ = 5, γ = 2 and θ generated from the sequential representation.All Metropolis–Hastings moves have proposal standard deviation 0.1 exceptthe sparse network αu move, which has standard deviation 0.03.Figs. 4.2 and 4.3 show histograms of 5,000 marginal posterior samples ofthe hyperparameters for the dense (true α = 0) and sparse (true α = 0.2)644.3. POSTERIOR INFERENCE EXPERIMENTS(a) True α = 0 (b) True γ = 1 (c) True λ = 2Figure 4.2: Posterior of α, γ, λ for the dense simulated network. Orangehistograms depict the first adaptation iteration, and blue histograms depictthe second and final iteration.(a) True α = 0.2 (b) True γ = 1 (c) True λ = 2Figure 4.3: Posterior of α, γ, λ for the sparse simulated network. Orangehistograms depict the first adaptation iteration, and blue histograms depictthe second and final iteration.networks. In both cases, the approximate posterior in the first round of adap-tation (orange histogram) does not concentrate on the true hyperparametervalues, despite the relatively large number of generative rounds N .Fig. 4.4—which displays the truncation error and predictive adaptationprocedure—shows why this is the case. In both networks, the first adaptationiteration identifies a large truncation error. After a single round of adapta-tion, the approximate posterior distributions (blue histograms) in Figs. 4.2and 4.3 concentrate more on the true values as expected, and the truncationerrors fall well below the desired threshold (log10(0.01) = −2). It is worthnoting that the predictive extrapolation appears to be quite conservative654.3. POSTERIOR INFERENCE EXPERIMENTS(a) Dense synthetic network (b) Sparse synthetic networkFigure 4.4: Visualization of the truncation expansion procedure. Black circlesdenote the truncation level and error in each iteration of the adaptation.The adaptation stops when the truncation error falls below the desiredthreshold (here log10(0.01) = −2). Grey dashed lines and circles visualizethe predicted truncation error using rates generated from the sequentialrepresentation. The vertical dotted line shows that the algorithm selects avalue of K that attempts to match the desired log10 error threshold of −2using the predictions.in these examples, and especially in the dense network. This happens be-cause the approximate posterior for the dense network (respectively, sparsenetwork) assigns mass to higher values of α (respectively, γ) than it should,which results in larger truncation error and thus a larger predicted requiredtruncation level.4.3.2 Real network dataNext, we apply the model to a Facebook-like Social Network6 (Panzarasaet al., 2009). The original source network contains a sequence of 61,734weighted, time-stamped edges, and 1,899 vertices. We preprocess the databy removing the edge weights, binning the edge sequence into 30-minuteintervals, and removing the initial transient of network growth, resulting in1,899 vertices and 10,435 edges over N = 6,427 rounds of generation. We6Available at https://toreopsahl.com/datasets/664.3. POSTERIOR INFERENCE EXPERIMENTSFigure 4.5: Posterior of α, γ, λ in the Facebook-like networkagain set a desired total variation error guarantee of 0.01 during inference.All Metropolis–Hastings moves have proposal standard deviation 0.1 exceptthe αu and degree-0 θu moves, which have standard deviation 0.04 and 0.03,respectively. We initialize the sampler to α = 0.1, γ = 2, λ = 20 and samplerates θ from the prior sequential representation.Fig. 4.5 shows the posterior marginal histograms for the hyperparametersα, λ, γ in both the first iteration (orange) and the second and final iteration(blue) of truncation adaptation. The posterior distribution suggests thatthe network is dense (i.e. α ≈ 0). This conclusion is supported both bythe close match of 100 samples from the posterior predictive distribution,shown in Figs. 4.6a and 4.6b, and the findings of past work using this dataCaron and Fox (2017). Further, as in the synthetic examples, the truncationadaptation terminates after two iterations; but in this case the histogramsdo not change very much between the two. This is essentially because thetruncation error in the first iteration is relatively low (≈ 0.02), leading to areasonably accurate truncated posterior and hence accurate predictions ofthe truncation error at higher truncation levels.Our posterior error bounds can be estimated using only samples fromthe truncated posterior and a generative simulation error bound; both aretypically straightforward to obtain. Thus our method is adaptive: given aninitial fixed truncation level, the practitioner can sample from the truncatedposterior, estimate the total variation error, and if it is deemed to high thepractitioner can increase the truncation level.674.3. POSTERIOR INFERENCE EXPERIMENTS(a) # Vertices vs. # edges (b) Degree density (c) Truncation expansionFigure 4.6: (4.6a,4.6b): The characteristics of the observed Facebook-likenetwork (black) and 100 samples from the posterior predictive distribution(blue). (4.6c): The truncation adaptation process, with truncation level anderror bound in each iteration (black circles) and predictive truncation errors(gray dashed lines).68Chapter 5Truncated inference for BNPmodels with latent labelsMany Bayesian nonparametric models have latent features and feature as-signments. Examples of such models include distribution approximation,clustering analysis using mixture models, and latent feature analysis. In thischapter we turn to the truncated inference for BNP models based on CRMswith latent feature assignments. We analyze the posterior approximationerror of the truncated model and propose a truncated inference algorithmthat has error guarantees. We extend the theoretical result in Chapter 4 andanalyze the approximation error using upper and lower bounds on the ratioπµ/πˆµ, where πµ and πˆµ are the unnormalized density of the target distribu-tion and the approximate distribution, respectively. Additionally, we find anupper and a lower bound needed to quantify the posterior approximationerror of the truncated inference algorithm. We also show that the posteriorapproximation error derived earlier in Theorem 4.1.2 can be recovered as aspecial case of the result in this chapter.5.1 Truncated BNP models with latent featureassignmentWe turn to Bayesian nonparametric models based on CRMs with latentfeature assignments. Examples of BNP models in this family include latentfeature allocation models (Griffiths and Ghahramani, 2011; Teh and Gorur,2009) and latent trait allocation models (Titsias, 2007; Zhou et al., 2012;695.1. TRUNCATED BNP MODELS WITH LATENT FEATURE ASSIGNMENTBroderick et al., 2014). The generative definition of such models is as followsµ =∞∑k=1θkδψk ∼ CRM(ν,H),Xn | µ =∞∑k=1Zn,kδψki.i.d.∼ LP(G,µ), n = 1, . . . N (5.1)Yn | Xn indep∼ F( ∞∑k=1zn,kψk). (5.2)Compared to the BNP models studied in Chapter 4, the features ψk andfeature assignments Zn,k are latent in the model specified in Eq. (5.1).In addition, we assume that the distribution F of each observation Ynconditioned on its latent features is parameterized by a linear combinationof its latent features, as specified in Eq. (5.2). Depending on the choice ofthe feature assignment distribution G, Zn,k could be in {0, 1} or take anypositive values in N. For example, G could be the Bernoulli distribution,which represents whether a feature is active or not in one observation, as inthe Indian buffet process; G can also be the Poisson distribution (Titsias,2007) or the negative binomial distribution (Broderick et al., 2014; Zhouet al., 2012), which allows one feature to be expressed in different magnitudesin one observation. In this chapter we restrict our discussion to BNP modelsbased on completely random measures where each observation can expressmultiple features. We leave BNP models based on normalized CRMs (e.g.Dirichlet process mixture model, etc.) for future study.The truncated BNP model can be obtained by truncating the CRM andthen constructing the feature assignment process accordingly, which takesthe following formµK =K∑k=1θkδψkXn | µK =K∑k=1Zn,kδψki.i.d.∼ LP(G,µK), i = 1, . . . , NYn | Xn indep∼ F(K∑k=1zn,kψk).Since features and feature assignments are latent in this case, they bothneed to be inferred during posterior inference. As we constrain the BNP705.1. TRUNCATED BNP MODELS WITH LATENT FEATURE ASSIGNMENTmodel to the finite dimensional space of the truncated part, the target of theposterior inference is the marginal distribution of the truncated part in theoriginal BNP model. Specifically, we are interested in the following posteriordistributionp(ϱ, θ1:K , ψ1:K , Z1:N,1:K | Y1:N ) (5.3)∝ p(ϱ) ·K∏k=11[θK ≤ θk]νϱ(θk)e−νϱ[θK ,∞) ·K∏k=1h(ψk) ·N∏n=1K∏k=1g(zn,k | θk)· f(Y1:N | ϱ, θ1:K , z1:N,1:K , ψ1:K),where we have assumed the inverse Le´vy representation and p(ϱ) representsthe probability density function of the prior distribution of the parametersof ν.As with many commonly adopted truncated inference algorithms, weperform posterior inference for only the truncated CRM, µK , and its corre-sponding feature assignments, Z1:N,1:K . However, in Eq. (5.3) we integrateout the tail, µK+ of the CRM, and the feature assignments Z1:N,K+ in the tailfrom the joint distribution, rather than discarding the whole part and usinga completely parametric model. Therefore, the posterior model specifiedin Eq. (5.3) does not eliminate the possibility of active features in tail foreach observation Yn. This indicates that Eq. (5.3) represents the marginaldistribution of the truncated part in the true posterior distribution, insteadof the posterior distribution of the truncated model. Similar approach istaken in the variational inference algorithm for the Dirichlet process modelproposed by Blei and Jordan (2006), where only the variational distributionis truncated under the stick-breaking representation of the DP (Sethuraman,1994).The last term f(Y1:N | ϱ, θ1:K , z1:N,1:K , ψ1:K) in Eq. (5.3) takes the formf(Y1:N | ϱ, θ1:K , ψ1:K , z1:N,1:K):=∫f(Y1:N | x1:N ) · g(z1:N,K+ | θ1:K , ϱ) · h(ψK+) dψK+ dz1:N,K+,where we abuse the notation and denote xn =∑∞k=1 zn,kδψk ,g(z1:N,K+ | θ1:K , ϱ) :=∫ν(θK+ | θK , ϱ)N∏n=1∞∏k=K+1g(zn,k | θk) dθK+,and f(Y1:N | x1:N ) is the joint likelihood of observations under the true715.2. APPROXIMATE INFERENCE VIA TRUNCATIONmodel, which can be decomposed into a product form as followsf(Y1:N | x1:N ) =N∏n=1f(Yn |K∑k=1zn,kψk +∞∑k=K+1zn,kψk).To facilitate the theoretical analysis of the truncated model, we furtherassume that the likelihood function of the observation Yn takes the formf(Yn | xn) = exp{f˜(Yn −K∑k=1zn,kψk −∞∑k=K+1zn,kψk)}, (5.4)for some negative valued function f˜ < 0. We note that this assumption ismainly for the convenience of theoretical analysis and the form of f(Yn | xn)specified in Eq. (5.4) are common in practice.5.2 Approximate inference via truncationInference for the marginal distribution of the truncated part of the posteriorBNP requires us to evaluate the unnormalized probability density specifiedin Eq. (5.3). These parts in the first line of Eq. (5.3) are priors of thetruncated model and feature assignments in the truncated part, which arestraightforward to evaluate. However, the (partially) marginal likelihoodf(Y1:N | ϱ, θ1:K , ψ1:K , z1:N,1:K) given only features in the truncated CRM µKis hard to evaluate, since the i.i.d. relationship among observations is lostwithout conditioning on features in the tail µK+.To develop a viable inference algorithm, we propose to perform pos-terior inference of the truncated part using an approximate distribution,Πˆµ(ϱ, θ1:K , ψ1:K , z1:N,1:K | Y[N ]), to the marginal posterior, which can beevaluated easily. In particular, we replace the intractable likelihood termf(Y1:N | ϱ, θ1:K , ψ1:K , z1:N,1:K) with an approximation fˆ defined as follows:fˆ(Y1:N | θ1:K , ψ1:K , z1:N,1:K):=f(Y1:N | θ1:K , ψ1:K , z1:N,1:K , z1:N,K+ = 0, ψK+)=N∏n=1exp{f˜(yn −K∑k=1zn,kψk)}. (5.5)It can be seen from Eq. (5.5) that the approximate likelihood fˆ(Y1:N |θ1:K , ψ1:K , z1:N,1:K) is simply the likelihood of observations conditioned onz1:N,K+ = 0. This is the likelihood of the fully parametric truncated model725.3. POSTERIOR TRUNCATION ERRORwidely adopted in posterior inference for BNP models in practice. Substi-tuting the approximate likelihood specified in Eq. (5.5) into the marginalposterior specified in Eq. (5.3) leads to the approximate distribution thatwill be used for approximate posterior inferencepˆ(ϱ, θ1:K , ψ1:K , z1:N,1:K | Y1:N )∝ p(ϱ) ·K∏k=11[θK ≤ θk]νϱ(θk)e−νϱ[θK ,∞) ·K∏k=1h(ψk) ·N∏n=1K∏k=1g(zn,k | θk)· fˆ(Y1:N | θ1:K , z1:N,1:K , ψ1:K).5.3 Posterior truncation errorThe truncated approximate distribution in Eq. (5.8) is fully parametric andeasy to evaluate. To make the approximate posterior useful, we analyzethe approximation error of the truncated model in this section, which couldalso be used to select the appropriate truncation level. In particular, wewill quantify the total variation distance between the marginal posterior ofthe truncated part of the true model in Eq. (5.3) and its approximation inEq. (5.8). The total variation distance is bounded by a function of an upperbound and a lower bound of ratio πˆµ/πµ of the unnormalized densities of thesetwo models. We will see that the truncation error derived in Theorem 4.1.2is a special case of the general result we will derive in this section.5.3.1 Approximation error between arbitrary distributionsBefore we proceed to derive the approximation error of the truncated model,we first show that for any pair of distributions Πµ and Πˆµ that are defined onthe same space, the total variation distance between them can be boundedas long as the ratio of their unnormalized density functions πµ/πˆµ can bebounded from above and below. Similarly, the forward and backward KLdivergence between these two distributions can also be bounded under thesame condition. Furthermore, we show in Theorem 5.3.1 that the totalvariation distance (and KL divergences) can be estimated using only samplesfrom Πˆµ, which is the distribution that is easier to sample from.Theorem 5.3.1. Let Πµ be a distribution with density proportional to πµ(s),and Πˆµ be another distribution with density proportional to πˆµ(s). If there arefunctions rℓ(s) and ru(s), such that rℓ(s) ≤ πµ(s)πˆµ(s) ≤ ru(s), the total variation735.3. POSTERIOR TRUNCATION ERRORFigure 5.1: An illustration of the posterior truncation errordistance between Πµ and Πˆµ can be bounded as followingDTV(Πˆµ,Πµ)≤ 12EΠˆµ[max{| ruEΠˆµ [rℓ]− 1|, | rℓEΠˆµ [ru]− 1|}], (5.6)where EΠˆµ [·] represents the expectation under the distribution Πˆ.Furthermore,DKL(Πˆµ||Πµ)≤ logEΠˆµ [ru]− EΠˆµ [log rℓ]andDKL(Πµ||Πˆµ)≤ − logEΠˆµ [rℓ] +1EΠˆµ [rℓ]Eˆ[ru log ru].We first note that the conclusion in Theorem 5.3.1 applies to any distri-butions Πµ and Πˆµ on the same space, not only to a posterior distributionand its approximation. We also note that the result in Theorem 5.3.1 willbe trivial if the ratio πµ(s)/πˆµ(s) can be evaluated. Therefore the result willonly be useful in practice when one distribution, say π(s), is hard to evaluateand then we can use the upper and the lower bound of the ratio, as well assamples from Πˆµ to estimate the total variation distance (or KL divergences)between the two distributions.745.3. POSTERIOR TRUNCATION ERRORProof of Theorem 5.3.1Here we derive the total variation distance between Πµ and Πˆµ. We denoteCµ and Cˆµ as the normalization constant of πµ and πˆµ correspondingly, thenCµ =∫πµ =∫πµπˆµπˆµ ≤∫ruπˆµCˆµCˆµ = Cˆµ EΠˆµ [ru];andCµ =∫πµ =∫πµπˆµπˆµ ≥∫rℓπˆµCˆµCˆµ = Cˆµ EΠˆµ [rℓ].So1EΠˆµ [ru]≤ CˆµCµ≤ 1EΠˆµ [rℓ].Therefore it can be shown thatDTV(Πˆµ,Πµ)=12∫| πµCµ− πˆµCˆµ| = 12∫πˆµCˆµ|πµπˆµCˆµCµ− 1|≤ 12∫πˆµCˆµmax{| ruEΠˆµ [rℓ]− 1|, | rℓEΠˆµ [ru]− 1|}=12EΠˆµ[max{| ruEΠˆµ [rℓ]− 1|, | rℓEΠˆµ [ru]− 1|}].The forward and the backward KL divergence between Πµ and Πˆµ canbe derived similarly. In particular,DKL(Πˆµ||Πµ)=∫πˆµCˆµlogCµCˆµπˆµπµ= − log CˆµCµ−∫πˆµCˆµlogπµπˆµ≤ − log 1EΠˆµ [ru]−∫πˆµCˆµlog rℓ= logEΠˆµ [ru]− EΠˆµ [log rℓ].Similarly,DKL(Πµ||Πˆµ)=∫πµCµlogCˆµCµπµπˆµ= logCˆµCµ+∫πµCµlogπµπˆµ= logCˆµCµ+∫πˆµCˆµπµπˆµCˆµCµlogπµπˆµ≤ − logEΠˆµ [rℓ] +1EΠˆµ [rℓ]EΠˆµ [ru log ru].755.3. POSTERIOR TRUNCATION ERRORFig. 5.1 illustrates the total variation distance between πµ and πˆµ. Herewe do not require that πˆµ is an upper bound of πµ as we did in Fig. 3.2.The assumption in Theorem 4.1.2 is equivalent to assume that ru = 1.Additionally, there we also assume that truncation level is large enoughsuch that πµ/πˆµ > 1− ϵ with high probability under the truncated modelΠˆ. This condition can be translated into rℓ = (1 − ϵ)1A for a subset A ofthe state space. By substituting the choice of ru and rℓ corresponding toassumptions of Theorem 4.1.2 into Theorem 5.3.1, we can obtain a similarposterior approximation error with that in Theorem 4.1.2. The result isstated in Corollary 5.3.2, and the proof is in Appendix A.4.Corollary 5.3.2. Suppose that A is a subset of the sample space such thatΠˆµ(A) = 1− η for some 0 < η < 1. If ru = 1, which indicates that πˆµ > πµ,and rℓ = (1 − ϵ)1A for some 0 < ϵ < 1 then the total variation distancebetween Πµ and Πˆµ is bounded,DTV(Πˆ,Π)≤ 12(ϵ1− ϵ + 2η).5.3.2 Approximation error for latent feature BNP modelsTheorem 5.3.1 quantifies the total variation distance between any two dis-tributions given an upper bound and a lower bound of the ratio of theirunnormalized probability densities. Also, the total variation distance (andthe KL divergences) can be estimated using samples from the distributionthat is easier to handle. To estimate the total variation distance betweenthe approximate truncated model Πˆµ specified in Eq. (5.8) and the marginalposterior of the truncated part Πµ in Eq. (5.3), it is necessary to find properru and rℓ for the ratio πµ/πˆµ.In this section we derive an upper bound ru and a lower bound rℓ for theratio πµ/πˆµ in truncated inference for the BNP model with latent featureassignments, and use it to quantify the posterior approximation error of thetruncated model.Recall that the marginal likelihood in the BNP model takes the formf(Y1:N | ϱ, θ1:K , ψ1:K , z1:N,1:K)=E[exp{N∑n=1f˜(yn −∞∑k=1zn,kψk)}| ϱ, θ1:K , ψ1:K , z1:N,1:K]. (5.7)For simplicity of notation, we denote rn = yn −∑Kk=1 zn,kψk which could beinterpreted as the “residual” of observation Yn that is unexplained by the765.3. POSTERIOR TRUNCATION ERRORfeatures from the truncated model. Lemma 5.3.3 finds an upper and a lowerbound for the ratio of unnormalized densities.Lemma 5.3.3. Let rn = yn −∑Kk=1 znkψk. DenoteLθK =∫ θK0(1− g(0 | θ))νϱ(dθ).Assume that f˜ ≤ 0 in Eq. (5.7), and that there are functions u0(θK), u1(θK)such that for all j = 1, . . . , N ,∫ θK0θj νϱ(dθ) ≤ u0(θK) · u1(θK)j .Denote UθK := u0(θK)u1(θK), then the ratio πµ/πˆµ can be bounded frombelow and above as follows:exp(−N · LθK ) ≤πµπˆµ≤N∏n=1(1 + e−f˜(rn) · 4N5 log(N + 1)· UθK).A closer look at LθK tells us that we can expect LθK converges to zeroapproximately in the order of θ1−αK for common CRMs. For example, LθK ∝∫ θK0 θ·θ−1−α(1−θ)λ+α−1dθ for the beta process. We have seen from Eq. (2.11)and Eq. (2.10) that θK itself typically converges to zero exponentially fordense model with α = 0, and it converges to zero polynomially for sparsemodel with α > 0. Therefore, we can expect that exp(−N · LθK ), the lowerbound of the ratio converges to one approximately in the rate of N · exp(−K)for dense models, and in the rate N ·K−(1−α)/α for sparse models, whichcould be very slow when α gets larger.Next we consider the convergence rate of the upper bound on πµ/πˆµ. Asthe truncation level K gets large, we can expect that all observations are wellexplained by existing features and that the approximate likelihood stabilizesand be bounded quickly with as K increases. Therefore, the tightness of theupper bound on the likelihood ratio in Lemma 5.3.3 will be mainly determinedby the rate that UθK converges to zero. We will see from examples belowthat for common CRMs including the beta process and the gamma process,UθK converges to one approximately in the rate N2/ log(N) · exp(−K) fordense processes and in the rate N2/ log(N) ·K−(1−α)/α for sparse ones. Thisimplies that to ensure that the truncated model Πˆµ has small approximationerror, the truncation level K should be in the order of log(N2/ log(N)) forthe dense model, and in the order of (N2 log(N))α/(1−α) for sparse models.775.3. POSTERIOR TRUNCATION ERRORBeta-Bernoulli processSuppose that ν is the beta process BP(γ, λ, α) specified in Eq. (2.2). In thiscase,∫ θK0θj νϱ(dθ) = γΓ(λ+ 1)Γ(1− α)Γ(λ+ α)∫ θK0θj · θ−1−α(1− θ)λ+α−1dθ≤ γ Γ(λ+ 1)Γ(1− α)Γ(λ+ α)∫ θK0θj−1−αdθ≤ γ Γ(λ+ 1)Γ(2− α)Γ(λ+ α)θ−αK · θjKTo fit it in Lemma 5.3.3, we can choose u0(θK) =γΓ(λ+1)Γ(2−α)Γ(λ+α)θ−αK , andu1(θK) = θK . Therefore UθK =γΓ(λ+1)Γ(2−α)Γ(λ+α)θ1−αK for the beta process .Similar to the analysis for the lower bound, we can expect that the upperbound on πµ/πˆµ converges to one approximately in the rate N2/ log(N) ·exp(−K) for dense processes, and in the rate N2/ log(N) · K−(1−α)/α forsparse ones.Gamma-Poisson processSuppose that ν is the gamma process GP(γ, λ, α) specified in Eq. (2.3). Inthis case, ∫ θK0θj νϱ(dθ) = γλ1−αΓ(1− α)∫ θK0θi · θ−1−αe−λθdθ≤ γ λ1−αΓ(1− α)∫ θK0θi−1−αdθ≤ γ λ1−αΓ(2− α)θ−αK · θjK .To fit it in Lemma 5.3.3, we set u0(θK) = γλ1−αΓ(2−α)θ−αK , and u1(θK) = θK .Therefore UθK =γλ1−αΓ(2−α)θ1−αK for the gamma process.In analogy to the beta process, we can also see that πµ/πˆµ convergesto one approximately in the rate N2/ log(N) · exp(−K) for dense processes,and in the rate N2/ log(N) ·K−(1−α)/α for sparse ones.785.4. ADAPTIVELY TRUNCATED INFERENCE5.4 Adaptively truncated inferenceGiven the error bound of the truncated model derived in Theorem 5.3.1, wecan use it to determine the truncation level K when performing approximateinference using the truncated model. Similar to the adaptive truncatedinference for BNP models with explicit labels in Chapter 4, we can take thefollowing adaptive truncation procedure:1. Pick a value of K and desired total variation error level ξ ∈ (0, 1).2. Obtain samples from the truncated posterior.3. Estimate the upper bound of the total variation error in Theorem 5.3.1.4. If the bound exceeds ξ, increase K and return to step 2.5. Otherwise, return the samples.In practice, we can start with a relatively small truncation level K. If thecurrent truncation level is not large enough, we choose the next truncationlevel by extrapolating the total variation bound using samples from the priordistributions, similar to what we did in the adaptively truncated inferencewith explicit labels.5.5 Adjusted Gibbs sampling algorithmsIn this section, we empirically analyze the properties of the posterior distribu-tions of the truncated model, as well as the properties of the total variationerror bound. We experiment with synthetic datasets generated from thebeta-Bernoulli processes together with the Gaussian likelihood. In particular,the synthetic datasets are generated as followsθk ∼ BP(γ = 1, λ = 1, α),ψki.i.d.∼ N (ψk | 0,Σψ),zn,k ∼ Bern(θk),yn ∼ N (yn |∑kzn,kψk,Σy), n = 1, . . . , N,where Σψ = 5× I2, Σy = 0.1× I2, and I2 is the two dimensional diagonalmatrix. To demonstrate the difference between the dense model and thesparse model, we perform truncated inference for two models, includinga dense one with α = 0 and a sparse one with α = 0.3. The number ofobservations is set to N = 1000 in both experiments.795.5. ADJUSTED GIBBS SAMPLING ALGORITHMSThe Gibbs algorithmWe apply the Gibbs sampling algorithm to perform posterior inference forthe truncated model. We fix the hyperparameters, ϱ = (γ, λ, α), of the betaprocess, as well as the variance parameter of the prior distribution of Ψ, andthe variance parameter of the distribution of the observation Y . We onlyperform posterior inference for Θ1:K , Z1:N,1:K , and Ψ1:K using the truncatedmodel. Given a truncation level K, the density of the truncated posteriortakes the form:πˆ(θ1:K , ψ1:K , z1:N,1:K | Y1:N ) (5.8)∝K∏k=1νϱ(θk)1[θk+1 ≤ θk ≤ θk−1] · e−νϱ[θK ,∞) ·K∏k=1h(ψk) ·N∏n=1K∏k=1g(zn,k | θk)·N∏n=1fˆ(Yn | θ1:K , zn,1:K , ψ1:K).Based on Eq. (5.8), the Gibbs algorithm will sample from the truncatedmodel by updating each term in πˆ(θ1:K , ψ1:K , z1:N,1:K | Y1:N ) iteratively. Inparticular, the conditional distribution of each term are as follows:1. Density of θk conditioning on other variable, for k = 1, . . . ,K − 1:p(θk | ·) ∝ 1[θk+1 ≤ θk ≤ θk−1] · νϱ(θk) · g(z1:N,k|θk); (5.9)2. Density of θK conditioning on others:p(θK | ·) ∝ 1[0 ≤ θK ≤ θK−1] · νϱ(θK) · e−νϱ[θK ,∞) · g(z1:N,K | θK);(5.10)3. Conditional probability of Zn,k for k = 1, . . . ,K, and n = 1, . . . N :p(zn,k | ·) ∝ g(zn,k|θk) · fˆ(Yn | zn,1:K , ψ1:K); (5.11)4. Density of Ψk conditioning on other variables, for k = 1, . . .K:p(ψk | ·) ∝ h(ψk) · fˆ(Y1:N | z1:N,1:K , ψ1:K); (5.12)where we use the notation p(x | ·) to denote the conditional probabilitydensity (or mass) function of a variable x conditioned on all other variables.More specifically, we apply the slice sampler (Neal, 2003) to sample θkfrom its conditional distribution. The conditional distribution of the binary805.5. ADJUSTED GIBBS SAMPLING ALGORITHMSfeature allocation variable, Zn,k, maintains a Bernoulli distribution. Addition-ally, the updates of all feature allocations of one one observation, Zn,1:K forn = 1, . . . , N , are embarrassingly parallel over n ∈ 1, . . . , N . As we use theconjugate prior and likelihood distribution for features Ψ1:K , the conditionaldistribution of ψk are Gaussian as well. In fact, the features Ψ1:K jointlyfollows a matrix Gaussian distribution conditioned on all other variables.Therefore it is also viable to update Ψ1:K together in each iteration of theGibbs sampling. However, experiment results indicate that updating Ψ1:Kjointly does not improve the mixing efficiency of the Gibbs sampler. Hencein our experiments we sample Ψk individually to maintain the generality ofthe Gibbs sampler.As the goal of the posterior inference is to discover the distribution offeatures in the population, we apply the Gibbs sampling and examine thesample means of (log10(θk))k∈1:K under different truncation levels.Fig. 5.2 presents the behavior of sample means of (log10(θk))k∈1:K forsamples from the Gibbs sampler. It is a characteristic of the inverse Le´vyrepresentation that the traces of θk for both models decrease monotonicallywith k. It is also as expected that θk from the sparse model with α = 0.3decrease much slower with k than that from the dense model with that hasα = 0.However, one obvious problematic pattern in Fig. 5.2 is that the traces ofθk of all truncation levels contain a lot flat regions, and this pattern happensfor both the sparse and the dense model. Some closer examination of theGibbs sampling algorithm indicates that the flat regions in each curve arethe result of that inactive features having Z1:N,k = 0 take places in front ofactive feature having∑N1 Zn,k > 0. Since the inverse Le´vy representationrequires that θk ≥ θk+1, the rates of those inactive features are forced tobe as large as the active features that come after them. As the likelihoodof observations are not affected by the order of features, it is hard for theGibbs sampler to correct the orders by itself especially when these disorderedfeatures fit the observations well. Similarly, those θk’s of active features arealso affected by this characteristic of the inverse Le´vy representation, and itcauses the flat parts at the beginning of all traces of log10(θk) in Fig. 5.2.The Gibbs algorithm with feature swapsThese problematic flat regions motivates us to add two extra steps that swapfeatures and their assignments in each iteration of the Gibbs sampler. In thefirst swap, we randomly pick an inactive feature and one active feature andpropose an Metropolis Hastings (MH) move that swaps their order. Similarly,815.5. ADJUSTED GIBBS SAMPLING ALGORITHMS(a) α = 0. (b) α = 0.(c) α = 0.3. (d) α = 0.3.Figure 5.2: Sample means of log10(θk) for samples generated from the regularGibbs sampler. Each curve represents a trace of log10(θk) for one truncationlevel K. Samples of different truncation levels are generated independentlyfrom the same model. Subplots (a) and (b) are samples from the densemodel with α = 0, and the subplots (c) and (d) are samples from the sparsemodel with α = 0.3. Subplots (b) and (d) are magnified parts of the subplots(a) and (c) respectively, for the purpose to show the pattern of samples ofθk’s clearly. It can be seen from the plots that sample means of log10(θk)for samples generated from the standard Gibbs sampler, demonstrating theproblematic flat regions that appear to due the inability of the sampler toreorder the columns.we randomly select two active features and swap their orders with an MHstep in the second swap.More specifically, we denote si =∑Nn=1 zn,i represents the number ofobservations that express the feature ψi. Let I represent the set of indices ofactive features under the inverse Le´vy representation, and let J be the setof indices of inactive features in the first swap and as the set of indices of825.5. ADJUSTED GIBBS SAMPLING ALGORITHMSactive features other than ψi in the second swap. In each step, the featuresare selected and swapped as follows:1. Randomly choose two features ψi, ψj and their feature assignmentvectors Z1:N,i, Z1:N,j :i ∼ Unif{i1, i2, . . .}i∈I , j ∼ Unif{j1, j2, . . .}j∈J ; (5.13)2. Swap the positions of ψi with ψj and Z1:N,i with Z1:N,j with theacceptance probability:pswap = min{1,θsji (1− θi)N−sjθsij (1− θj)N−si}. (5.14)The acceptance probability in Eq. (5.14) does not include the likelihood termsince the orders of features do not affect the likelihood of observations. Themotivation for these swap steps is to put active features in front of inactiveones, and put more-active features in front of less-active ones. While thesampling distribution is not changed by these two MH steps, we find in ourexperiments that the order of features are usually very stable once moreactive features are swapped to the front places. It is very unlikely that moreactive features be swapped to the back of less active ones in these MH steps.Fig. 5.3 plots traces of sample means of log10(θk) generated from theadjusted Gibbs sampler with two additional column swap steps. It canbe seen that the extra swap steps improve the performance Gibbs samplersignificantly for both the dense and sparse model. As we are working withthe synthetic data and know the true values of θ’s, we further compare thesample means of log10(θ) to their true values that are used to generate thesynthetic data. It can be seen that the estimates of the rates of featuresare not accurate enough for the Gibbs sampler even with these additionalfeature swap steps.Annealed Gibbs algorithm with feature swapsGiven the high dimension of the truncated model, it is likely that theinsufficiency of the Gibbs algorithm originates from the multimodal nature ofthe target model. To address this problem, we add one additional warm-upstage using a fixed annealing scheme at the beginning of the Gibbs sampler.This is achieved by replacing the approximate likelihood fˆ(Yn | zn,1:K , ψ1:K)with the annealed likelihood fˆ(Yn | zn,1:K , ψ1:K)at , using an annealing factorat increases from 0 to 1 in the warm-up stage.835.5. ADJUSTED GIBBS SAMPLING ALGORITHMS(a) α = 0. (b) α = 0.(c) α = 0.3. (d) α = 0.3.Figure 5.3: Sample means of log10(θk) for samples generated from the Gibbssampler with column swaps. Each curve represents a trace of log10(θk) forone truncation level K. Samples of different truncation levels are generatedindependently from the same model. The black dotted lines are curves ofthe true values of log10(θk) used to generate the synthetic data. Subplots (a)and (b) are samples from the dense model with α = 0, and the subplots (c)and (d) are samples from the sparse model with α = 0.3. Subplots (b) and(d) are magnified parts of the subplots (a) and (c) respectively. It can beseen that these flat regions in the original Gibbs sampler is alleviated a lotby these two additional column swap steps. However, the estimates of theθk’s are still not close enough to the true values, implying that there mightbe multimodal issue due to the high dimension of the truncated model.Only the conditional distributions of Z1:N,1:K and Ψ1:K are influenced bythe likelihood term. Therefore, in each iteration of the Gibbs sampling, wesubstitute Eq. (5.11) and Eq. (5.12) with the following annealed conditionaldistributions:845.5. ADJUSTED GIBBS SAMPLING ALGORITHMS1. Annealed conditional probability of Zn,k:p(zn,k | ·) ∝ g(zn,k|θk) · fˆ(Yn | zn,1:K , ψ1:K)at ; (5.15)2. Annealed conditional probability of Ψk for k = 1, . . .K:p(ψk | ·) ∝ h(ψk) · fˆ(Y1:N | z1:N,1:K , ψ1:K)at . (5.16)Note that the annealed posterior of ψk is still conjugate to the Gaussianprior. We choose the quadratic annealing scheme in the experiments, whichhelped the sampler converges to the correct position without in a relativelyfast speed. Some preliminary results indicated that the linear annealingscheme increases the temperature at too fast and provides little improvementto the Gibbs sampler, while polynomial schemes with orders greater thantwo provide significant improvement.Fig. 5.4 plots traces of means of log10(θk) sampled with the Gibbs samplerwith two additional column swap steps and the warm-up stage using a fixedannealing scheme. It can be seen from the plots that the warm-up stageimproves the efficiency of the Gibbs sampler significantly. Even thoughcorrectly discovering θ1:K alone is not enough to determine the correctness ofthe whole posterior, the distribution of features in the population is usuallythe most interesting characteristics in applying the Bayesian nonparametriclatent feature models. Therefore we will use this adjusted Gibbs sampler toestimate posterior approximation errors of truncated models. Algorithm 1presents the details of the adjusted Gibbs sampling algorithm with featureswap steps and a warm-up stage using the quadratic annealing scheme.Analysis of the truncation errorUsing samples from the Algorithm 1, we estimate the lower and the upperbound of πµ/πˆµ and show empirically how the estimated truncation errorbound derived in Theorem 5.3.1 decreases with the truncation levels.It can be seen from Eq. (5.6) that the upper bound on DTV(Πˆµ,Πµ) ismainly determined by the tightness of the lower bound rℓ, and the upperbound ru, of πµ/πˆµ, which are determined by the rate that θK converges to0. So we first show empirically how θK converges to 0 during the truncatedinference. Fig. 5.5 demonstrates the decaying trend of θK with the truncationlevel K, for both the dense beta-Bernoulli model with α = 0 and the sparsemodel with α = 0.3. As expected, it can be seen that log10(θK) decreaseslinearly when α = 0, while it decreases much slower when α = 0.3.855.5. ADJUSTED GIBBS SAMPLING ALGORITHMS(a) α = 0. (b) α = 0.(c) α = 0.3. (d) α = 0.3.Figure 5.4: Sample means of log10(θk) for samples generated from the Gibbssampler with column swaps and warm-up stage with a fixed annealing scheme.Each curve represents a trace of log10(θk) for one truncation level K. Samplesof different truncation levels are generated independently from the samemodel. The black dotted lines are curves of the true values of log10(θk) usedto generate the synthetic data. Subplots (a) and (b) are samples from thedense model with α = 0, and the subplots (c) and (d) are samples from thesparse model with α = 0.3. Subplots (b) and (d) are magnified parts of thesubplots (a) and (c) respectively. It can be seen that the annealed stagegenerates significant improvement than the Gibbs algorithm with featureswap steps.We then show the decreasing trend of the estimated total variationbound estimated using samples from Algorithm 1. Fig. 5.6 shows that whenα = 0 the estimated bound of DTV(Πˆµ,Πµ) decreases very fast with thetruncation level K. In the specific example, we see that the estimatedDTV(Πˆµ,Πµ) is smaller than 10−4 when K > 30, indicating that in thecase of the dense model, even the truncated model with a relatively small865.5. ADJUSTED GIBBS SAMPLING ALGORITHMSAlgorithm 1: Gibbs sampling algorithm (with fixed annealing)Set the truncation level K and hyperparameters ϱ = (α, γ, λ) of BP,and hyperparameters of H and F ;Se the number of iterations Ts, and the warm-up iterations Tburn;Initialize θ1:K , Z1:N,1:K , and Ψ1:K ;for t ∈ 1, . . . , Tburn + Ts do// Set annealing factorat = min(1, (t/Tburn)2)// Update θ1:K :for k ∈ 1, . . . ,K − 1 doUpdate θk via the slice sampler targeting Eq. (5.9);Update θK via the slice sampler targeting Eq. (5.10);// Update Z1:N,1:K ;for k ∈ 1, . . . ,K and n ∈ 1, . . . , N doSample Zn,k from the Bernoulli distribution with probabilityproportional to Eq. (5.15);// Update ψk;for k ∈ 1, . . .K doSample Ψk from Gaussian distribution targeting Eq. (5.16);// Swap features and feature assignments:Swap an active feature with an inactive one following Eq. (5.13)and Eq. (5.14);Swap two active features following Eq. (5.13) and Eq. (5.14).return Samples of θ1:K , Z1:K , Ψ1:K after the warm-up stage.truncation level provides a good enough approximation to the true modelfor posterior inference. The experimental results also implies that the boundof DTV(Πˆµ,Πµ) in Theorem 5.3.1 is tight enough to be used in selecting thetruncation level when the model is dense.Fig. 5.7 shows that when α = 0.3 the estimated bound of DTV(Πˆµ,Πµ)decrease much slower than that when α = 0. In fact, the estimated totalvariation error bound greater than 1 even with the truncation level K = 1000,indicating that this bound can only be used to choose truncation level whenK > 1000. However, the fast decreasing rate of the estimated bound ofDTV(Πˆµ,Πµ) indicates that we can still use it to quantify the approximationquality of the sparse truncated model when K is large enough, even thoughthe estimation of the approximation error might not be accurate enough.875.5. ADJUSTED GIBBS SAMPLING ALGORITHMS(a) α = 0. (b) α = 0.3.Figure 5.5: Sample mean of log10(θK) varying with the truncation level K,for both the dense model with α = 0 and the sparse model with α = 0.3.Samples are generated from the Gibbs sampler with column swaps andwarm-up stage with a fixed annealing scheme.(a) α = 0. (b) α = 0.Figure 5.6: TV error bound of the truncated model varying with the trunca-tion level K. Samples are generated from the Gibbs sampler with columnswaps and warm-up stage with a fixed annealing scheme. The subplot (b)is a magnified part of subplot (a) where the error bound is smaller than 1.The shaded regions are 95% Bootstrap confidence intervals.885.5. ADJUSTED GIBBS SAMPLING ALGORITHMS(a) α = 0.3.Figure 5.7: TV error bound of the truncated model varying with the trunca-tion level K. Samples are generated from the Gibbs sampler with columnswaps and warm-up stage with a fixed quadratic annealing scheme. Theshaded regions (which is almost invisible) is the 95% Bootstrap confidenceintervals. When K = 1000, log10(DTV(Πˆµ,Πµ)) ≈ 0.1.89Chapter 6ConclusionModel selection is a persistent problem in statistical inference, and selectinga model with appropriate complexity is necessary to avoid under-fittingor over-fitting the dataset and to discover the underlying structures andpatterns in the observations. However, model selection is not trivial. Inparametric modeling, this usually involves fitting multiple models of differentcomplexities from a limited family of parametric models, which is limitedand inefficient. The problem is only harder in the face of exploding data invarious fields.Bayesian nonparametric models provide a modeling framework wheremodels lie in an infinite dimensional parameter space but only invoke afinite subset of the parameters, which grows with the amount of data beingfitted. The flexibility of BNP models has made them popular in various fields,including biology, genomics, language models, finance, and recommendersystems. It is known that many popular BNP models have completely randommeasures (CRMs) or normalized completely random measures (NCRMs) astheir prior distributions. Therefore, in this thesis, we focus on the posteriorinference for BNP models with CRM priors.The posterior inference for BNP models is notoriously hard becauseBNP models have an infinite number of parameters. To avoid the infinitenumber of parameters that come with BNP models and keep only the activeparameters in the posterior distribution that are related to the data, previousinference methods either marginalize out the prior distributions of the BNPmodels or truncate the BNP model into a parametric model with only finiteparameters. The first approach is constrained to BNP models with conjugatepriors and likelihoods, breaking the i.i.d. relationship between observationsand causing inefficient posterior inference. The second approach does nothave limitations over the choice of CRM priors or the likelihood process andprovides efficient algorithms that can be easily applied to general BNP models.However, truncation induces truncation errors, which should be quantifiedand controlled to ensure that the truncated model is a good approximationto the true posterior. Previous work on truncated inference for BNP modelsbased on CRMs determined truncation levels either according to the prior90truncation error or based on statistics of the truncated model. Neitherapproach measures the approximation error between the truncated modeland the target posterior distribution.In this thesis, we first extend BNP models by proposing BNP mod-els based on transformations of CRMs. As an example, we propose theedge-exchangeable network model built on the self-product of CRMs. Wedeveloped truncated simulation algorithms for edge-exchangeable networkmodels and analyzed their prior truncation error. We then analyzed theposterior approximation error of the truncated edge-exchangeable networkmodel. The approximation error between the truncated model and theposterior network model can be estimated using prior approximation er-ror and samples from the truncated model. With the posterior truncationerror, we developed an adaptively truncated inference algorithm for theedge-exchangeable network model. Starting with a small truncated model,the algorithm samples from the truncated model and estimates the trun-cation error in each round, gradually increasing the truncation level untilthe desired posterior truncation error is achieved. The expression of theposterior truncation error developed for the network model can be appliedto general BNP models when feature assignment labels are observed.In Chapter 5, we studied truncated inference methods for BNP modelswith latent feature assignments. We derived an upper bound for the totalvariation distance between any distributions, given any upper bound andlower bound of the ratio of the unnormalized densities of the two distribu-tions. The error can be estimated using only samples from the approximatedistribution that is easy to sample from. We also showed that the posteriorapproximation error derived in Chapter 4 is a special example of this approx-imation error. We derived the ratio of the density of the posterior of thetruncated BNP models with latent features to that of the true posterior andused it to find the posterior approximation error of the truncated model.The goal of this thesis is to develop truncated inference methods for BNPmodels that can be easily applied to general BNP models and also provideerror guarantees. We first note that there are many BNP models not coveredin this work, and the theoretical analysis in this thesis can be extended tomore BNP models. Such examples include BNP models based on normalizedCRMs and hierarchical BNP models based on CRMs. The analysis ofposterior approximation error using ratios of unnormalized densities couldalso be extended to BNP models built on other representations of completelyrandom measures. One example is the Poisson-Kingman process, which isa generalized CRM obtained by allowing the total mass to come from anindependent distribution, as we have seen in Chapter 2. Another example is91BNP models based on stick-breaking constructions of CRMs and NCRMs.While BNP models based on stick-breaking representation are easy to simulateand truncate, the posterior approximation error is not studied in previouswork.We note that our proposed truncated inference algorithms only involvestandard parametric inference algorithms. Additionally,the algorithms arevery general, and do not require any user-specified mathematical analysis.These properties make it possible to extend the work of this thesis to providegeneral BNP inference algorithms in popular computational frameworks,leveraging their existing computational engines.92Appendix AProofsA.1 Equivalence between rejection and inverseLe´vy representationsProof of Theorem 2.5.1. Denote Tk1 be the first nonzero element that is gen-erated from the rejection representation from Eq. (2.12) and correspondingly,denote Γk1 be the jump of the unit-rate homogeneous Poisson process on R+such that Tk1 = νˆ←(Γk1), where νˆ is the proposal measure in the rejectionrepresentation. Let f be a bounded continuous function. ThenE[f(Tk1)]=E[f(νˆ←(Γk1))]=E ∞∑j=1f(νˆ←(Γj))1[dνdνˆ(νˆ←(Γj)) ≥ Uj)] j−1∏i=11[dνdνˆ(νˆ←(Γi)) < Ui]=E ∞∑j=1f(νˆ←(Γj))dνdνˆ(νˆ←(Γj))j−1∏i=1(1− dνdνˆ(νˆ←(Γi)))=∞∑j=1E[f(νˆ←(Γj))dνdνˆ(νˆ←(Γj))E[j−1∏i=1(1− dνdνˆ(νˆ←(Γi)))∣∣∣∣∣Γj]]Note that given Γj , Γii.i.d.∼ Unif(0,Γj), for i = 1, · · · , j − 1, soE[j−1∏i=1(1− dνdνˆ(νˆ←(Γi)))∣∣∣∣∣Γj]= E[1− dνdνˆ(νˆ←(U)) | Γj]j−1,where Ui.i.d.∼ Unif(0,Γj). Using the change of variable y = νˆ←(u), we obtainE[1− dνdνˆ(νˆ←(U)) | Γj]= 1− 1Γj∫ Γj0dνdνˆ(νˆ←(u))du = 1− 1Γj∫ ∞νˆ←(Γj)dν.93A.1. EQUIVALENCE BETWEEN REJECTION AND INVERSE LE´VY REPRESENTATIONSTherefore, using the same change of variable trick,E[f(Tk1)]=∞∑j=1Ef(νˆ←(Γj))dνdνˆ(νˆ←(Γj))(1− 1Γj∫ ∞νˆ←(Γj)dν)j−1=∞∑j=1∫ ∞0f(νˆ←(γ))dνdνˆ(νˆ←(γ))(1− 1γ∫ ∞νˆ←(γ)dν)j−1γj−1(j − 1)!e−γdγ=∞∑j=1∫ ∞0f(y)(1− νˆ[y,∞)−1ν[y,∞))j−1 νˆ[y,∞)j−1(j − 1)! e−νˆ[y,∞)ν(dy)=∫ ∞0f(y)∞∑j=1(νˆ[y,∞)− ν[y,∞))j−1(j − 1)! e−νˆ[y,∞)ν(dy)=∫ ∞0f(y)e−ν[y,∞)ν(dy).Suppose that θ1 is the first rate generated using the inverse Le´vy representa-tion. ThenE[f(θ1)] =∫ ∞0f(ν←(γ))e−γdγ.Making the change of variable y = ν←(γ), we obtainE[f(θ1)] =∫ ∞0f(y)e−ν[γ,∞)ν(dy) = E[f(Tk1)].Therefore, the first nonzero rate θk1 from the rejection representation hasthe same marginal distribution as the first rate θ1 from the inverse Le´vyrepresentation.We now employ an inductive argument. Suppose that we have shown thatthe marginal distribution of first nonzero M elements ΞM = (Tk1 , · · · , TkM )from the rejection representation has the same marginal distribution as thefirst M elements ΘM = (θ1, · · · , θM ) from the inverse Le´vy representation.To prove the same for M +1 elements, it suffices to show that the conditionaldistribution of TkM+1 given ΞM is equal to the conditional distribution ofθM+1 given ΘM when ΞM = ΘM .94A.1. EQUIVALENCE BETWEEN REJECTION AND INVERSE LE´VY REPRESENTATIONSDenote Γ′j =∑ji=1 e′i, where e′ii.i.d.∼ Exp(1), and U ′i i.i.d.∼ Unif[0, 1]. ThenE[f(TkM+1)|ΞM ] = E[f(νˆ←(ΓkM+1))|ΞM ]= E ∞∑j=1f(νˆ←(ΓkM + Γ′j))1[. . . ]∣∣∣∣∣∣ΞM ,where 1[. . . ] is shorthand for1[. . . ] = 1[dνdνˆ(νˆ←(ΓkM + Γ′j)) ≥ U ′j] j−1∏i=11[dνdνˆ(νˆ←(ΓkM + Γi)) < U′i].Using steps similar to the base case,E[f(TkM+1)|ΞM ]=∞∑j=1E[f(νˆ←(ΓkM + Γ′j))dνdνˆ(νˆ←(ΓkM + Γ′j))E[. . . ]j−1∣∣∣∣ΞM] ,whereE[. . . ] = E[1− dνdνˆ(νˆ←(ΓkM + U))∣∣∣∣Γ′j ,ΞM] U ∼ Unif[0,Γ′j ].Making the change of variable y = νˆ←(ΓkM + u) as before, we obtainE[1− dνdνˆ(νˆ←(ΓkM + U))∣∣∣∣Γ′j ,ΞM] = 1− ∫ Γ′j01Γ′jdνdνˆ(νˆ←(ΓkM + u))du= 1− 1Γ′j∫ νˆ←(ΓkM )νˆ←(ΓkM+Γ′j)dν.Making another change of variables y = νˆ←(ΓkM + γ) in the original integral,95A.2. TRUNCATION ERROR BOUNDS FOR SELF-PRODUCT MEASURESand hence γ = νˆ[y,∞)− ΓkM = νˆ[y, νˆ←(ΓkM )), yieldsE[f(TkM+1)|ΞM ]=∞∑j=1∫ ∞0f(νˆ←(ΓkM + γ))dνdνˆ(νˆ←(ΓkM + γ))E[. . . ]j−1 γj−1(j − 1)!e−γdγ=∞∑j=1∫ νˆ←(ΓkM )0f(y)(1− 1γ∫ νˆ←(ΓkM )ydν)j−1γj−1(j − 1)!e−γν(dy)=∞∑j=1∫ νˆ←(ΓkM )0f(y)(νˆ[y, νˆ←(ΓkM ))− ν[y, νˆ←(ΓkM )))j−1(j − 1)! e−γν(dy)=∫ νˆ←(ΓkM )0f(y)e−ν[y,νˆ←(ΓkM ))ν(dy)=∫ TkM0f(y)e−ν[y,TkM )ν(dy).On the other hand,E[f(θM+1)|ΘM ] = E[f(ν←(ΓM + Γ′1))|ΘM ] =∫ ∞0f(ν←(ΓM + γ))e−γdγ=∫ ν←(ΓM )0f(y)e−ν[y,ν←(ΓM ))ν(dy) =∫ θM0f(y)e−ν[y,θM )ν(dy).Thus the distribution of the (M + 1)th nonzero rate in the rejection repre-sentation TkM+1 given ΞM is equal to the distribution of the (M + 1)th ratefrom the inverse Le´vy representation θM+1 given ΘM when ΞM = ΘM .A.2 Truncation error bounds for self-productmeasuresProof of Lemma 3.4.1. Denote µ˜ = {µ(d), µ(d)K }. Denote the marginal proba-bility mass function (PMF) of EN ∈ Nd and EN,K ∈ Nd as PN and PN,K ,96A.2. TRUNCATION ERROR BOUNDS FOR SELF-PRODUCT MEASURESand denote their PMFs given µ˜ as f(x|µ˜) and fK(x|µ˜) respectively.DTV (PN , PN,K)=12∑x∈Nd∣∣∣PN (x)− PN,K(x)∣∣∣=12∑x∈Nd|E[f(x|µ˜)]− E[fK(x|µ˜)]|≤ 12P(IN ≤ K)∑x∈Nd|E[f(x|µ˜)|IN ≤ K]− E[fK(x|µ˜)|IN ≤ K]|+12P(IN > K)∑x∈Nd|E[f(x|µ˜)|IN > K]− E[fK(x|µ˜)|IN > K]|Conditioned on IN ≤ K, f(x|µ˜) = fK(x|µ˜) under both the independent andcategorical likelihood. SoE[f(x|µ˜)|IN ≤ K]− E[fK(x|µ˜)|IN ≤ K] = 0.By Fubini’s Theorem,DTV (PN , PN,K)≤12P(IN > K)∑x∈Nd( E[f(x|µ˜)|IN > K] + E[fK(x|µ˜)|IN > K])=12P(IN > K)E ∑x∈Ndf(x|µ˜)∣∣∣∣∣∣ IN > K+ E ∑x∈NdfK(x|µ˜)∣∣∣∣∣∣ IN > K=P(IN > K) = 1− P (IN ≤ K) .Proof of Theorems 3.4.4 and 4.1.1. Denote the set of indicesIℓ,K = {i ∈ Nd̸= : 1 ≤ i1, · · · , iℓ ≤ K, K + 1 ≤ iℓ+1, · · · , id <∞}such that i ∈ Iℓ,K indicates that the first ℓ elements of i belong to thetruncation, and the remaining d− ℓ elements belong to the tail. By Jensen’s97A.2. TRUNCATION ERROR BOUNDS FOR SELF-PRODUCT MEASURESinequality,P (IN ≤ K | U1:K ,Γ1:K)= EexpNd−1∑ℓ=0(dℓ) ∑i∈Iℓ,Klog g0 (ϑi) | U1:K ,Γ1:K≥ expNd−1∑ℓ=0(dℓ)E ∑i∈Iℓ,Klog g0 (ϑi) | U1:K ,Γ1:K . (A.1)This equation arises by noting that IN ≤ K if and only if for all i involvingan index ij > K, the count of edge i is 0 after N rounds; the factor(dℓ)accounts for the fact that ϑi =∏dj=1 θij is independent of the ordering ofthe ij .Consider a single term∑i∈Iℓ,K log g0(ϑi) in the above sum. Since weare conditioning on U1:K ,Γ1:K , we have that θi1 , . . . , θiℓ are fixed in theexpectation, and the remaining steps ΓK+1,ΓK+2, · · · are the ordered jumpsof a unit rate homogeneous Poisson process on [ΓK ,∞). By the markingproperty of the Poisson process Kingman (1992), conditioned on ΓK , we havethat (Ui,Γi)∞i=K+1 is a Poisson process on R+ × [ΓK ,∞) with rate measuregu(du)dγ. Thus we apply the Slivnyak-Mecke theorem Last and Penrose(2017) to the remaining d− ℓ indices to obtainE ∑i∈Iℓ,Klog g0 (ϑi)∣∣∣∣∣∣Θ[K]E ∑i∈Iℓ,Klog g0 ℓ∏j=1θijd∏j=ℓ+1τ(Uij ,Γij )∣∣∣∣∣∣Θ[K]=∑i1 ̸=···̸=iℓ≤KE∫Rd−ℓ+log g0 ℓ∏j=1θijd∏j=ℓ+1τ(Uij ,ΓK + γj) d∏j=ℓ+1dγj | Θ[K]=∑L⊆[K]|L|=ℓ∫[ΓK ,∞)d−ℓElog g0∏j∈Lθj ·d−ℓ∏j=1τ(U ′j , γj) | θ1:Kdγ1:d−ℓ,where U ′ji.i.d.∼ gu, Θ[K] := {U1:K ,Γ1:K}. Substitution of this expression intoEq. (A.1) yields the result of Theorem 4.1.1. Next, we consider the boundon the marginal probability P (IN ≤ K). By Jensen’s inequality applied to98A.2. TRUNCATION ERROR BOUNDS FOR SELF-PRODUCT MEASURESEq. (A.1) and following the previous derivation, we find thatP (IN ≤ K) ≥ expNd−1∑ℓ=0(dℓ)E ∑i∈Iℓ,Klog g0 (ϑi) ,andE ∑i∈Iℓ,Klog g0 (ϑi)=∑i1 ̸=···̸=iℓ≤KE∫[ΓK ,∞)d−ℓlog g0 ℓ∏j=1θijd−ℓ∏j=1τ(U ′j , γj) dγ1:d−ℓ .Using the fact that conditioned on ΓK , Γ1:K−1 are uniformly distributed on[0,ΓK ], and that at most one ij can be equal to K,∑i1 ̸=···̸=iℓ≤KE∫[ΓK ,∞)d−ℓlog g0 ℓ∏j=1θijd−ℓ∏j=1τ(U ′j , γj) dγ1:d−ℓ=(K − 1)!(K − 1− ℓ)!EΓ−ℓK ∫SdΓlog g0 d∏j=1τ(Uj , γj) dγ1:d+ ℓ(K − 1)!(K − ℓ)!EΓ−(ℓ−1)K ∫SdΓδγ1=ΓK log g0 d∏j=1τ(Uj , γj) dγ1:d ,where SdΓ := [0,ΓK ]ℓ × [ΓK ,∞)d−ℓ, the first and second terms arise fromportions of the sum where all indices satisfy ij ̸= K and one index satisfiesij = K, respectively.To complete the result, we study the asymptotics of the marginal proba-bility bound. It follows from Eq. (3.4) that IN <∞ almost surely. ThereforelimK→∞EexpNd−1∑ℓ=0(dℓ) ∑i∈Iℓ,Klog g0 (ϑi) = limK→∞P(IN ≤ K) = 1.It then follows from (Campbell et al., 2019, Lemma B.1) and continuousmapping theorem thatd−1∑ℓ=0(dℓ) ∑i∈Iℓ,Klog g0 (ϑi) p→ 0 as K →∞.99A.2. TRUNCATION ERROR BOUNDS FOR SELF-PRODUCT MEASURESSince this sequence is monotonically increasing in K, we have thatBK = Ed−1∑ℓ=0(dℓ) ∑i∈Iℓ,K− log g0 (ϑi)→ 0 as K →∞.A.2.1 Proof of Theorem 3.4.5We first state an useful results which states that if one perturbs the prob-abilities of a countable discrete distribution by i.i.d. Gumbel(0, 1) randomvariables, the arg max of the resulting set is a sample from that distribution.Lemma A.2.1. (Campbell et al., 2019, Lemma 5.2) Let (pj)∞j=1 be a count-able collection of positive numbers such that∑j pj <∞ and let p¯j = pj∑k pk. If(Wj)∞j=1 are i.i.d. Gumbel(0, 1) random variables, then argmaxj∈NWj+log pjexists, is unique a.s., and has distributionargmaxj∈NWj + log pj ∼ Categorical((p¯j)∞j=1).Proof of Theorem 3.4.5. Since the N edges from the categorical likelihoodprocess are i.i.d. categorical random variables, by Jensen’s inequality,P(IN ≤ K) = E [P(IN ≤ K|µ)] = E[P(I1 ≤ K|µ)N] ≥ E [P(I1 ≤ K|µ)]N .Next, since ϑi =∏dj=1 θij , we can simulate a categorical variable withprobabilities proportional to ϑi, i ∈ Nd̸= by sampling d indices (J1, · · · , Jd)independently from a categorical distribution with probabilities proportionalto (θ1, θ2, . . . ) and discarding samples where Jj = Jk for some 1 ≤ j, k ≤ d.Denote θ′k = θk/∑k θk to be the normalized rates, PJ,K =∑Kj=J θ′j , and theevent Q = {Jj ̸= Jk,∀1 ≤ j, k ≤ d}. ThenP(I1 ≤ K|µ) = P(1 ≤ Jj ≤ K, 1 ≤ j ≤ d | Q, µ).Since the normalized rates θ′k are generated from the inverse Le´vy represen-tation, they are monotone decreasing. ThereforeP(1 ≤ Jj ≤ K, 1 ≤ j ≤ d | Q, µ) ≥ P1,K1− 0 ·P2,K1− P1,1 · · ·Pd,K1− P1,d−1≥(Pd,K1− P1,d−1)d.100A.2. TRUNCATION ERROR BOUNDS FOR SELF-PRODUCT MEASURESBy Jensen’s inequality,E[P(1 ≤ Jj ≤ K, 1 ≤ j ≤ d | Q, µ)] ≥ E[Pd,K1− P1,d−1]d.Note that for a categorical random variable J with class probabilities givenby θ′j/(1− P1,d−1), j ≥ d, the quantity Pd,K/(1− P1,d−1) is the probabilitythat d ≤ J ≤ K. So by the infinite Gumbel-max sampling lemma,E[Pd,K1− P1,d−1]= P(d ≤ argmaxj≥d(log θj +Wj) ≤ K), Wji.i.d.∼ Gumbel(0, 1),where we can replace θ′j with the unnormalized θj because the normalizationdoes not affect the argmax. DenotingMK = maxd≤k≤Klog ν←(Γk) +Wk, MK+ = supk>Klog ν←(Γk) +Wk,we have thatP(IN ≤ K) ≥ (1− E[P(MK < MK+|ΓK)])N ·d ,and so the remainder of the proof focuses on the conditional expectation.Conditioned on ΓK ,MKd= max{log ν←(ΓK) +WK , maxd≤k≤Klog ν←(uk) +Wk}.The cumulative distribution function and the probability density function ofthe Gumbel distribution Gumbel(0, 1) isF (x) = e−e−x, f(x) = e−(x+e−x).SoP(log ν←(ΓK) +WK ≤ x | ΓK) = e−ν←(ΓK)e−x ,andP(log ν←(uk) +Wk ≤ x | ΓK) =∫ 10e−ν←(ΓKu)e−xdu.Therefore,P(MK ≤ x | ΓK) =(∫ 10e−ν←(ΓKu)e−xdu)K−de−ν←(ΓK)e−x .101A.2. TRUNCATION ERROR BOUNDS FOR SELF-PRODUCT MEASURESDenote Q(u, t) = e−ν←(u)e−t , andP(MK ≤ x | ΓK) =(∫ 10Q(ΓKu, x)du)K−dQ(ΓK , x).Conditioned on ΓK , the tail ΓK+1,ΓK+2, · · · is a unit rate homogeneous Pois-son process on [ΓK ,∞) that is independent of Γ1, · · · ,ΓK−1. So conditionedon ΓK ,MK+d= supk≥1log ν←(ΓK + Γ′k) +Wk,where Γ′k is unit rate of homogeneous Poisson process on R+. Since Γ′k is aPoisson process, log ν←(ΓK + Γ′k) +Wk is also a Poisson process with therate measure (∫ ∞0e−(t−log ν←(ΓK+γ))−e−(t−log ν←(ΓK+γ))dγ)dt.P(MK+ ≤ x | ΓK) is the probability that no atom of the Poisson processis greater than x. For a Poisson process with rate measure ν(dt), thisprobability is e−∫∞x ν(dt),P(MK+ ≤ x | ΓK) = e− ∫∞x (∫∞0 e−(t−log ν←(ΓK+γ))−e−(t−log ν←(Γ+γ))dγ)dt= e∫∞0(e−ν←(ΓK+γ)e−x−1)dγ= e∫∞0 Q(ΓK+γ, x)−1 dγ ,where the second equation comes from the fact that the inner integrand isthe probability density function of a Gumbel distribution. Therefore,P(MK < MK+ | ΓK)=∫ ∞−∞P(MK ≤ x | ΓK) ddxP(MK+ ≤ x | ΓK)dx=∫ ∞−∞Q(ΓK , x)(∫ 10Q(ΓKu, x)du)K−d(ddxe∫∞0 Q(ΓK+γ, x)−1 dγ)dx.For the categorical variable J with class probabilities given by θ′j/(1−P1,d−1),j ≥ d, it holds that P(d ≤ J ≤ K) ↑ 1 as K → ∞. By the monotoneconvergence theoremBK = P(argmaxj≥d(log θj +Wj) > K)= 1− E[Pd,K1− P1,d]→ 0.102A.3. ERROR BOUNDS FOR EDGE-EXCHANGEABLE NETWORKSA.3 Error bounds for edge-exchangeable networksA.3.1 Rejection representationWe first derive the specific form of BK in Corollary 3.4.3 for the rejectionrepresentation from Eq. (2.12), whereτ(U,Γ) = νˆ←(Γ)1[dνdνˆ(νˆ←(Γ)) ≥ U],and U ∼ Unif[0, 1]. So in Corollary 3.4.3,BK,1=E[∫R4+− log g0 (τ(u1, γK,1)τ(u2, γK,2)) dγ1dγ2gu(u1)du1gu(u2)du2]=E[∫ ∞ΓK∫ ∞ΓK∫ 10∫ 10− log g0 (τ(u1, γ1)τ(u2, γ2)) du1du2dγ1dγ2],where γK,1 := ΓK + γ1 and γK,2 := ΓK + γ2. Since g0(0) = 1, log g0(0) = 0,we can take the indicator in τ out of the function log g0(τ(u1, γ1)τ(u2, γ2))to obtain∫ 10∫ 10log g0 (τ(u1, γ1)τ(u2, γ2)) du1du2=∫ 10∫ 10log g0 (νˆ←(γ1)νˆ←(γ2))1 [rγ1 ≥ u1]1 [rγ2 ≥ u2] du1du2= log g0 (νˆ←(γ1)νˆ←(γ2))dνdνˆ(νˆ←(γ1))dνdνˆ(νˆ←(γ2)),where rγ1 :=dνdνˆ (νˆ←(γ1)) and rγ2 :=dνdνˆ (νˆ←(γ2)). Transforming variables viax1 = νˆ←(γ1) and x2 = νˆ←(γ2), and noting that νˆ←(ΓK) ≥ x if and only ifΓK ≤ νˆ[x,∞), we have:103A.3. ERROR BOUNDS FOR EDGE-EXCHANGEABLE NETWORKSBK,1=E[∫ νˆ←(ΓK)0∫ νˆ←(ΓK)0− log g0(x1x2)dνdνˆ(x1)dνdνˆ(x2)νˆ(dx1)νˆ(dx2)]=E[∫R2+− log g0(x1x2)1[x1 ≤ νˆ←(ΓK)]1[x2 ≤ νˆ←(ΓK)]ν(dx1)ν(dx2)]=∫R2+− log g0(x1x2)E [1[ΓK ≤ νˆ[max{x1, x2},∞)]] ν(dx1)ν(dx2)=∫R2+− log g0(x1x2)FK (νˆ[max{x1, x2},∞)) ν(dx1)ν(dx2), (A.2)where FK is the cumulative distribution function of ΓK . Using the samevariable transformation again,12(K − 1)BK,2=E[∫ ΓK0∫R3+− 1ΓKlog g0 (τ(u1, γ1)τ(u2, γK,2)) g(u1)du1g(u2)du2dγ1dγ2]=E[∫ ∞νˆ←(ΓK)∫ νˆ←(ΓK)0− 1ΓKlog g0(x1x2)dνdνˆ(x1)dνdνˆ(x2)νˆ(dx1)νˆ(dx2)]=∫R2+−E[1ΓK1[νˆ[x1,∞) ≤ ΓK ≤ νˆ[x2,∞)]]log g0(x1x2)ν(dx1)ν(dx2),where γK,2 = γK + γ2.BK,3 = 2E[∫R3+− log g0 (τ(u1,ΓK)τ(u2, γ2 + ΓK)) g(u1)du1g(u2)du2dγ2]= 2E[∫ ∞ΓK− log g0 (νˆ←(ΓK)νˆ←(γ2)) dνdνˆ(νˆ←(ΓK))dνdνˆ(νˆ←(γ2))dγ2]= 2E[∫ νˆ←(ΓK)0− log g0 (νˆ←(ΓK)x) dνdνˆ(νˆ←(ΓK))dνdνˆ(x)νˆ(dx)]= 2E[∫R+− log g0 (νˆ←(ΓK)x) dνdνˆ(νˆ←(ΓK))1[ΓK ≤ νˆ[x,∞)]ν(dx)].104A.3. ERROR BOUNDS FOR EDGE-EXCHANGEABLE NETWORKSA.3.2 Beta-independent Bernoulli process networkDense network When α = 0,ν(dθ) = γλθ−1(1− θ)λ−1dθ, νˆ(dθ) = γλθ−1dθ,anddνdνˆ= (1− θ)λ−1, νˆ ([x, 1]) = −γ′ log x, νˆ←(u) = e−u/γ′ .Substituting ν(dx), νˆ(dx) and h0(θ) into Eq. (A.2) and noting that theintegrand is symmetric around the line x1 = x2,BK,1 = 2∫ 10∫ 1x1− log(1− x1x2)FK (νˆ[x2, 1]) γλx−12 (1− x2)λ−1dx2ν(dx1).Next, note that 0 ≤ FK (νˆ[x2, 1]) ≤ FK (νˆ[x1, 1]) when x2 ≥ x1 and0 ≥ log(1− x1x2) ≥ −(11− x1x2 − 1)= −(x1x21− x1x2)≥ −(x1x21− x1).SoBK,1 ≤ 2γλ∫ 10FK (νˆ([x1, 1]))∫ 1x1x1x21− x1x−12 (1− x2)λ−1dx2ν(dx1)= 2γ2λ∫ 10(1− x)2(λ−1)FK (νˆ([x, 1])) dx.For any a ∈ (0, 1), dividing the integral into two parts and bounding eachpart separately,BK,1 ≤ 2γ2λ[∫ a0(1− x)2(λ−1)dx+∫ 1aFK(−γ′ log x) (1− x)2(λ−1)dx] .Assume for the moment that λ ̸= 0.5. Use the fact that FK(t) ≤ (3t/K)Kand note also that when x ∈ [a, 1], − log x ≤ [− log a/(1− a)] (1− x),12γ2λBK,1≤1− (1− a)2λ−12λ− 1 +∫ 1a(3γ′ [− log a/(1− a)] (1− x)K)K(1− x)2λ−2dx=1− (1− a)2λ−12λ− 1 +(−3γ′ log aK)K 1K + 2λ− 1(1− a)2λ−1.105A.3. ERROR BOUNDS FOR EDGE-EXCHANGEABLE NETWORKSIt can be seen that there exists a constant value c > 1 such that 3γ′ log c = 1/c.Setting a = c−K and using the first order Taylor’s expansion to approximatethe first term, it can be seen that1− (1− a)2λ−12λ− 1 ∼ c−K ,(−3γ′ log aK)K 1K + 2λ− 1(1− a)2λ−1 ∼ 1Kc−K .Therefore, there exists K0 such that when K > K0,BK,1 ≤ 4γ2λc−K .If λ = 0.5,BK,1 ≤ 2γ2λ∫ 10(1− x)−1FK (νˆ([x, 1])) dx≤ 2γ2λ[− log(1− a) +(3γ′K)K 1K(− log a)K],which can be bounded similarly by choosing the same c and setting a = c−K .Therefore we can find a constant K0 and for K > K0,BK,1 ≤ 4γ2λc−K .Next, the term BK,2 is bounded via12BK,2=∫ 10∫ 1x1− log(1− x1x2) [FK−1 (νˆ([x1, 1]))− FK−1 (νˆ([x2, 1]))] ν(dx2)ν(dx1)≤∫ 10∫ 1x1− log(1− x1x2)FK−1 (νˆ[x1, 1]) ν(dx2)ν(dx1).This has exactly the same form as BK,1, except that the CDF of ΓK isreplaced with that of ΓK−1. Therefore, it can be shown that for large K,BK,2 ≤ 4γ2λc−(K−1).Finally, BK,3 may be expressed as12BK,3=E[∫ 10− log (1− νˆ←(ΓK)x) dνdνˆ(νˆ←(ΓK))1 [x ≤ νˆ←(ΓK)]]ν(dx)=γλE[∫ νˆ←(Γk)0− log (1− νˆ←(ΓK)x) (1− νˆ←(ΓK))λ−1 x−1(1− x)λ−1dx].106A.3. ERROR BOUNDS FOR EDGE-EXCHANGEABLE NETWORKSSince log (1− νˆ←(ΓK)x) ≥ −νˆ←(ΓK)x/ (1− νˆ←(ΓK)),BK,3 ≤ 2γλE[(1− νˆ←(ΓK))λ−2 νˆ←(ΓK)∫ νˆ←(ΓK)0(1− x)λ−1dx]≥ 2γE[(1− νˆ←(ΓK))λ−2 νˆ←(ΓK)]= 2γ∫ ∞0(1− e−x/γ′)λ−2e−x/γ′ xK−1Γ(K)e−xdx.We split the analysis of this term into two cases. In the first case, assumingλ ≥ 2, we have thatBK,3 ≤ 2γ∫ ∞0xK−1Γ(K)e−x(1+γ′)/γ′dx = 2γ(γ′1 + γ′)K.On the other hand, if λ < 2, we bound the integral over [0, γ′] and over[γ′,∞) separately. Since 1− e−x ≥ x2 for x ∈ [0, 1],12γBK,3≤∫ γ′0(xγ′)2(λ−2) xK−1Γ(K)e−x(1+γ′)/γ′dx+ (1− e−1)λ−2∫ ∞γ′xK−1Γ(K)e−x(1+γ′)/γ′dx≤γ′2(2−λ)∫ γ′0xK−1+2(λ−2)Γ(K)dx+ (1− e−1)λ−2∫ ∞0xK−1Γ(K)e−x(1+γ′)/γ′dx=1Γ(K)γ′KK + 2(λ− 2) + (1− e−1)λ−2(γ′1 + γ′)K.As K →∞, the second term will dominate the first term, so when λ− 2 < 0,the following inequality holds for large K,BK,3 ≤ 4γ(1− e−1)λ−2(γ′1 + γ′)K.BK = BK,1 + BK,2 + BK,3 and as K → ∞, BK,3 will dominate BK,1 andBK,2. So there exists a K0 ∈ N such that for K > K0,BK ≤ 12γ(1− e−1)λ−2(γ′1 + γ′)K−→ 0.107A.3. ERROR BOUNDS FOR EDGE-EXCHANGEABLE NETWORKSSparse network When α > 0,ν(dθ) = γ′θ−1−α(1− θ)λ+α−1dθ, νˆ(dθ) = γ′θ−1−αdθ,anddνdνˆ= (1− θ)λ+α−1, νˆ ([x, 1]) = γ′α−1(x−α − 1),νˆ←(u) =(1 +αuγ′)−1/α.Similar to the case when α = 0,BK,1=2∫ 10∫ 1x1− log(1− x1x2)FK (νˆ([x2, 1])) γ′x−1−α2 (1− x2)λ+α−1dx2ν(dx1).Since log(1− x1x2) ≥ −x1x2/(1− x1),BK,1 ≤ 2γ′2∫ 10FK(γ′α−1(x−α − 1))x−2α(1− x)2[(λ+α)−1]dx≤ 2γ′2∫ 10(γ′α−1x−α(1− xα))KΓ(K + 1)x−2α(1− x)2[(λ+α)−1]dx≤ 2γ′2∫ 10(γ′α−1x−α)KΓ(K + 1)x−2α(1− x)2[(λ+α)−1]dx≤ 2α2∫ 10(γ′α−1x−α)K+2Γ(K + 1)dx.Denoting t = γ′α−1x−α, thenBK,1 ≤ 2α(γ′α−1)1/α∫ ∞γ′α−1tK+1−1/αΓ(K + 1)dt≤ 2α(γ′α−1)1/α eγ′α−1 Γ(K + 2− 1/α)Γ(K + 1)∫ ∞γ′α−1tK+1−1/αΓ(K + 2− 1/α)e−t dt≤ 2α(γ′α−1)1/α eγ′α−1 Γ(K + 2− 1/α)Γ(K + 1).By Stirling’s formula,BK,1 ≤ 2α(γ′α−1)1/α eγ′α−1 Kα−1α .108A.3. ERROR BOUNDS FOR EDGE-EXCHANGEABLE NETWORKSNow we consider the error bound for BK,2. As we have shown in the examplewhen α = 0, here we can obtain thatBK,2 ≤ 2α(γ′α−1)1/α eγ′α−1 (K − 1)α−1α .Similar to BK,1 and BK,2,BK,3=2E[∫ 10− log (1− νˆ←(ΓK)x) dνdνˆ(νˆ←(ΓK))1 (ΓK ≤ νˆ←[x, 1]) ν(dx)]≤2γ′E[(1− νˆ←(ΓK))λ+α−1∫ νˆ←(ΓK)0νˆ←(ΓK)x1− x x−1−α(1− x)λ+α−1dx]=2γ′E[(1− νˆ←(ΓK))λ+α−1 νˆ←(ΓK)∫ νˆ←(ΓK)0x−α(1− x)λ+α−2dx].We again split our analysis into two cases. First, suppose that λ+ α− 2 ≥ 0.ThenBK,3≤ 2γ′1− αE[(1− νˆ←(ΓK))λ+α−1 νˆ←(ΓK)2−α]=2γ′1− α∫ ∞0[1−(1 +αxγ′)−1/α]λ+α−1(1 +αxγ′)−(2−α)/α xK−1Γ(K)e−xdx.Since α < 1, (1 +αxγ′)−(2−α)/α≤(αxγ′)−(2−α)/α,and soBK,3 ≤ 2γ′1− α∫ ∞0(αγ′)−(2−α)/α xK−1−(2−α)/αΓ(K)e−xdx≤ 2γ′1− α(αγ′)−(2−α)/α Γ(K − (2− α)/α)Γ(K)∼ 2γ′1− α(αγ′)−(2−α)/αK−(2−α)/α,where the last equation is obtained from Stirling’s formula.109A.3. ERROR BOUNDS FOR EDGE-EXCHANGEABLE NETWORKSOn the other hand, if λ+ α− 2 < 0,BK,3≤2γ′E[(1− νˆ←(ΓK))λ+α−1 νˆ←(ΓK)∫ νˆ←(ΓK)0x−α(1− x)λ+α−2dx]≤ 2γ′1− αE[(1− νˆ←(ΓK))2(λ+α)−3 νˆ←(ΓK)2−α]=2γ′1− α∫ ∞0[1−(1 +αxγ′)−1/α]2(λ+α)−3(1 +αxγ′)−(2−α)/α xK−1Γ(K)e−xdx.If 2(λ+α)− 3 ≥ 0, we get the same result as in the case when λ+α− 2 ≥ 0.When 2(λ+α)−3 ≤ 0, note that we can find an x0 such that when x ∈ [0, x0],1−(1 +αxγ′)−1/α≥ x2.SoBK,3≤ 2γ′1− α{∫ x00x4(λ+α)−6(1 +αxγ′)−(2−α)/α xK−1Γ(K)e−xdx+[1−(1 +αx0γ′)−1/α]2(λ+α)−3 ∫ ∞x0(1 +αxγ′)−(2−α)/α xK−1Γ(K)e−xdx≤ γ′1− α(αγ′)−(2−α)/α{∫ x00xK−1+4(λ+α)−6−(2−α)/αΓ(K)e−xdx+[1−(1 +αx0γ′)−1/α]2(λ+α)−3 Γ(K − (2− α)/α)Γ(K) .Because here we assume that 2(λ+α)−3 < 0, the second term will dominate.So in this case we obtain that for large K,BK,3 ≤ 4γ′1− α(αγ′)−(2−α)/α [1−(1 +αx0γ′)−1/α]2(λ+α)−3K−(2−α)/α.Asymptotically, BK,2 will dominate BK,1 and BK,3, so there exists K0 ∈ Nsuch that when K > K0,BK ≤ 6α(γ′α−1)1/α eγ′α−1 (K − 1)α−1α −→ 0.110A.3. ERROR BOUNDS FOR EDGE-EXCHANGEABLE NETWORKSA.3.3 Gamma-independent Poisson networkDense network When α = 0,ν(dθ) = γλθ−1e−λθdθ, νˆ(dθ) = γλθ−1(1 + λθ)−1dθ.In this case,dνdνˆ= (1 + λθ)e−λθ, νˆ[x,∞) = γλ log(1 + (λx)−1),νˆ←(u) =1λ(e(γλ)−1u − 1) .For Poisson distribution, h0(θ) = e−θ, soBK,1 =∫R2+x1x2FK (νˆ[max{x1, x2},∞)) ν(dx1)ν(dx2).Note that the integrand is symmetric around the line x1 = x2, so we only needto compute the integral above the line x1 = x2. In this region, e−λx2 ≤ e−λx1and 0 ≤ FK (νˆ[x2,∞)) ≤ FK (νˆ[x1,∞)) ≤ 1. SoBK,1 = 2∫R+x1∫ ∞x1x2FK (νˆ[x2,∞))λγx−12 e−λx2dx2ν(dx1)≤ 2∫R+x1FK (νˆ[x1,∞)) γe−λx1ν(dx1)= 2γ2λ∫R+FK (νˆ[x,∞)) e−2λxdx.For any a > 0, we divide the integral into two parts and bound eachpart separately. We denote b = log(1 + (λa)−1) and use the fact that∫ a0 e−2λxdx ≤ a and FK(t) ≤ tK/K! ≤ (3t/K)K . SoBK,1 ≤ 2γ2λ[∫ a0e−2λxdx+ FK (νˆ[a,∞))∫ ∞ae−2λxdx]≤ 2γ2λ[a+12λFK(γλ log(1 + (λa)−1))]≤ 2γ2λ[(λ(eb − 1))−1+1λ(3γλbK)K].Setting two terms equal and use the fact that (eb − 1)−1 ≈ e−b when b islarge, we get b = KW0((3γλ)−1) and W0 is defined byW0(y) = x ⇐⇒ xex = y.111A.3. ERROR BOUNDS FOR EDGE-EXCHANGEABLE NETWORKSTherefore,BK,1 ≤ 4γ2eKW0((3γλ)−1) − 1 ∼ 4γ2e−KW0((3γλ)−1).Similarly,12(K − 1)BK,2=∫R2+x1x2E[1ΓK1 (νˆ[x2,∞) ≤ ΓK ≤ νˆ[x1,∞))]ν(dx1)ν(dx2).Note thatE[1ΓK1 (νˆ[x2,∞) ≤ ΓK ≤ νˆ[x1,∞))]=1K − 1 [FK−1 (νˆ[x1,∞))− FK−1 (νˆ[x2,∞))] .Keeping only the positive part,BK,2 ≤ 2∫R+x1FK−1 (νˆ[x1,∞))∫ ∞x1γλe−λx2dx2ν(dx1)= 2γ2λ∫R+FK−1 (νˆ[x,∞)) e−2λxdx.This has the same form as BK,1, soBK,2 ≤ 4γ2e(K−1)W0((3γλ)−1) − 1 ∼ 4γ2e−(K−1)W0((3γλ)−1).Next,BK,3=2E[∫R+xνˆ←(ΓK)dνdνˆ(νˆ←(ΓK))1 [ΓK ≤ νˆ[x,∞)]]ν(dx)=2E[∫ νˆ←(ΓK)0xνˆ←(ΓK) (1 + λνˆ←(ΓK)) e−λνˆ←(ΓK)γλx−1e−λxdx]=2γE[νˆ←(ΓK) (1 + λνˆ←(ΓK)) e−λνˆ←(ΓK)(1− e−λνˆ←(ΓK))]≤2γλE[νˆ←(ΓK)2 (1 + λνˆ←(ΓK)) e−λνˆ←(ΓK)].112A.3. ERROR BOUNDS FOR EDGE-EXCHANGEABLE NETWORKSNote that (1 + x)e−x ≤ 1, soBK,3 ≤ 2λγE[νˆ←(ΓK)2]≤ 2γλE[e−(γλ)−1ΓK]= 2γλ(γλ1 + γλ)K−1.Since BK,3 will dominate BK,1 and BK,2 asymptotically, there exists K0 ∈ Nsuch that for K > K0,BK ≤ 3BK,3 ≤ 6γλ(γλ1 + γλ)K−1.Sparse network When α > 0,ν(dθ) = γλ1−αΓ(1− α)θ−α−1e−λθdθ, νˆ(dθ) = γλ1−αΓ(1− α)θ−α−1dθ,anddνdνˆ= e−λθ, νˆ[x,∞) = γ′x−α, νˆ←(u) = (γ′u−1)1/α, γ′ = γ λ1−ααΓ(1− α) .Similar to the example when α = 0,BK,1 =∫R2+x1x2FK (νˆ[max{x1, x2},∞)) ν(dx1)ν(dx2)= 2∫R+x1∫ ∞x1FK (νˆ[x2,∞)) γ λ1−αΓ(1− α)x−α2 e−λx2dx2ν(dx1)≤ 2γ∫R+x1FK (νˆ[x1,∞))∫ ∞x1λ1−αΓ(1− α)x−α2 e−λx2dx2ν(dx1).Note that the integrand with respect to x2 is the density function of thegamma distribution with shape α and rate λ, so the integral is less than 1.We partition the outer integral into two parts and bound them separately,BK,1 ≤ 2γ∫R+FK (νˆ[x,∞)) γ λ1−αΓ(1− α)x−αe−λxdx≤ 2γ2 λ1−αΓ(1− α)[∫ a0x−αdx+ FK (νˆ[a,∞))∫ ∞ax−αeλxdx]≤ 2γ2 1Γ(1− α)[λ1−α1− αa1−α + Γ(1− α)(3γ′a−αK)K].113A.3. ERROR BOUNDS FOR EDGE-EXCHANGEABLE NETWORKSBy setting the two terms in the brackets equal, we geta =[(1− α)Γ(1− α)λ1−α] 1(K−1)α+1(3γ′K) K(K−1)α+1∼(3γ′K) 1α.SoBK,1 ≤ 4γ2λ1−α(1− α)Γ(1− α)(3γ′K) 1−αα.Similar to the last example where α = 0, hereBK,2=2(K − 1)∫R2+x1x2E[1ΓK1 (νˆ[x2,∞) ≤ ΓK ≤ νˆ[x1,∞))]ν(dx1)ν(dx2)≤2∫R+x1FK−1 (νˆ[x1,∞))∫ ∞x1x2ν(dx2)ν(dx1)≤2γ∫R+x1FK−1 (νˆ[x1,∞)) ν(dx1).This has the same form as BK,1, and thereforeBK,2 ≤ 4γ2λ1−α(1− α)Γ(1− α)(3γ′K − 1) 1−αα.Finally, since both e−λx < 1 and e−λνˆ←(ΓK) < 1,BK,3=2E[∫R+νˆ←(ΓK)xe−λνˆ←(ΓK)1 [x ≤ νˆ←(ΓK)] γ λ1−αΓ(1− α)x−1−αe−λxdx]≤2γ λ1−αΓ(1− α)E[νˆ←(ΓK)∫ νˆ←(ΓK)0x−αdx]≤2γ λ1−α(1− α)Γ(1− α)E[νˆ←(ΓK)2−α]=2γλ1−α(1− α)Γ(1− α)(γ′) 2−ααΓ(K − 2−αα )Γ(K).By Stirling’s formula,Γ(K − (2− α)/α)Γ(K)∼√2π(K − 2−αα )(K−(2−α)/αe)K−(2−α)/α√2πK(Ke)K ∼ K− 2−αα .114A.4. ERRORS OF TRUNCATED BNP MODEL WITH LATENT FEATURESSoBK,3 ≤ 2γ λ1−α(1− α)Γ(1− α)(γ′) 2−αα K−2−αα .BK,2 dominates BK,1 and BK,3 asymptotically, so there exists K0 ∈ N suchthat for K > K0,BK ≤ 12γ2λ1−α(1− α)Γ(1− α)(3γ′K − 1) 1−αα.A.4 Errors of truncated BNP model with latentfeaturesProof of Corollary 5.3.2. Since ru = 1,DTV(Πˆµ,Πµ)≤ 12EΠˆµ[max{| ruEΠˆµ [rℓ]− 1|, | rℓEΠˆµ [ru]− 1|}]=12EΠˆµ[max{ 1EΠˆµ [rℓ]− 1, 1− rℓ}]=12EΠˆµ[max{ 1EΠˆµ [rℓ]− 1, ϵ1A + 1Ac}]By definitionEΠˆµ [rℓ] = (1− ϵ)(1− η),1EΠˆµ [rℓ]− 1 = ϵ+ η − ηϵ1− ϵ− η + ηϵ115A.4. ERRORS OF TRUNCATED BNP MODEL WITH LATENT FEATURESIf we further assume that η and ϵ are both small such that Eˆ[rℓ] > 1/2, thenDTV(Πˆµ,Πµ)≤12EΠˆµ[(ϵ1A + 1Ac)1[ϵ1A + 1Ac >1EΠˆµ [rℓ]− 1]]+12EΠˆµ[(1EΠˆµ [rℓ]− 1)1[ϵ1A + 1Ac ≤ 1EΠˆµ [rℓ]− 1]]=12EΠˆµ [(ϵ1A + 1Ac)1Ac ] +12EΠˆµ[(1EΠˆµ [rℓ]− 1)1A]=12η +12(1− η)1− (1− ϵ)(1− η)(1− ϵ)(1− η)=12η +121− (1− ϵ)(1− η)(1− ϵ)=12(ϵ1− ϵ + 2η)Proof of Lemma 5.3.3. In this section derive the upper and the lower boundon the ratiof(Y1:N | ϱ, θ1:K , ψ1:K , z1:N,1:K)fˆ(Y1:N | θ1:K , ψ1:K , z1:N,1:K).We first derive the upper bound. Let rn = yn−∑Kk=1 znkψk, and assumef˜ ≤ 0. We denoteP˜n0 := P( ∞∑k=K+1zn,k = 0 | θ1:∞)and P˜n+ := P( ∞∑k=K+1zn,k > 0 | θ1:∞).116A.4. ERRORS OF TRUNCATED BNP MODEL WITH LATENT FEATURESThen f(Y1:N | ϱ, θ1:K , ψ1:K , z1:N,1:K) can be written as:E[e∑Nn=1 f˜(yn−∑∞k=1 zn,kψk) | ϱ, θ1:K , ψ1:K , z1:N,1:K]=E[e∑Nn=1 f˜(rn−∑∞k=K+1 zn,kψk) | ϱ, θ1:K , ψ1:K , z1:N,1:K]=E[E[e∑Nn=1 f˜(rn−∑∞k=K+1 zn,kψk) | θ1:∞, ψ∞, z1:N,1:K]| ϱ, θ1:K , ψ1:K , z1:N,1:K]=E[N∏n=1E[ef˜(rn−∑∞k=K+1 zn,kψk) | θ1:∞, ψ1:∞, z1:N,1:K]| ϱ, θ1:K , ψ1:K , z1:N,1:K]≤E[N∏n=1(ef˜(rn)P˜n0 + P˜n+)| ϱ, θ1:K , ψ1:K , z1:N,1:K]≤E[N∏n=1(ef˜(rn) + P˜n+)| ϱ, θ1:K , ψ1:K , z1:N,1:K]=e∑n f˜(rn)E[N∏n=1(1 + e−f˜(rn)P˜n+)| ϱ, θ1:K , ψ1:K , z1:N,1:K]=e∑n f˜(rn)E[N∏n=1(1 + e−f˜(rn)∑k>Kθk)| ϱ, θ1:K , ψ1:K , z1:N,1:K]=fˆ(Y1:N | ·) · E[N∏n=1(1 + e−f˜(rn)∑k>Kθk)| ϱ, θ1:K , ψ1:K , z1:N,1:K].We denote Sθ =∑∞k=K+1 θk and consider the conditional expectation ofthe product∏Nn=1(1 + e−f˜(rn)Sθ). According to the binomial theorem,E[N∏n=1(1 + e−f˜(rn)∑k>Kθk)]=E[ ∑In⊂1:Ne−∑n∈In f˜(rn)S|In|θ]=∑In⊂1:Ne−∑n∈In f˜(rn)E[S|In|θ ]We denote f(x) = ex as the exponential function, and g(t) =∫(etθ −1)νθK (dθ). By definition, the nth derivative of g(t), evaluated at t = 0, i.e.,mi = g(i)(t)|t=0 =∫θiνθKd(θ). In deriving the following upper bound, weassume that we can find an upper bound for mi such that mi ≤ u0 · ui1 fori = 1, . . . , N , where both u0 and u1 are functions of θK .We apply the Campbell’s theorem and the Faa´ di Bruno’s formula thatgeneralizes the chain rule and gives an expression for the nth derivative117A.4. ERRORS OF TRUNCATED BNP MODEL WITH LATENT FEATURESdnf(g(t))/dnt:E[Snθ | θK ] =dnE[etSθ ]dtn∣∣∣∣t=0=dn exp{∫ (etθ − 1)νθK (dθ)}dtn∣∣∣∣t=0=Bn(g′, g′′, . . . , g(n)) =n∑i=1Bn,i(g′, g′′, . . . , g(n−i+1))=n!∑1j1+...+njn=nn∏i=1mjii(i!)jiji!≤n!∑1j1+...+njn=nn∏i=1(u0 · ui1)ji(i!)jiji!≤ (u0 · u1)nn!∑1j1+...+njn=nn∏i=11(i!)jiji!=UnθKBn(1, . . . , 1) = UnθKn∑i=1Bn,i(1, . . . , 1)≤UnθK(0.792 · nlog(n+ 1))n≤ UnθK(0.792 ·Nlog(N + 1))n,where UθK := u0 × u1, Bn,i(x1, . . . , xn−k+1) are the incomplete exponentialBell polynomials and Bn(x1, . . . , xn) is the nth complete exponential Bell’spolynomial. Bn(1, . . . , 1) is called the nth Bell number.Place back E[Siθ|θK ], we getE[N∏n=1(1 + e−f˜(rn)∑k>Kθk)| ϱ, θ1:K]≤∑In⊂1:Ne−∑n∈IN f˜(rn)( 4N5 log(N + 1)UθK)|In|=N∏n=1(1 + e−f˜(rn)4N5 log(N + 1)UθK)118A.4. ERRORS OF TRUNCATED BNP MODEL WITH LATENT FEATURESWe then derive the lower bound of the ratio:E[N∏n=1E[ef˜(rn−∑∞k=K+1 zn,kψk) | θ1:∞, ψ1:∞, z1:N,1:K]| ϱ, θ1:K , ψ1:K , z1:N,1:K]≥E[N∏n=1(ef˜(rn)P( ∞∑k=K+1zn,k = 0 | θ1:∞))| ϱ, θ1:K , ψ1:K , z1:N,1:K]=e∑n f˜(rn) P(Z1:N,K+ = 0 | θ1:K)We follow the approach in Li and Campbell (2021) to derive the posteriorapproximation error of the truncated model. (Campbell et al., 2019, Theorem4.2) derived the truncated error for general series representationP(Z[N ],[K+] = 0 | θK) ≥ exp(−N ·B(θK))≥ exp{−N∫ ∞0(1− π(τ(V, u+GK)))du}where GK ∼ Gam(K, 1). In particular, for the inverse Le´vy representation,B(θK) ≤∫ θK0(1− g(0 | θ))ν(dθ).Therefore, with the inverse Le´vy representation,P(Z[N ],[K+] = 0 | θK) ≥ exp[−N ·∫ θK0(1− g(0 | θ))ν(dθ)]=: exp(−N · LθK ).119BibliographyAldous, D. J. (1981). Representations for partially exchangeable arrays ofrandom variables. Journal of Multivariate Analysis, 11(4):581–598.Antoniak, C. E. (1974). Mixtures of Dirichlet processes with applicationsto Bayesian nonparametric problems. The annals of statistics, pages1152–1174.Argiento, R., Bianchini, I., and Guglielmi, A. (2016). A blocked Gibbssampler for NGG-mixture models via a priori truncation. Statistics andComputing, 26(3):641–661.Argiento, R., Cremaschi, A., and Vannucci, M. (2019). Hierarchical normal-ized completely random measures to cluster grouped data. Journal of theAmerican Statistical Association.Arlot, S. and Celisse, A. (2010). A survey of cross-validation procedures formodel selection. Statistics Surveys, 4(none):40 – 79.Blackwell, D. and MacQueen, J. B. (1973). Ferguson distributions via po´lyaurn schemes. The annals of statistics, 1(2):353–355.Blei, D. M. and Jordan, M. I. (2006). Variational inference for Dirichletprocess mixtures. Bayesian Analysis, 1(1):121 – 143.Blei, D. M., Jordan, M. I., and Ng, A. Y. (2003a). Hierarchical Bayesianmodels for applications in information retrieval. Bayesian statistics, 7(2):5–43.Blei, D. M., Ng, A. Y., and Jordan, M. I. (2003b). Latent Dirichlet allocation.Journal of machine Learning research, 3(Jan):993–1022.Bondesson, L. (1982). On simulation from infinitely divisible distributions.Advances in Applied Probability, 14(4):855–869.Borgs, C., Chayes, J. T., Cohn, H., and Holden, N. (2018). Sparse exchange-able graphs and their limits via graphon processes. Journal of MachineLearning Research, 18(210):1–71.120BibliographyBrix, A. (1999). Generalized gamma measures and shot-noise Cox processes.Advances in Applied Probability, 31(4):929–953.Broderick, T., Jordan, M. I., and Pitman, J. (2012). Beta Processes, Stick-Breaking and Power Laws. Bayesian Analysis, 7(2):439 – 476.Broderick, T., Jordan, M. I., and Pitman, J. (2013). Cluster and FeatureModeling from Combinatorial Stochastic Processes. Statistical Science,28(3):289 – 312.Broderick, T., Mackey, L., Paisley, J., and Jordan, M. I. (2014). Combinato-rial clustering and the beta negative binomial process. IEEE Transactionson Pattern Analysis and Machine Intelligence, 37(2):290–306.Broderick, T., Wilson, A. C., and Jordan, M. I. (2018). Posteriors, conjugacy,and exponential families for completely random measures. Bernoulli,24(4B):3181 – 3221.Cai, D., Campbell, T., and Broderick, T. (2016). Edge-exchangeable graphsand sparsity. Advances in Neural Information Processing Systems, 29.Camerlenghi, F., Lijoi, A., Orbanz, P., and Pru¨nster, I. (2019). Distributiontheory for hierarchical processes. The Annals of Statistics, 47(1):67 – 92.Campbell, T., Cai, D., and Broderick, T. (2018). Exchangeable trait alloca-tions. Electronic Journal of Statistics, 12(2):2290 – 2322.Campbell, T., Huggins, J., How, J., and Broderick, T. (2019). Truncatedrandom measures. Bernoulli, 25(2):1256–1288.Caron, F. and Fox, E. (2017). Sparse graphs using exchangeable randommeasures. Journal of the Royal Statistical Society: Series B (StatisticalMethodology), 79(5):1295–1366.Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betan-court, M., Brubaker, M. A., Guo, J., Li, P., and Riddell, A. (2017). Stan:A probabilistic programming language. Journal of statistical software, 76.Crane, H. and Dempsey, W. (2018). Edge exchangeable models for interactionnetworks. Journal of the American Statistical Association, 113(523):1311–1326. PMID: 30467447.De Finetti, B. (1929). Funzione caratteristica di un fenomeno aleatorio. InAtti del Congresso Internazionale dei Matematici: Bologna del 3 al 10 desettembre di 1928, pages 179–190.121BibliographyDoshi, F., Miller, K., Van Gael, J., and Teh, Y. W. (2009). Variationalinference for the Indian buffet process. In Artificial Intelligence andStatistics, pages 137–144.Doshi-Velez, F. and Ghahramani, Z. (2009). Accelerated sampling for theIndian buffet process. In Proceedings of the 26th annual internationalconference on machine learning, pages 273–280.Doshi-Velez, F., Mohamed, S., Ghahramani, Z., and Knowles, D. (2009).Large scale nonparametric Bayesian inference: Data parallelisation in theIndian buffet process. Advances in Neural Information Processing Systems,22.Dubey, K. A., Zhang, M., Xing, E., and Williamson, S. (2020). Distributed,partially collapsed MCMC for Bayesian nonparametrics. In Interna-tional Conference on Artificial Intelligence and Statistics, pages 3685–3695.PMLR.Efron, B. and Thisted, R. (1976). Estimating the number of unsen species:How many words did Shakespeare know? Biometrika, 63(3):435–447.Ferguson, T. S. (1973). A Bayesian analysis of some nonparametric problems.The annals of statistics, pages 209–230.Ferguson, T. S. and Klass, M. J. (1972). A representation of independentincrement processes without Gaussian components. The Annals of Mathe-matical Statistics, 43(5):1634–1643.Fisher, R. A., Corbet, A. S., and Williams, C. B. (1943). The relationbetween the number of species and the number of individuals in a randomsample of an animal population. Journal of Animal Ecology, 12(1):42–58.Gelman, A., Carlin, J., Stern, H., Dunson, D., Vehtari, A., and Rubin, D.(2013). Bayesian Data Analysis, Third Edition. Chapman & Hall/CRCTexts in Statistical Science. Taylor & Francis.Ghahramani, Z. (1994). Factorial learning and the EM algorithm. Advancesin neural information processing systems, 7.Gnedin, A., Hansen, B., and Pitman, J. (2007). Notes on the occupancyproblem with infinitely many boxes: general asymptotics and power laws.Goldwater, S., Johnson, M., and Griffiths, T. (2005). Interpolating be-tween types and tokens by estimating power-law generators. In Weiss,122BibliographyY., Scho¨lkopf, B., and Platt, J., editors, Advances in Neural InformationProcessing Systems, volume 18. MIT Press.Griffin, J. E. (2016). An adaptive truncation method for inference in Bayesiannonparametric models. Statistics and Computing, 26(1):423–441.Griffin, J. E. and Walker, S. G. (2011). Posterior simulation of normal-ized random measure mixtures. Journal of Computational and GraphicalStatistics, 20(1):241–259.Griffiths, T. L. and Ghahramani, Z. (2011). The Indian buffet process: Anintroduction and review. Journal of Machine Learning Research, 12(4).Herlau, T., Schmidt, M. N., and Mørup, M. (2016). Completely randommeasures for modelling block-structured sparse networks. Advances inNeural Information Processing Systems, 29.Hjort, N. L. (1990). Nonparametric Bayes estimators based on beta processesin models for life history data. the Annals of Statistics, pages 1259–1294.Hoover, D. N. (1979). Relations on probability spaces and arrays of randomvariables. Preprint, Institute for Advanced Study, Princeton, NJ, 2.Hougaard, P. (1986). Survival models for heterogeneous populations derivedfrom stable distributions. Biometrika, 73(2):387–396.Ishwaran, H. and James, L. F. (2001). Gibbs sampling methods forstick-breaking priors. Journal of the American statistical Association,96(453):161–173.Ishwaran, H. and Zarepour, M. (2000). Markov chain Monte Carlo inapproximate Dirichlet and beta two-parameter process hierarchical models.Biometrika, 87(2):371–390.Ishwaran, H. and Zarepour, M. (2002). Exact and approximate sum repre-sentations for the Dirichlet process. The Canadian Journal of Statistics /La Revue Canadienne de Statistique, 30(2):269–283.James, L. F., Lijoi, A., and Pru¨nster, I. (2009). Posterior analysis fornormalized random measures with independent increments. ScandinavianJournal of Statistics, 36(1):76–97.Janson, S. (2018). On edge exchangeable random graphs. Journal of statisticalphysics, 173(3):448–484.123BibliographyKingman, J. (1967). Completely random measures. Pacific Journal ofMathematics, 21(1):59–78.Kingman, J. F. (1975). Random discrete distributions. Journal of the RoyalStatistical Society: Series B (Methodological), 37(1):1–15.Kingman, J. F. (1978). The representation of partition structures. Journalof the London Mathematical Society, 2(2):374–380.Kingman, J. F. C. (1992). Poisson processes, volume 3. Clarendon Press.Last, G. and Penrose, M. (2017). Lectures on the Poisson process, volume 7.Cambridge University Press.Lee, J., Miscouridou, X., and Caron, F. (2023). A unified construction forseries representations and finite approximations of completely randommeasures. Bernoulli, 29(3):2142–2166.Li, X. and Campbell, T. (2021). Truncated simulation and inference in edge-exchangeable networks. Electronic Journal of Statistics, 15(2):5117–5157.Lijoi, A., Mena, R. H., and Pru¨nster, I. (2005). Hierarchical mixture mod-eling with normalized inverse-Gaussian priors. Journal of the AmericanStatistical Association, 100(472):1278–1291.Lijoi, A., Mena, R. H., and Pru¨nster, I. (2007a). Bayesian nonparametricestimation of the probability of discovering new species. Biometrika,94(4):769–786.Lijoi, A., Mena, R. H., and Pru¨nster, I. (2007b). Controlling the reinforcementin Bayesian non-parametric mixture models. Journal of the Royal StatisticalSociety Series B: Statistical Methodology, 69(4):715–740.Lijoi, A., Nipoti, B., and Pru¨nster, I. (2014). Bayesian inference withdependent normalized completely random measures. Bernoulli, 20(3):1260– 1291.Lijoi, A., Pru¨nster, I., et al. (2010). Models beyond the Dirichlet process.Bayesian nonparametrics, 28(80):342.Muliere, P. and Tardella, L. (1998). Approximating distributions of randomfunctionals of Ferguson-Dirichlet priors. Canadian Journal of Statistics,26(2):283–297.124BibliographyNaik, C., Caron, F., Rousseau, J., Teh, Y. W., and Palla, K. (2022). Bayesiannonparametrics for sparse dynamic networks. In Joint European Conferenceon Machine Learning and Knowledge Discovery in Databases, pages 191–206. Springer.Neal, R. M. (2000). Markov chain sampling methods for Dirichlet processmixture models. Journal of computational and graphical statistics, 9(2):249–265.Neal, R. M. (2003). Slice sampling. The Annals of Statistics, 31(3):705 –767.Orbanz, P. and Roy, D. (2014). Bayesian models of graphs, arrays and otherexchangeable random structures. IEEE Transactions on Pattern Analysisand Machine Intelligence, 37(2):437–461.Paisley, J., Blei, D., and Jordan, M. (2012). Stick-breaking beta processesand the Poisson process. In Lawrence, N. D. and Girolami, M., editors,Proceedings of the Fifteenth International Conference on Artificial In-telligence and Statistics, volume 22 of Proceedings of Machine LearningResearch, pages 850–858, La Palma, Canary Islands. PMLR.Panzarasa, P., Opsahl, T., and Carley, K. M. (2009). Patterns and dynamicsof users’ behavior and interaction: Network analysis of an online community.Journal of the American Society for Information Science and Technology,60(5):911–932.Perman, M. (1990). Random discrete distributions derived from subordinators.University of California, Berkeley.Perman, M., Pitman, J., and Yor, M. (1992). Size-biased sampling of Poissonpoint processes and excursions. Probability Theory and Related Fields,92(1):21–39.Pitman, J. (1996). Some developments of the Blackwell-Macqueen urnscheme. Lecture Notes-Monograph Series, pages 245–267.Pitman, J. (2003). Poisson-kingman partitions. Lecture Notes-MonographSeries, pages 1–34.Pitman, J. (2006). Combinatorial stochastic processes: Ecole d’ete´ de proba-bilite´s de saint-flour xxxii-2002. Springer.125BibliographyPitman, J. and Yor, M. (1997). The two-parameter Poisson-Dirichlet dis-tribution derived from a stable subordinator. The Annals of Probability,pages 855–900.Regazzini, E., Lijoi, A., and Pru¨nster, I. (2003). Distributional results formeans of normalized random measures with independent increments. TheAnnals of Statistics, 31(2):560–585.Rosinski, J. (1990). On series representations of infinitely divisible randomvectors. The Annals of Probability, 18(1):405–430.Rosin´ski, J. (2001). Series representations of Le´vy processes from the per-spective of point processes. In Le´vy processes: Theory and applications,pages 401–415. Springer.Roweis, S. and Ghahramani, Z. (1999). A unifying review of linear Gaussianmodels. Neural computation, 11(2):305–345.Roy, D. M. (2014). The continuum-of-urns scheme, generalized betaand Indian buffet processes, and hierarchies thereof. arXiv preprintarXiv:1501.00208.Roychowdhury, A. and Kulis, B. (2015). Gamma processes, stick-breaking,and variational inference. In Artificial Intelligence and Statistics, pages800–808. PMLR.Schervish, M. J. (2012). Theory of statistics. Springer Science & BusinessMedia.Sethuraman, J. (1994). A constructive definition of Dirichlet priors. Statisticasinica, pages 639–650.Sudderth, E. and Jordan, M. (2008). Shared segmentation of natural scenesusing dependent Pitman-Yor processes. Advances in neural informationprocessing systems, 21.Teh, Y. W. (2006). A hierarchical Bayesian language model based on Pitman-Yor processes. In Proceedings of the 21st International Conference onComputational Linguistics and 44th Annual Meeting of the Association forComputational Linguistics, pages 985–992.Teh, Y. W. and Gorur, D. (2009). Indian buffet processes with power-lawbehavior. In Advances in Neural Information Processing Systems, pages1838–1846.126BibliographyTeh, Y. W., Gru¨r, D., and Ghahramani, Z. (2007). Stick-breaking construc-tion for the Indian buffet process. In Artificial intelligence and statistics,pages 556–563. PMLR.Teh, Y. W. and Jordan, M. I. (2010). Hierarchical Bayesian nonparametricmodels with applications. Bayesian nonparametrics, 1:158–207.Teh, Y. W., Jordan, M. I., Beal, M. J., and Blei, D. M. (2006). Hierarchi-cal Dirichlet processes. Journal of the American Statistical Association,101(476):1566–1581.Thibaux, R. and Jordan, M. I. (2007). Hierarchical beta processes andthe Indian buffet process. In Artificial Intelligence and Statistics, pages564–571. PMLR.Thibaux, R. J. (2008). Nonparametric Bayesian models for machine learning.University of California, Berkeley.Tipping, M. E. and Bishop, C. M. (1999). Probabilistic principal componentanalysis. Journal of the Royal Statistical Society Series B: StatisticalMethodology, 61(3):611–622.Titsias, M. (2007). The infinite gamma-Poisson feature model. Advances inNeural Information Processing Systems, 20.Todeschini, A., Miscouridou, X., and Caron, F. (2020). Exchangeable randommeasures for sparse and modular graphs with overlapping communities.Journal of the Royal Statistical Society Series B: Statistical Methodology,82(2):487–520.Tversky, A. (1972). Elimination by aspects: A theory of choice. Psychologicalreview, 79(4):281.Veitch, V. and Roy, D. (2015). The class of random graphs arising fromexchangeable random measures. arXiv:1512.03099.Walker, S. G. (2007). Sampling the Dirichlet mixture model with slices.Communications in Statistics—Simulation and Computation®, 36(1):45–54.Williamson, S. A. (2016). Nonparametric network models for link prediction.Journal of Machine Learning Research, 17(202):1–21.127BibliographyWood, F., Archambeau, C., Gasthaus, J., James, L., and Teh, Y. W. (2009).A stochastic memoizer for sequence data. In Proceedings of the 26th AnnualInternational Conference on Machine Learning, ICML ’09, page 1129–1136,New York, NY, USA. Association for Computing Machinery.Zarepour, M. and Al Labadi, L. (2012). On a rapid simulation of the Dirichletprocess. Statistics & Probability Letters, 82(5):916–924.Zemel, R. and Hinton, G. E. (1993). Developing population codes by min-imizing description length. Advances in neural information processingsystems, 6.Zhang, J. and Dassios, A. (2023). Truncated Poisson–Dirichlet approxima-tion for Dirichlet process hierarchical models. Statistics and Computing,33(1):30.Zhou, M., Hannah, L., Dunson, D., and Carin, L. (2012). Beta-negativebinomial process and Poisson factor analysis. In Artificial Intelligence andStatistics, pages 1462–1471. PMLR.Zhu, P., Bouchard-Coˆte´, A., and Campbell, T. (2020). Slice sampling forgeneral completely random measures. In Conference on Uncertainty inArtificial Intelligence, pages 699–708. PMLR.Zipf, G. K. (2013). The psycho-biology of language: An introduction todynamic philology. Routledge.128
Thesis/Dissertation
eng
Statistics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Campbell, Trevor (Professor of statistics)
Approximate posterior inference for Bayesian nonparametrics with total variation error bounds
Text