Chen Chen

Adrien Ugon

Xun Zhang

Amara Amara

Patrick Garda

Jean-Gabriel Ganascia

Carole Philippe

Andrea Pinna

Personalized Sleep Staging System using Evolutionary Algorithm and Symbolic Fusion

This paper presents a novel system for automatic sleep staging based on evolutionary technique and symbolic intelligence. Proposed system mimics decision making process of clinical sleep staging using Symbolic Fusion and considers personal singularity with an adaptive thresholds setting up system using Evolutionary Algorithm. It proved to be an effective and promising system in personalizing sleep staging. This system can also be integrated with other medical systems to realize remote sleep monitoring or home-care.

I. INTRODUCTION

Sleep is an essential part that contributes to refreshment and self-repairs in our daily life. However, the prevalence of sleep disorders is approximate to 10% of the population [START_REF] Ram | Prevalence and impact of sleep disorders and sleep habits in the united states[END_REF], including sleep apnea, insomnia and so on. Sleep disorders deteriorate quality of life [START_REF] D'ambrosio | Quality of life in patients with obstructive sleep apnea : Effect of nasal continuous positive airway pressure-a prospective study[END_REF] and become a significant cause of morbidity and mortality [START_REF] Punjabi | The epidemiology of adult obstructive sleep apnea[END_REF], [START_REF] Arzt | Association of sleep-disordered breathing and the occurrence of stroke[END_REF]. Therefore, it is crucial to diagnose and treat sleep disorders in time. That's why a clinical sleep analysis is needed and sleep staging is its fundamental step.

Polysomnography (PSG) is the gold standard diagnostic test for sleep disorders. It consists in recording different physiological signals during a sleep period. Resulting curves are visual analysed by sleep experts who assess sleep cycles using a subset of the PSG signals (brain activity (3 derivations EEG), eyes movement (2 EOG) and chin muscle tone (1 EMG)). Guidelines for the visual scoring of sleep and associated events have been defined by the American Academy of Sleep Medicine (AASM) [START_REF] Of | AASM manual for the scoring of sleep and associated events[END_REF]. Five sleep stages are defined : W (Wakefulness), N1, N2, N3 and REM (Rapid Eye Movements). The whole recording is divided in 30seconds epochs. Each epoch should be assigned to a sleep stage. Each sleep stage is defined by rules based on a list of criteria, combining the observation of specific visual patterns on the PSG signals. AASM guidelines provide definitions of sleep events and rules to score sleep stages. Visual patterns to identify sleep stages are defined by frequency ranges, amplitude, or other parameters that can be measured on the signals using standard signal processing methods. Although guidelines have been published to explain how to score sleep and Marie Curie University, Paris, France. {first name.last name}@lip6.fr 2 Xun Zhang and Amara Amara are with Institut Supérieur d'Electronique de Paris, MINARC, Paris, France. {first name.last name}@isep.fr 3

Carole

Philippe is with Hôpital Universitaire Pitié-Salpêtrière, Unité des Pathologie du Sommeil, Paris, France.

carole.philippe@aphp.fr in general situation, it has been shown that a variability between subjects in the shape of recording waves exist [START_REF] Jo | Genetic fuzzy classifier for sleep stage identification[END_REF].

Visual scoring of sleep stages is a tedious and timeconsuming task. The inter and intra scorers variability is a real issue [START_REF] Kuna | Agreement in computer-assisted manual scoring of polysomnograms across sleep centers[END_REF], although partially solved by the publication of international guidelines. However, it is still acknowledged that scoring practices need to be improved [START_REF] Penzel | Inter-scorer reliability between sleep centers can teach us what to improve in the scoring rules[END_REF]. To support it, automatic sleep staging methods have gained a wide spread attention in researches. However, the automated systems that on the market are not accurate enough to make the sleep expert sufficiently confident in their use. How to win the trust of the sleep expert ? The symbolic fusion is a decision method to abstract scattered raw data to build a rich knowledge. Symbolic fusion can be conducted through the 3-levels framework designed by B. Dasarathy [START_REF] Dasarathy | Sensor fusion potential exploitation-innovative architectures and illustrative applications[END_REF] : data level, features level and decision level. As a decision support method, symbolic fusion is faithful to the decision-making process of AASM manual. For each sleep stage, visual patterns should be extracted from the data, that correspond to one or more of the criteria of the AASM. Visual patterns are then combined, according to AASM rules, to take a decision about the sleep stage to assign to an epoch. An automatic sleep staging approach based on symbolic fusion has been proposed in [START_REF] Ugon | How to use symbolic fusion to support the sleep apnea syndrome diagnosis[END_REF]. The solution proposed by the authors still have some limitations. Some AASM criteria need to be added and fusion strategies need to be adapted to the new patterns within the symbolic fusion structure. This is due to the fact that the default symbolic fusion is still not enough « smart » to take into account the uniqueness of the person. In order to do that, an additional module should be developed.

Differential Evolution (DE) is a popular and efficient method of evolutionary algorithm which has been successfully applied to solve multi-parameter problems in diverse domains like mechanical engineering design or chemical engineering [START_REF] Storn | Differential evolution-a simple and efficient heuristic for global optimization over continuous spaces[END_REF]. It owns many advantages : 1) DE can mimic natural biological evolution and provide a fast and stable convergence ; 2) It is less sensitive to initial population ; 3) It is a parallel search method ; and 4) It can improve fitness function value iteratively.

In this paper we proposed a method to personalize the sleep staging using the symbolic fusion and differential evolution. We modified the Dasarathy symbolic structure by introducing an additional personalization module based on evolutionary algorithm between the Data and Feature level. This paper is organized as follows. Next section describes the limitations of symbolic fusion and presents a method to personalize sleep staging. Followed by results and comparison with the symbolic fusion proposed in [START_REF] Ugon | How to use symbolic fusion to support the sleep apnea syndrome diagnosis[END_REF]. Lastly, discussion and conclusion are provided.

II. METHODS

Before describing the evolutionary algorithm and explaining its implementation in the symbolic fusion, a quick overview on the limitations of symbolic fusion to personalize the sleep staging is presented.

A. Sleep Staging by Symbolic Fusion

An existing system based on symbolic fusion has been designed to realize classification of sleep stages [START_REF] Ugon | How to use symbolic fusion to support the sleep apnea syndrome diagnosis[END_REF]. This method uses the three-level architecture defined by B. Dasarathy to abstract data. Figure 1 presents the architecture of symbolic fusion : data fusion, feature fusion and decision fusion. In data fusion, time-domain analysis, frequency-domain analysis, and non-linear analysis have been applied to maximize the useful information and to minimize noise and artifacts. Nine digital parameters were extracted as shown in data fusion section in Figure 1. In feature fusion, digital parameters were transformed into feature parameters. In decision fusion, inference method is used to fulfill sleep staging on the basis of feature parameters. The link between digital and feature parameters is explained via thresholds in Table I.

However, visually interpretation is needed to perform thresholds setting-up before the feature fusion. These thresholds are used to differentiate the boundary of different feature parameters. The boundary for low, middle or high feature is determined by a visual setting of threshold for digital parameter. Total 15 thresholds are defined for translating nine digital parameters into 24 feature parameters. In this paper, an Adaptive Thresholds Setting Up (ATSU) system is proposed to automatically provide optimal values of thresholds.

B. Personalized Sleep Staging Workflow

In order to personalize the automatic sleep staging, we propose an ATSU system to set up the thresholds automatically using Evolutionary Algorithm. This algorithm is THRESHOLD SETTING UP adopted because of its robustness and ability in searching and discovering combinatorial thresholds [START_REF] Fogel | The advantages of evolutionary computation[END_REF]. With the feasibility of evolutionary technique and the elasticity of symbolic intelligence, steps of the new workflow of sleep stages detection is shown in Figure 2 ; details of ATSU system architecture are shown in Figure 3. The ATSU has the digital parameters extracted from the PSG signals as input, a first sleep detection by feature and decision fusion is made only on a subset of the epochs with a default value of the thresholds regardless of the person. The same epochs subset is also scored by the sleep expert. Then we compare the result of the first partial symbolic fusion with the sleep expert scoring and an assessment is made in the « Assessment of Automatic Sleep Staging » block. If the assessment is positive, the thresholds can be used for full sleep staging analysis by the feature and decision fusion. If the assessment is negative, we need to find new values of the thresholds in order to fit as much as possible to the partial sleep expert scoring. This is performed in the « Threshold Setting Up » block by an Evolutionary Algorithm. Once the new values of thresholds are generated, a new threshold adaptation cycle will start. Details are described in the following paragraph.

1) Differential Evolution Algorithm and Assessment of Automatic Sleep Staging Specification: Before using an evolutionary algorithm to optimize a set of values, it is necessary to define the population that will be used. Our population is composed of NP individuals, each described by a D-dimensional vector. In this work we used the classical strategy that is commonly noted DE/rand/1/bin. « rand » specifies that the vector to be mutated is selected randomly, « 1 » is the number of pairs of individuals used in construction of mutation vector, « bin » is the crossover scheme. Details of DE for ATSU are described as follows.

Initial population generation The initial population x i,G , i = 1, 2, ..., N P is generated randomly and is composed of 15-dimensional individuals which corresponds to 15 thresholds. Three different population sizes (NP=60, 90, 120) are compared with the analysis provided in Section III-B. Generation G is set to 50, which means if termination condition does not reach within 50 loops, it will terminate and provide the best threshold combination of the 50 th generation.

Fitness function evaluation F-Measure is used as the fitness function in this study and is evaluated in the Assessment Automatic Sleep Staging block. F-Measure balances both precision and recall as shown in Equation 1. Pr and Re represents Precision and Recall, respectively. Precision equals to (Tp/(Tp + Fp)) and Recall equals to (Tp/(Tp + Fn)). The parameters TP, FN, and FP are respectively True Positives, False Negatives, and False Positives, which are used to quantify the quality of classification results.

F -M easure = 2 * P recision * Recall P recision + Recall (1)
Check termination conditions The process terminates if the fitness function reaches a desired value(F-Measure=0.92, e.g as shown in Figure 4). Otherwise, the process continues on the mutation, crossover and selection.

Mutation A mutated population is generated, that will be used for crossover. Mutation uses the Equation 2. v i,G+1 is the i th individual of the mutated population of G+1 generation. The formula explains how to generate mutated individual v i,G+1 from 3 randomly selected individual in previous generation. Mutation Factor F is used to control the amplification of the differential variation (x r2,G -x r3,G ). r1, r2, r3 are randomly selected in [START_REF] Ram | Prevalence and impact of sleep disorders and sleep habits in the united states[END_REF]NP].

v i,G+1 = x r1,G + F * (x r2,G -x r3,G ) (2)
Crossover Crossover is used to decide if the j th parameter of i th individual of the crossed-over population should be mutated or not. Two kinds of random values are generated : random value j r (i) is a position generated for i th individual indicating that the j r (i)

th parameter of the i th individual will be mutated ; random value rb (i, j) is used to decide if the j th parameter of i th individual is mutated or not as shown in Equation 3.

u j,i,G+1 = v j,i,G+1 , if rb (i, j) ≤ CR or j = j r (i) x j,i,G , otherwise (3) where :
u j,i,G+1 is the j th parameter of i th individual in crossed-over population of G+1 generation ; -v j,i,G+1 is the j th parameter of i th individual in mutated population of G+1 generation ; -x j,i,G is the j th parameter of i th individual in initial (non-mutated) population of G generation ; -rb (i, j) is a uniform random number (rb

(i, j) ∈ [0, 1]) ; -CR is the crossover constant (CR ∈ [0, 1]) ; -j r (i) is the randomly generated parameter of i th individual, (j r (i) ∈ [1, 15]) ; -i is the individual (i ∈ [1, N P ]) ; -j is the parameter (i ∈ [1, 15]) ; -G is the generation.
Selection Greedy criterion is used to decide whether an individual will be a member of next generation. The crossedover individual is selected if the value of its fitness function is higher in comparison to fitness value of the initial individual, otherwise the initial individual is selected. In our case, the fitness function is the F-Measure f. This is formalized in Equation 4.

x i,G+1 = u i,G+1 , if f (u i,G+1 ) > f (x i,G ) x i,G+1 , otherwise (4)
III. RESULTS

A. Subjects and Recordings

Overnight PSG signals were recorded from 12 subjects (2 males and 10 females) ranging from 26 to 67 years old (mean = 53.42 years, std = 14 years) in Tenon hospital. PSG recordings of EEG, EOG and EMG were segmented into 30s epoch and manually scored by experts into five different stages : W, N1, N2, and N3 and REM according to AASM manual.

B. ATSU Evaluation

For each subject, 20% of the total epochs were randomly selected as a set and uniformly distributed in different stages to be the take into account during the ATSU. With this epochs selection, a population size of 90 for DE is used for adaptive thresholds selection. Figure 4 shows maximum F-Measure of Stage N3 in each generation with population sizes of 60, 90 and 120. Faster convergence is achieved with increase in the population size at the cost of greater computational complexity. Population size of 90 provides optimum between accuracy and computational complexity.

C. Comparison between Personalized and Non-Personalized Sleep Staging

We compared the F-Measure between the two sleep staging systems using symbolic fusion, one with a manual setting-up of thresholds and another with an automatic method based on the DE which can realize personalized sleep staging. In Table II, F-Measure for all the stages is increased for ATSU based personalized sleep staging in comparison to non-personalized sleep staging. For stage REM, F-Measure increases over 6% ; for stage W, N2 and N3, F-Measure increases 11%, 10% and 15 % respectively. These results are encouraging and prove the feasibility of a Personalized Sleep Staging System to help sleep experts in the PSG analysis. In this paper, we don't evaluate the accuracy of the symbolic fusion itself ( as already mentioned the symbolic fusion for sleep staging presented in [START_REF] Ugon | How to use symbolic fusion to support the sleep apnea syndrome diagnosis[END_REF] is still incomplete). With these results we prove how to make automatic sleep staging using the symbolic fusion according to AASM and taking into account the personal singularity based on differential evolution.

IV. DISCUSSION AND CONCLUSION

In the literature there are few works that try to validate an automated sleep staging system according to the AASM and the experience of sleep experts. Symbolic fusion and evolutionary algorithm proved to be a promising approach to help the experts in the PSG analysis. A new way to do PSG is proposed, Figure 5 shows the new workflow for personalized sleep staging and its possible use to help experts. Only a few number of epochs need to be analysed by expert instead of overall PSG examen ; and the remaining epochs are ananlysed by the automated sleep staging system. Moreover, an evolutionary algorithm is implemented in order to take into account the personal singularity. The combination of evolutionary technique and symbolic intelligence proposed in Personalized Sleep Staging System supplements the shortcomings of symbolic fusion and overcomes personal singularity limitation. This is a first proof of concept, for the further work, we need to define an effective way to choose the good number of generations in order to optimize the computational time for the thresholds setting up (e.g. for G=50 we have 37.81 sec with a Intel Core 5@1.6 GHz and 4 GB of RAM). We also need to understand, with the sleep experts, more efficient way to choose the optimal epochs for the partial clinical sleep staging.

Fig. 1 .

 1 Fig. 1. Sleep staging workflow using symbolic fusion in accordance with AASM

Fig. 2 .Fig. 3 .

 23 Fig. 2. Personalized sleep staging workflow

Fig. 4 .

 4 Fig. 4. F-Measure in each generation for ATSU

Fig. 5 .

 5 Fig. 5. The new workflow for personalized sleep staging

TABLE II

 II

	F-MEASURE COMPARISON

Chen Chen, Adrien Ugon, Patrick Garda, Jean-Gabriel Ganascia and Andrea Pinna are with Laboratoire d'Informatique de Paris 6, Pierre

ACKNOWLEDGMENT This work is financially supported by China Scholarship Council (CSC) and under collaboration of Laboratoire d'Informatique de Paris 6 and Institut Supérieur d'Electronique de Paris.


Julyan Arbel
email: julyan.arbel@carloalberto.org

Igor Prünster

A moment-matching Ferguson & Klass algorithm *

Keywords: Bayesian Nonparametrics, Completely random measures, Ferguson & Klass algorithm, Moment-matching, Normalized random measures, Posterior sampling

. This can immediately be turned into an algorithm for sampling realizations of CRMs or more elaborate models involving transformed CRMs. However, concrete implementation requires to truncate the random series at some threshold resulting in an approximation error. The goal of this paper is to quantify the quality of the approximation by a moment-matching criterion, which consists in evaluating a measure of discrepancy between actual moments and moments based on the simulation output. Seen as a function of the truncation level, the methodology can be used to determine the truncation level needed to reach a certain level of precision. The resulting moment-matching Ferguson & Klass algorithm is then implemented and illustrated on several popular Bayesian nonparametric models.

INTRODUCTION

Independent increment processes or, more generally, completely random measures (CRMs) are ubiquitous in modern stochastic modeling and inference. They form the basic building block of countless popular models in, e.g., Finance, Biology, Reliability, Survival Analysis. Within Bayesian nonparametric statistics they play a pivotal role. The Dirichlet process, the cornerstone of the discipline introduced in [START_REF] Ferguson | A Bayesian analysis of some nonparametric problems[END_REF], can be obtained as normalization or exponentiation of suitable CRMs (see [START_REF] Ferguson | Prior distributions on spaces of probability measures[END_REF]. Moreover, as shown in [START_REF] Lijoi | Models beyond the Dirichlet process[END_REF], CRMs can be seen as the unifying concept of a wide variety of Bayesian nonparametric models. See also [START_REF] Jordan | Hierarchical models, nested models and completely random measures[END_REF]. The concrete implementation of models based on CRMs often requires to simulate their realizations. Given they are discrete infinite objects, i≥1 J i δ Z i , some kind of truncation is required, producing an approximation error i≥M +1 J i δ Z i . Among the various representations useful for simulating realizations of CRMs the method due to [START_REF] Ferguson | A representation of independent increment processes without Gaussian components[END_REF] and popularized by [START_REF] Walker | Miscellanea. Representations of Lévy processes without Gaussian components[END_REF] stands out in that, for each realization, the weights J i 's are sampled in decreasing order. This clearly implies that for a given truncation level M the approximation error over the whole sample space is minimized. The appealing feature of decreasing jumps has lead to a huge literature exploiting the Ferguson & Klass algorithm. Limiting ourselves to recall contributions within Bayesian Nonparametrics we mention, among others, [START_REF] Argiento | A blocked Gibbs sampler for NGG-mixture models via a priori truncation[END_REF][START_REF] Argiento | A priori truncation method for posterior sampling from homogeneous normalized completely random measure mixture models[END_REF]; [START_REF] Barrios | Modeling with normalized random measure mixture models[END_REF][START_REF] De Blasi | A class of neutral to the right priors induced by superposition of beta processes[END_REF]; [START_REF] Epifani | Exponential functionals and means of neutral-tothe-right priors[END_REF]; [START_REF] Griffin | Posterior simulation of normalized random measure mixtures[END_REF]; [START_REF] Griffin | An adaptive truncation method for inference in bayesian nonparametric models[END_REF]; Nieto-Barajas and [START_REF] Nieto-Barajas | Markov Beta and Gamma Processes for Modelling Hazard Rates[END_REF]; Nieto- Barajas et al. (2004); Nieto-Barajas and Walker (2004); Nieto- [START_REF] Nieto-Barajas | A sensitivity analysis for Bayesian nonparametric density estimators[END_REF]; [START_REF] Nieto-Barajas | Bayesian semiparametric analysis of short-and long-term hazard ratios with covariates[END_REF]. General references dealing with the simulation of Lévy processes include [START_REF] Rosiński | Series representations of Lévy processes from the perspective of point processes[END_REF] and [START_REF] Cont | Financial modelling with jump processes[END_REF], who review the Ferguson & Klass algorithm and the compound Poisson process approximation to a Lévy process.

However, the assessment of the quality of the approximation due to the truncation for general CRMs is limited to some heuristic criteria. For instance, [START_REF] Barrios | Modeling with normalized random measure mixture models[END_REF] implement the Ferguson & Klass algorithm for mixture models by using the so called relative error index. The corresponding stopping rule prescribes to truncate when the relative size of an additional jump is below a pre-specified fraction of the sum of sampled jumps. The inherent drawbacks of such a procedure and related heuristic threshold-type procedures employed in the several of the above references is two-fold. On the one hand the threshold is clearly arbitrary without quantifying the total mass of the ignored jumps. On the other hand the total mass of the jumps beyond the threshold, i.e. the approximation error, can be very different for different CRMs or, even, for the same CRM with different parameter values; this implies that the same threshold can produce very different approximation errors in different situations. Starting from similar concerns about the quality of the approximation, the recent paper by [START_REF] Griffin | An adaptive truncation method for inference in bayesian nonparametric models[END_REF] adopts an algorithmic approach and proposes an adaptive truncation sampler based on sequential Monte Carlo for infinite mixture models based on normalized random measures and on stick-breaking priors. The measure of discrepancy that is used in order to assess the convergence of the sampler is based on the effective sample size (ESS) calculated over the set of particles: the algorithm is run until the absolute value of the difference between two consecutive ESS gets under a pre-specified threshold. Also motivated by the same concerns, [START_REF] Argiento | A blocked Gibbs sampler for NGG-mixture models via a priori truncation[END_REF][START_REF] Argiento | A priori truncation method for posterior sampling from homogeneous normalized completely random measure mixture models[END_REF] adopt an interesting approach to circumvent the problem of truncation by changing the model in the sense of replacing the CRM part of their model with a Poisson process approximation, which having an (almost surely) finite number of jumps can be sampled exactly. However, this leaves the question of the determination of the quality of approximation for truncated CRMs open. Another line of research, originated by [START_REF] Ishwaran | Gibbs sampling methods for stick-breaking priors[END_REF], is dedicated to validating the trajectories from the point of view of the marginal density of the observations in mixture models. In this context, the quality of the approximation is measured by the L 1 distance between the marginal densities under truncated and non-truncated priors. Recent interesting contributions in this direction include bounds for a Ferguson & Klass representation of the beta process [START_REF] Doshi | Variational inference for the Indian buffet process[END_REF] and bounds for the beta process, the Dirichlet process as well as for arbitrary CRMs in a size biased representation [START_REF] Paisley | Stick-breaking beta processes and the Poisson process[END_REF][START_REF] Campbell | Truncated completely random measures[END_REF].

This paper faces the problem by a simple yet effective idea. In contrast to the above strategies, our approach takes all jumps of the CRMs into account and hence leads to select truncation levels in a principled way, which vary according to the type of CRM and its parameters. The idea is as follows: given moments of CRMs are simple to compute, one can quantify the quality of the approximation by evaluating some measure of discrepancy between the actual moments of the CRM at issue (which involve all its jumps) and the "empirical" moments, i.e. the moments computed based on the truncated sampled realizations of the CRM. By imposing such a measure of discrepancy not to exceed a given threshold and selecting the truncation level M large enough to achieve the desired bound, one then obtains a set of "validated" realizations of the CRM, or, in other terms, satisfying a moment-matching criterion. An important point to stress is that our validation criterion is all-purpose in spirit since it aims at validating the CRM samples themselves rather than samples of a transformation of the CRM. Clearly the latter type of validation would be ad hoc, since it would depend on the specific model. For instance, with the very same set of moment-matching realizations of a gamma process, one could obtain a set of realizations of the Dirichlet process via normalization and a set gamma mixture hazards by combination with a suitable kernel. Moreover, given moments of transformed CRMs are typically challenging to derive, a moment-matching strategy would not be possible in most cases. Hence, while the quantification of the approximation error does not automatically translate to transformed CRMs, one can still be confident that the momentmatching output at the CRM level produces good approximations. That this is indeed the case is explicitly shown in some practical examples both for prior and posterior quantities in Section 3.

The outline of the paper is as follows. In Sections 2.1-2.2 we recall the main properties of CRMs and provide expressions for their moments. In Sections 2.3-2.4 we describe the Ferguson & Klass algorithm and introduce the measure of discrepancy between moments used to quantify the approximation error due to truncation. Section 3 illustrates the moment-matching Ferguson & Klass algorithm for some popular CRMs and CRM-based Bayesian nonparametric models, namely normalized CRMs and the beta-stable Indian buffet process. Some probabilistic results, discussed in Section 2.3, are given in the Appendix.

COMPLETELY RANDOM MEASURES

Definition and main properties

Let M X be the set of boundedly finite measures on X, which means that if µ ∈ M X then µ(A) < ∞ for any bounded set A. X is assumed to be a complete and separable metric space and both X and M X are equipped with the corresponding Borel σ-algebras. See [START_REF] Daley | An introduction to the theory of point processes[END_REF] for details.

Definition 1. A random element μ, defined on (Ω, F , P) and taking values in M X , is called a completely random measure (CRM) if, for any collection of pairwise disjoint sets A 1 , . . . , A n in X, the random variables μ(A 1 ), . . . , μ(A n ) are mutually independent.

An important feature is that a CRM μ selects (almost surely) discrete measures and hence can be represented as

(1) μ = i≥1 J i δ Z i
where the jumps J i 's and locations Z i 's are random and independent. In (1) and throughout we assume there are no fixed points of discontinuity a priori. The main technical tool for dealing with CRMs is given by their Laplace transform, which admits a simple structural form known as Lévy-Khintchine representation.

In fact, the Laplace transform of μ(A), for any A in X, is given by

(2) L A (u) = E e -λμ(A) = exp - R + ×A
1 -e -λv ν(dv, dx)

for any λ > 0. The measure ν is known as Lévy intensity and uniquely characterizes μ. In particular, there corresponds a unique CRM μ to any measure ν on R + × X satisfying the integrability condition

(3) B R + min{v, 1}ν(dv, dx) < ∞
for any bounded B in X. From an operational point of view this is extremely useful, since a single measure ν encodes all the information about the jumps J i 's and the locations Z i 's. The measure ν will be conveniently rewritten as

(4) ν(dv, dx) = ρ(dv|x)α(dx),
where ρ is a transition kernel on R + × X controlling the jump intensity and α is a measure on X determining the locations of the jumps. If ρ does not depend on x, the CRM is said homogeneous, otherwise it is non-homogeneous. We now introduce two popular examples of CRMs that we will serve as illustrations throughout the paper.

Example 1. The generalized gamma process introduced by [START_REF] Brix | Generalized gamma measures and shot-noise Cox processes[END_REF] is characterized by a Lévy intensity of the form

(5) ν(dv, dx) = e -θv Γ(1 -γ)v 1+γ dv α(dx),
whose parameters θ ≥ 0 and γ ∈ [0, 1) are such that at least one of them is strictly positive. Notable special cases are: (i) the gamma CRM which is obtained by setting γ = 0; (ii) the inverse-Gaussian CRM, which arises by fixing γ = 0.5; (iii) the stable CRM which corresponds to θ = 0. Moreover, such a CRM stands out for its analytical tractability. In the following we work with θ = 1, a choice which excludes the stable CRM. This is justified in our setting because the moments of the stable process do not exist. See Remark 1.

Example 2. The stable-beta process, or three-parameter beta process, was defined by [START_REF] Teh | Indian buffet processes with power-law behavior[END_REF] as an extension of the beta process [START_REF] Hjort | Nonparametric bayes estimators based on beta processes in models for life history data[END_REF]. Its jump sizes are upper-bounded by 1 and its Lévy intensity on [0, 1] × X is given by

(6) ν(dv, dx) = Γ(c + 1) Γ(1 -σ)Γ(c + σ) v -σ-1 (1 -v) c+σ-1 dv α(dx),
where σ ∈ [0, 1) is termed discount parameter and c > -σ concentration parameter. When σ = 0, the stable-beta process reduces to the beta CRM of [START_REF] Hjort | Nonparametric bayes estimators based on beta processes in models for life history data[END_REF]. Moreover, if c = 1 -σ, it boils down to a stable CRM where the jumps larger than 1 are discarded.

Moments of a CRM

For any measurable set A of X, the n-th (raw) moment of μ(A) is defined by

m n (A) = E μn (A) .
In the sequel the multinomial coefficient is denoted by

n k 1 •••kn = n! k 1 !...kn! .
In the next proposition we collect known results about moments of CRMs which are crucial for our methodology.

Proposition 1. Let μ be a CRM with Lévy intensity ν(dv, dx). Then the i-th cumulant of μ(A), denoted by κ i (A), is given by

κ i (A) = R + ×A v i ν(dv, dx),
which, in the homogeneous case ν(dv, dx) = ρ(dv)α(dx), simplifies to

κ i (A) = α(A) ∞ 0 v i ρ(dv).
The n-th moment of μ(A) is given by

m n (A) = ( * ) ( n k 1 •••kn ) n i=1 κ i (A)/i! k i ,
where the sum ( * ) is over all n-tuples of nonnegative integers (k 1 , . . . , k n ) satisfying the constraint

k 1 + 2k 2 + • • • + nk n = n.
A proof is given in the Appendix.

In the following we focus on (almost surely) finite CRMs i.e. μ(X) < ∞. This is motivated by the fact that most Bayesian nonparametric models, but also models in other application areas, involve finite CRMs. Hence, we assume that the measure α in (3) is finite i.e. α(X) := a ∈ (0, ∞). This is a sufficient condition for μ(X) < ∞ in the non-homogeneous case and also necessary in the homogeneous case (see e.g. [START_REF] Regazzini | Distributional results for means of normalized random measures with independent increments[END_REF]. A common useful parametrization of α is then given as aP * with P * a probability measure and a a finite constant. Note that, if μ(X) = ∞, one could still identify a bounded set of interest A and the whole following analysis carries over by replacing μ(X) with μ(A).

As we shall see in Section 2.3, the key quantity for evaluating the truncation error is given by the random total mass of the CRM, μ(X). Proposition 1 shows how the moments m n = m n (X) can be obtained from the cumulants κ i = κ i (X) and, in particular, the relations between the first four moments and the cumulants are

m 1 = κ 1 , m 2 = κ 2 1 + κ 2 , m 3 = κ 3 1 + 3κ 1 κ 2 + κ 3 , m 4 = κ 4 1 + 6κ 2 1 κ 2 + 4κ 1 κ 3 + 3κ 2 2 + κ 4 . CRM Cumulants Moments κi m1 m2 m3 m4 G a(i -1)! a a (2) a (3) a (4) IG a(1/2) (i-1) a a 2 + 1 2 a a 3 + 3 2 a 2 a 4 + 3a 3 + 3 4 a + 15 4 a 2 + 15 8 a GG a(1 -γ) (i-1) a a 2 + a(1-γ) a 3 + 3a 2 (1-γ) a 4 + 6a 3 (1-γ) +a(1-γ) (2) +a 2 (1-γ)(11-7γ) + a(1-γ) (3) B a (i-1)! (c+1) (i-1) a a 2 + a c+1 a 3 + 3a 2 c+1 a 4 + 6a 3 c+1 + 8a 2 (c+1) (2) + 2a (c+1) (2) + 3a 2 (c+1) 2 + 6a (c+1) (3) SB a (1-σ) (i-1) (c+1) (i-1) a a 2 + a 1-σ c+1 a 3 + 3a 2 1-σ c+1 a 4 + 6a 3 1-σ c+1 + 4a 2 (1-σ) (2) (c+1) (2) +a (1-σ) (2) (c+1) (2) +3a 2 (1-σ) 2 (c+1) 2 + a (1-σ) (3) (c+1) (3)
Table 1 Cumulants and first four moments of the random total mass μ(X) for the gamma (G), inverse-Gaussian (IG), generalized gamma (GG), beta (B) and stable-beta (SB) CRMs.

With reference to the two examples considered in Section 2.1, in both cases the expected value of μ(X) is a, which explains the typical terminology total mass parameter attributed to a. For the generalized gamma CRM the variance is given by Var(μ(X)) = a(1-γ), which shows how the parameter γ affects the variability. Moreover,

κ i = a(1 -γ) (i-1) with x (k) = x(x + 1) . . . (x + k -1)
denoting the ascending factorial. As for the stable-beta CRM, we have Var(μ(X)) = a 1-σ c+1 with both discount and concentration parameter affecting the variability, and also

κ i = a (1-σ) (i-1) (1+c) (i-1)
. Table 1 summarizes the cumulants κ i and moments m n for the random total mass μ(X) for the generalized gamma (assuming as in Example 1 θ = 1), stable-beta CRMs and some of their special cases.

Remark 1. The stable CRM, which can be derived from the generalized gamma CRM by setting θ = 0, does not admit moments. Hence, it cannot be included in our moment-matching methodology. However, the stable CRM with jumps larger than 1 discarded, derived from the stable-beta process by setting c = 1 -σ, has all moments. Moreover, even when working with the standard stable CRM, posterior quantities typically involve an exponential updating of the Lévy intensity (see [START_REF] Lijoi | Models beyond the Dirichlet process[END_REF], which makes the corresponding moments finite. This then allows to apply the moment matching methodology to the posterior.

Ferguson & Klass algorithm

For notational simplicity we present the Ferguson & Klass algorithm for the case X = R. However, note that it can be readily extended to more general Euclidean spaces (see e.g. [START_REF] Orbanz | Unit-rate poisson representations of completely random measures[END_REF]. Given a CRM

(7) μ = ∞ i=1 J i δ Z i ,
the Ferguson & Klass representation consists in expressing random jumps J i occurring at random locations Z i in terms of the underlying Lévy intensity.

In particular, the random locations Z i , conditional on the jump sizes J i , are obtained from the distribution function F Z i |J i given by

F Z i |J i (s) = ν(dJ i , (-∞, s]) ν(dJ i , R) .
In the case of a homogeneous CRM with Lévy intensity ν(dv, dx) = ρ(dv) aP * (dx), the jumps are independent of the locations and, therefore F Z i |J i = F Z i implying that the locations are i.i.d. samples from P * . As far as the random jumps are concerned, the representation produces them in decreasing order, that is,

J 1 ≥ J 2 ≥ • • • . Indeed, they are obtained as ξ i = N (J i ),
where ∼ Exp(1). Therefore, the J i 's are obtained by solving the equations ξ i = N (J i ). In general, this is achieved by numerical integration, e.g., relying on quadrature methods (see, e.g. [START_REF] Burden | Numerical Analysis[END_REF]. For specific choices of the CRM, it is possible to make the equations explicit or at least straightforward to evaluate. For instance, if μ is a generalized gamma process (see Example 1), the function N takes the form ( 8)

N (v) = ν([v, ∞), R) is a decreasing function,
N (v) = a Γ(1 -γ) ∞ v e -u u -(1+γ) du = a Γ(1 -γ) Γ(v; -γ),
with Γ( • ; • ) indicating an incomplete gamma function. If μ is the stable-beta process, one has ( 9)

N (v) = a Γ(c + 1) Γ(1 -σ)Γ(c + σ) 1 v u -σ-1 (1-u) c+σ-1 du = a Γ(c + 1) Γ(1 -σ)Γ(c + σ) B(1-v; c+σ, -σ),
where B( • ; • , • ) denotes the incomplete beta function.

Hence, the Ferguson & Klass algorithm can be summarized as follows.

Algorithm 1 Ferguson & Klass algorithm

1: Sample ξi ∼ PP for i = 1, . . . , M 2: Define Ji = N -1 (ξi) for i = 1, . . . , M 3: Sample Zi ∼ P * for i = 1, . . . , M 4: Approximate μ by M i=1 JiδZ i
Since it is impossible to sample an infinite number of jumps, approximate simulation of μ is in order. This becomes a question of determining the number M of jumps to sample in (7) leading to the truncation

(10) μ ≈ μM = M i=1 J i δ Z i ,
with approximation error in terms of the un-sampled jumps equal to ∞ i=M +1 J i . The Ferguson & Klass representation has the key advantage of generating the jumps in decreasing order implicitly minimizing such an approximation error. Then, the natural path to determining the truncation level M would be the evaluation of the Ferguson & Klass tail sum (11)

∞ i=M +1 N -1 (ξ i ).
Brix (1999, Theorem A.1) provided an upper bound for (11) in the generalized gamma case. In Proposition 4 of Appendix B we derive also an upper bound for the tail sum of the stable-beta process. However, both bounds are far from sharp and therefore of little practical use as highlighted in Appendix B. This motivates the idea of looking for a different route and our proposal consists in the moment-matching technique detailed in the next section.

Moment-matching criterion

Our methodology for assessing the quality of approximation of the Ferguson & Klass algorithm consists in comparing the actual distribution of the random total mass μ(X) with its empirical counterpart, where by empirical distribution we mean the distribution obtained by the sampled trajectories, i.e. by replacing random quantities by Monte Carlo averages of their sampled trajectories. In particular, based on the fact that the first K moments carry much information about a distribution, theoretical and empirical moments of μ(X) are compared.

The infinite vector of jumps is denoted by J = (J i ) ∞ i=1 and a vector of jumps sampled by the Ferguson & Klass algorithm by

J (l) = (J (l) 1 , . . . , J (l) M ).
Here, l = 1, . . . , N FK stands for the l-th iteration of the algorithm, i.e. for the l-th sampled realization. We then approximate the expectation E of a statistic of the jumps, say S(J ), by the following empirical counterpart, denoted by E FK ,

(12) E S(J )] ≈ E FK S(J )] := 1 N FK N FK l=1
S J (l) .

Note that there are two layers of approximation involved in (12): first, only a finite number of jumps M is used; second, the actual expected value is estimated through an empirical average which typically conveys on Monte Carlo error. The latter is not the focus of the paper, so we take a large enough number of trajectories, N FK = 10 4 , in order to insure a limited Monte Carlo error of the order of 0.01. We focus on the first approximation inherent to the Ferguson & Klass algorithm.

More specifically, as far as moments are concerned, m K = (m 1 , . . . , m K ) denotes the first K moments of the random total mass μ(X) = ∞ i=1 J i provided in Section 2.2 and mK = ( m1 , . . . , mK ) indicates the first K empirical moments given by

(13) mn = E FK   M i=1 J i n   .
As measure of discrepancy between theoretical and empirical moments, a natural choice is given by the mean squared error between the vectors of moments or, more precisely, between the n-th roots of theoretical and empirical moments

(14) = (m K , mK ) = 1 K K n=1 m 1/n n -m1/n n 2 1/2 .
When using the Ferguson & Klass representation for computing the empirical moments the index depends on the truncation level M and we highlight such a dependence by using the notation M . Of great importance is also a related quantity, namely the number of jumps necessary for achieving a given level of precision, which essentially consists in inverting M and is consequently denoted by M ( ).

The index of discrepancy ( 14) clearly also depends on K, the number of moments used to compute it and 1/K in ( 14) normalizes the indices in order to make them comparable as K varies. A natural question is then about the sensitivity of ( 14) w.r.t. K. It is desirable for M to capture fine variations between the theoretical and empirical distributions, which is assured for large K. In extensive simulation studies not reported here we noted that increasing K in the range {1, . . . , 10} makes the index increase and then plateau and this holds for all processes and parameter specifications used in the paper. Recalling also the whole body of work by Pearson on eponymous curves, which shows that the knowledge of four moments suffices to cover a large number of known distributions, we adhere to his rule of thumb and choose K = 4 in our analyses. On the one hand it is a good compromise between targeted precision of the approximation and speed of the algorithm. On the other hand it is straightforward to check the results as K varies in specific applications; for the ones considered in the following sections the differences are negligible.

In the literature several heuristic indices based on the empirical jump sizes around the level of truncation have been discussed (cf Remark 3 in [START_REF] Barrios | Modeling with normalized random measure mixture models[END_REF]. Here, in order to compare such procedures with our moment criterion, we consider the relative error index which is based on the jumps themselves. It is defined as the expected value of the relative error between two consecutive partial sums of jumps. Its empirical counterpart is denoted by e M and given by ( 15)

e M = E FK J M M i=1 J i .

APPLICATIONS TO BAYESIAN NONPARAMETRICS

In this section we concretely implement the proposed moment-matching Ferguson & Klass algorithm to several Bayesian nonparametric models. The performance in terms of both a priori and a posteriori approximation is evaluated. A comparison of the quality of approximation resulting from using (15) as benchmark index is provided.

A priori simulation study

We start by investigating the performance of the proposed moment-matching version of the Ferguson & Klass algorithm w.r.t. the CRMs defined in Examples 1 and 2, namely the generalized gamma and stable-beta processes. Figure 1 displays the behaviour of both the moment-matching distance M (left panel) and the relative jumps' size index e M (right panel) as the truncation level M increases.

The plots, from top to bottom, correspond to: the generalized gamma process with varying γ and a = 1 fixed; the inverse-Gaussian process with varying total mass a (which is a generalized gamma process process with γ = 0.5); the stablebeta process with varying discount parameter σ and a = 1 fixed.

First consider the behaviour of the indices as the parameter specifications vary. It is apparent that, for any fixed truncation level M , the indices M and e M increase as each of the parameters a, γ or σ increases. For instance, roughly speaking, a total mass parameter a corresponds to sampling trajectories defined on the interval [0, a] (see [START_REF] Regazzini | Distributional results for means of normalized random measures with independent increments[END_REF], and a larger interval worsens the quality of approximation for any given truncation level. Also it is natural that γ and σ impact in similar way M and e M given they stand for the "stable" part of the Lévy intensity. See first and third rows of Figure 1.

As far as the comparison between M and e M is concerned, it is important to note that e M consistently downplays the error of approximation related to the truncation. This can be seen by comparing the two columns of Figure 1. M is significantly more conservative than e M for both the generalized gamma and the stable-beta processes, especially for increasing values of the parameters γ, a or σ. This indicates quite a serious issue related to e M as a measure for the quality of approximation and one should be cautious when using it. In contrast, the moment-matching index M matches more accurately the known behaviour of these processes as the parameters vary.

By reversing the viewpoint and looking at the truncation level M ( ) needed for achieving a certain error of approximation in terms of moment-match, the results become even more intuitive. We set = 0.1 and computed M ( ) on a grid of size 20 × 20 with equally-spaced points for the parameters (a, γ) ∈ (0, 2) × (0, 0.8) for the generalized gamma process and (a, c) ∈ (0, 2) × (0, 30) for the beta process. Figure 2 displays the corresponding plots. In general, it is interesting to note that a limited number of jumps is sufficient to achieve good precision levels. Analogously to Figure 1, larger values of the parameters require a larger number of jumps to achieve a given precision level. In particular, when γ > 0.5, one needs to sample a significantly larger number of jumps. For instance, in the generalized gamma process case, with a = 1, the required number of jumps increases from 28 to 53 when passing from γ = 0.5 to γ = 0.75. It is worth noting that for the normalized version of the generalized gamma process, to be discussed in Section 3.2 and quite popular in applications, the estimated value of γ rarely exceeds 0.75 in species sampling, whereas it is typically in the range [0.2, 0.4] in mixture modeling.

Normalized random measures with independent increments

Having illustrated the behaviour of the moment-matching methodology for plain CRMs we now investigate it on specific classes of nonparametric priors, which typically involve a transformation of the CRM. Moreover, given their posterior distributions involve updated CRMs it is important to test the momentmatching Ferguson & Klass algorithm also on posterior quantities. The first class of models we consider are normalized random measures with independent increments (NRMI) introduced by [START_REF] Regazzini | Distributional results for means of normalized random measures with independent increments[END_REF]. Such nonparametric priors .6 q q q q q q q q q q q q q q q q q q q q q q q q q q GG : a .3 q q q q q q q q q q q q q q q q q q q q q q q q q GG : a .6 q q q q q q q q q q q q q q q q q q q q q q q q q q IG : γ = 0.5

a = 0.5 a = 1 a = 2 M eM 0 5 10 15 20 25 0 .1 .2
.3 q q q q q q q q q q q q q q q q q q q q q q q q q IG : γ = 0.5

a = 0.5 a = 1 a = 2 M l M 0 5 10 15 20 25 0 .2 .4
.6 q q q q q q q q q q q q q q q q q q q q q q q q q SBP : a = 1, c = 0.5 σ = 0 σ = 0.25 σ = 0.5 .3 q q q q q q q q q q q q q q q q q q q q q q q q q SBP : a = 1, c = 0.5 σ = 0 σ = 0.25 σ = 0.5 M . Left panel: generalized gamma process for a ∈ (0, 2) and γ ∈ (0, 0.8). Right panel: beta process for a ∈ (0, 2) and c ∈ (0, 30).

(b) eM

then an NRMI is defined as

(16) P = μ μ(X)
.

Particular cases of NRMI are then obtained by specifying the CRM in (16). For instance, by picking the generalized gamma process defined in Example 1 one obtains the normalized generalied gamma process, denoted by NGG, and first used in a Bayesian context by [START_REF] Lijoi | Controlling the reinforcement in Bayesian non-parametric mixture models[END_REF].

Posterior Distribution of an NRMI

The basis of any Bayesian inferential procedure is represented by the posterior distribution. In the case of NRMIs, the determination of the posterior distribution is a challenging task since one cannot rely directly on Bayes' theorem (the model is not dominated) and, with the exception of the Dirichlet process, NRMIs are not conjugate as shown in [START_REF] James | Conjugacy as a distinctive feature of the Dirichlet process[END_REF]. Nonetheless, a posterior characterization has been established in [START_REF] James | Posterior analysis for normalized random measures with independent increments[END_REF] and it turns out that, even though NRMIs are not conjugate, they still enjoy a sort of "conditional conjugacy." This means that, conditionally on a suitable latent random variable, the posterior distribution of an NRMI coincides with the distribution of an NRMI having fixed points of discontinuity located at the observations. Such a simple structure suggests that when working with a general NRMI, instead of the Dirichlet process, one faces only one additional layer of difficulty represented by the marginalization with respect to the conditioning latent variable.

Before stating the posterior characterization to be used with our algorithm, we need to introduce some notation and basic facts. Let (Y n ) n≥1 be an exchangeable sequence directed by an NRMI, i.e. 18)

Y
f U |Y (u) ∝ u n-1 exp -ψ(u) k j=1 τ n j u|Y * j ,
where τ n j (u|Y * j ) = ∞ 0 v n j e -uv ρ(dv|Y * j ) and ψ is the Laplace exponent of μ defined by ψ(u) = -log L X (u) , cf (2). Finally, assume P * = E[ P ] to be nonatomic.

Proposition 2 [START_REF] James | Posterior analysis for normalized random measures with independent increments[END_REF]

(21) f * j (v) ∝ v n j e -uv ρ dv|Y * j ,
(c) μ * and J * j , j = 1, . . . , k, are independent. Moreover, the posterior distribution of the NRMI P , conditional on U , is given by

(22) [ P |U, Y] d = w μ * μ * (X) + (1 -w) k k=1 J * j δ Y * j k l=1 J * l ,
where w = μ * (X)/(μ * (X) + k l=1 J * l ).

In order to simplify the notation, in the statement we have omitted explicit reference to the dependence on [U |Y] of both μ * and {J * j : j = 1, . . . , k}, which is apparent from ( 20) and ( 21). A nice feature of the posterior representation of Proposition 2 is that the only quantity needed for deriving explicit expressions for particular cases of NRMI is the Lévy intensity (4). For instance, in the case of the generalized gamma process, the CRM part μ * in ( 19) is still a generalized gamma process characterized by a Lévy intensity of the form of ( 5)

(23) ν * (dv, dy) = e -(1+u)v Γ(1 -γ)v 1+γ dv aP * (dy).
Moreover, the distribution of the jumps ( 21) corresponding to the fixed points of discontinuity Y * j 's in ( 19) reduces to a gamma distribution with density

(24) f * j (v) = (1 + u) n j -γ Γ(n j -γ) v n j -γ-1 e -(1+u)v .
Finally, the conditional distribution of the non-negative latent variable U given Y (18) is given by ( 25)

f U |Y (u) ∝ u n-1 (u + 1) kγ-n exp - a γ (u + 1) γ .
The availability of this posterior characterization makes it then possible to determine several important quantities such as the predictive distributions and the induced partition distribution. See [START_REF] James | Posterior analysis for normalized random measures with independent increments[END_REF] for general NRMI and Lijoi et al. (2007) for the subclass of normalized generalized gamma processes. See also [START_REF] Argiento | A blocked Gibbs sampler for NGG-mixture models via a priori truncation[END_REF] for another approach to approximate the NGG with a finite number of jumps. From (19) it is apparent that the posterior of the unnormalized CRM μ, conditional on the latent variable U , is composed of the independent sum of a CRM μ * and fixed points of discontinuity at the distinct observations Y * j . The part which is at stake here is obviously μ * for which only approximate sampling is possible. As for the fixed points of discontinuities, they are independent from μ * and can be sampled exactly, at least in special cases.

Moment-matching for posterior NRMI

We focus on the case of the NGG process. By (20) the Lévy intensity of μ * is obtained by exponentially tilting the Lévy intensity of the prior μ. Hence, the Ferguson & Klass algorithm applies in the same way as for the prior. The sampling of the fixed points jumps is straightforward from the gamma distributions (24). As far as the moments are concerned, key ingredient of our algorithm, the cumulants of μ * are equal to

κ * i = a (1-γ) (i-1)
(u+1) i-γ and the corresponding moments are then obtained via Proposition 1.

Our simulation study is based on a sample of size n = 10. Such a small sample size is challenging in the sense that the data provide rather few information and the CRM part of the model is still prevalent. We examine three possible clustering configurations of the observations Y * i s: (i) k = 1 group, with n 1 = 10, (ii) k = 3 groups, with n 1 = 1, n 2 = 3, n 3 = 6, and (iii) k = 10 groups, with n j = 1 for j = 1, . . . , 10. First let us consider the behaviour of f U |Y , which is illustrated in Figure 3 for n = 10 and k ∈ {1, 2, . . . , 10}. It is clear that the smaller the number of clusters, the more f U |Y is concentrated on small values, and vice versa. Now we consider μ * (X), the random total mass corresponding to the CRM part of the posterior only given in ( 22). Such a quantity depends on U whose distribution is driven by the data Y. In order to keep the presentation as neat as possible, and in the same time to remain consistent with the data, we choose to condition on U = u for u equal to the mean of f U |Y , the most natural representative value. Given this, it is possible to run the Ferguson & Klass algorithm on the CRM part μ * of the posterior and compute moment-matching index M as the number of jumps varies. Figure 4 shows these results for the inverse-Gaussian CRM, a special case of the generalized gamma process corresponding to γ = 0.5. Such posteriors were sampled under the above mentioned Y clustering configuration scenarios (i)-(iii), which led to mean values of U |Y of, respectively, 6.3, 8.9 and 25.1. The plot also displays a comparison to the prior values of M and indicates that for a given number of jumps the approximation error, measured in terms of M , is smaller for the posterior CRM part μ * w.r.t. to the prior CRM μ.

Additionally, instead of considering only the CRM part μ * of the posterior, one may be interested in the quality of the full posterior which includes also the fixed discontinuities. For this purpose we consider an index which is actually of interest in its own. In particular, we evaluate the relative importance of the CRM part w.r.t. the part corresponding to the fixed points of discontinuity in terms of the ratio E k j=1 J * j /E μ * (X) . Loosely speaking one can think of the numerator as the expected weight of the data and the denominator as the expected weight of the prior. Recall that in the NGG case, for a given pair (n, k) and conditional on U = u, the sum of fixed location jumps is a gamma(n -kγ, u + 1). Hence, the index becomes

(26) E k j=1 J * j |U = u E μ * (X)|U = u = (n -kγ)/(u + 1) a/(u + 1) 1-γ = n -kγ a(u + 1) γ .
By separately mixing the conditional expected values in (26) over f U |Y (we use an adaptive rejection algorithm to sample from f U |Y ) we obtained the results summarized in the table of Figure 4. We can appreciate that the fixed part typically overcomes (or is at least of the same order than) the CRM part, a phenomenon which uniformly accentuates as the sample size n increases. Returning to the original problem of measuring the quality of approximation in terms of moment matching, these findings make it apparent that the comparative results of Figure 4 between prior and posterior are conservative. In fact, if performing the moment-match on the whole posterior, i.e. including the fixed jumps which can be sampled exactly, the corresponding moment-matching index would, for any given truncation level M , indicate a better quality of approximation w.r.t. the index based solely on μ * . Note that computing the moments of μ * (X) + k i=1 J i straightforward given the independence between μ * and the fixed jumps J i 's and also among the jumps themselves. From a practical point of view the findings of this section suggest that a given quality of approximation in terms of momentmatch for the prior represents an upper bound for the quality of approximation in the posterior. .6 q q q q q q q q q q q q q q q q q q q q q q q q q q IG : γ = 0.5, a = 1

Prior k = 1 k = 3 k = 10 E k j=1 J * j /E μ

A note on the inconsistency for diffuse distributions

In the context of Gibbs-type priors, of which the normalized generalized gamma process is a special case, De [START_REF] De Blasi | An asymptotic analysis of a class of discrete nonparametric priors[END_REF] showed that, if the data are generated from a "true" P 0 , the posterior of P concentrates at a point mass which is the linear combination bP * (•) + (1 -b)P 0 (•) of the prior guess P * = E( P ) and P 0 . The weight b depends on the prior and, indirectly, on P 0 , since P 0 dictates the rate at which the distinct observations k are generated. For a diffuse P 0 , all observations are distinct and k = n (almost surely). In the NGG case this implies that b = γ and hence the posterior is inconsistent since it does not converge to P 0 . For the inverse-Gaussian process, i.e. with γ = 0.5, the posterior distribution gives asymptotically the same weight to P * and P 0 . The last row of the table of Figure 4, which displays the ratio E k j=1 J * j /E μ * (X) for k = n, is an illustration of this inconsistency result since the ratio gets close to 1 as n grows. In contrast, when P 0 is discrete, which implies that k increases at a slower rate than n, one always has consistency. This is illustrated by the first two rows of the table of Figure 4, where one can appreciate that the ratio E k j=1 J * j /E μ * (X) increases as n increases, giving more and more weight to the data. These findings suggest that consistency issues for general NRMI could be explored from new perspectives based on the study of the asymptotic behavior of f U |Y , which will be subject to future work.

Stable-beta Indian buffet process

The Indian buffet process (IBP), introduced in [START_REF] Ghahramani | Infinite latent feature models and the Indian buffet process[END_REF], is one of the most popular models for feature allocation and is closely connected to the beta process discussed in Example 2. In fact, when marginalizing out the Dirichlet process and considering the resulting partition distribution one obtains the well known Chinese restaurant process. Likewise, as shown in [START_REF] Thibaux | Hierarchical beta processes and the Indian buffet process[END_REF], when integrating out a beta process in a Bernoulli process (BeP) model one obtains the IBP. Recall that a Bernoulli process, with an atomic base measure μ, is a stochastic process whose realizations are collections of atoms of mass 1, with possible locations given by the atoms of the base measure μ. Such an atom is element of the collection with probability given by the jump size in μ. Later, [START_REF] Teh | Indian buffet processes with power-law behavior[END_REF] generalized the construction and defined the stable-beta Indian buffet process as . . . , n, (27) μ|c, σ, aP * ∼ SBP(c, σ, aP * ).

Y i |μ i.i.d. ∼ BeP(μ) for i = 1,
Given the construction involves a CRM, it is clear that any conditional simulation algorithm will need to rely on some truncation for which we use our moment-matching Ferguson & Klass algorithm.

Posterior distribution in the IBP

Let us consider a conditional iid sample Y = (Y 1 , . . . , Y n ) as in (27). Note that due to the discreteness of μ, ties appear with positive probability. We adopt the same notations for the ties Y * j and frequencies n j as in Section 3.2. Then we can state the following result which highlights the posterior structure of the stable-beta process in the Indian buffet process.

Proposition 3 [START_REF] Teh | Indian buffet processes with power-law behavior[END_REF]. Let (Y n ) n≥1 be as in (27). Then the posterior distribution of μ conditional on Y is given by the distribution of

μ * + k j=1 J * j δ Y * j
where (a) μ * is a stable-beta process characterized by the Lévy intensity

ν * (dv, dx) = (1 -v) n ν(dv, dx), (b) the jump height J * j corresponding to Y * j is beta distributed J * j ∼ beta(n j -σ, c + σ + n -n j ), (
c) μ * and J * j , j = 1, . . . , k, are independent.

Note that due to the polynomial tilting of ν by (1 -u) n in (a) above, the CRM part μ * is still a stable-beta process with updated parameters

c * = c + n and a * = a (c + σ) (n) (c + 1) (n) ,
while the discount parameter σ remains unchanged.

3.3.2 Moment-matching for the IBP In order to implement the moment-matching methodology we first need to evaluate the posterior moments of the random total mass. For this purpose, we rely on the moments characterization in terms of the cumulants provided in Proposition 1. The cumulants κ * i of the CRM part μ * (X) are obtained from Table 1 with the appropriate parameter updates which leads to

κ * i = a * (1 -σ) (i-1) (1 + c * ) (i-1) = a (1 -σ) (i-1) (c + σ) (n) (1 + c) (n+i-1) .
We consider two stable-beta processes: the beta process prior μ ∼ SBP(c = 1, σ = 0, a = 1) and the stable-beta process prior μ ∼ SBP(c = 1, σ = 0.5, a = 1). We let n vary in {5, 10, 20}. In contrast to the NRMI case, there is no need to work under different scenarios for the clustering profile of the data, since the posterior CRM μ * is not affected by them with only the sample size entering the updating scheme. We compare the prior moment-match for μ with the posterior moment-match for μ * in terms of our discrepancy index M and the results are displayed in Figure 5. The comparison shows that there is a gain in precision between prior and posterior distributions in terms of M suggesting that the a priori error level represents an upper bound for the posterior approximation error.

As in Section 3.2, we also evaluate the relative weights of fixed jumps and posterior CRM or, roughly, of the data w.r.t. the prior. Recalling that fixed location jumps J * j are independent and beta(n j -σ, c + σ + n -n j ) and some algebra allow to re-write the ratio of interest as

E k j=1 J * j E μ * (X) = (n -kσ)(c + 1) (n-1) a(c + σ) (n) .
Table 2 displays the corresponding values for different choices of n and k. As in the NRMI case, the fixed part overcomes the CRM part, which means that the data dominate the prior, and, moreover, their relative weight increases as n increases. In terms of moment-matching this shows that, if one looks at the overall posterior structure, the approximation error connected to the truncation is further dampened. .6 q q q q q q q q q q q q q q q q q q q q q q q q q q Beta : σ = 0, a = .6 q q q q q q q q q q q q q q q q q q q q q q q q q q SBP : σ = 0.

Practical use of the moment-matching criterion

We illustrate the use of the moment-matching strategy by implementing it within location-scale NRMI mixture models, which can be represented in hierarchical form as

Y i |µ i , σ i ind ∼ k(•|µ i , σ i ), i = 1, . . . , n, (µ i , σ i )| P i.i.d. ∼ P , i = 1, . . . , n, P ∼ NRMI,
where k is a kernel parametrized by (µ, σ) ∈ R × R + and the NRMI P is defined in (16). Under this framework, density estimation is carried out by evaluating the posterior predictive density. Specifically, we consider the Gaussian kernel k(x|µ, σ) = N (x|µ, σ) and NGG on locations and scales with a normal base measure P 0 , parameter θ = 1 in Equation ( 5), and varying stability parameter γ ∈ {0, 0.25, 0.5, 0.75}.

The dataset we consider is the popular Galaxy dataset, which consists of velocities of 82 distant galaxies diverging from our own galaxy. Since the data are M , Fe M ) between estimated cdfs F M and Fe M under, respectively, the moment-match (with M = 0.01) and the relative error (with eM = 0.1, 0.05, 0.01) criteria. The mixing measure of normal mixture is the normalized generalized gamma process with varying γ ∈ {0, 0.25, 0.5, 0.75}.

Lemma 1. Let function N ( • ) be as in (9) for the stable-beta process. Then for any ξ > 0

N -1 (ξ) ≤        e 1-ξ/a c if σ = 0, (αξ + β) -1/σ if σ ∈ (0, 1),
where α = σΓ(1 -σ) Γ(c+σ) aΓ( c+1) and β = 1 -σ c+σ Γ(1 -σ).

Proof. For σ = 0, from u -1 (1 -u) c-1 ≤ u -1 + (1 -u) c-1 one obtains

1 v u -1 (1 -u) c-1 du ≤ 1/c -log v.
Hence, N (v)/a ≤ 1 -c log v and N -1 (ξ) ≤ e (1-ξ/a)/c . The argument for σ = 0 follows along the same lines starting from u -1-σ (1 -u) σ+c-1 ≤ Γ(1 -σ)u -σ-1 + (1 -u) σ+c-1 .

Proposition 4. Let (ξ j ) j≥1 be the jump times for a homogeneous Poisson process on R + with unit intensity. Define the tail sum of the stable-beta process as

T M = ∞ j=M +1 N -1 (ξ j ),
where N ( • ) is given by (9). Then for any ∈ (0, 1),

P T M ≤ t M ≥ 1 -, for t M =        C 1 e 1 c -M C 1 if σ = 0, σ 1-σ (C 2 ) 1/σ (M +βC 2 / ) 1/σ-1 if σ ∈ (0, 1),
where C 1 = 2ace and C 2 = 2e/α do not depend on .

Proof. The proof follows along the same lines as the proof of Theorem A.1. in [START_REF] Brix | Generalized gamma measures and shot-noise Cox processes[END_REF]. Let q j denote the 2 M -j quantile, for j = M + 1, M + 2, . . ., of a gamma distribution with mean and variance equal to j. Then

P ∞ j=M +1 N -1 (ξ j ) ≤ ∞ j=M +1
N -1 (q j ) ≥ 1 -.

An upper bound on t M = ∞ j=M +1 N -1 (q j ) is then found by resorting to Lemma 1 along with the inequality q j ≥ 2e j. If σ = 0

t M ≤ e 1/c ∞ j=M +1 e - q j ac ≤ e 1/c ∞ j=M +1 e -j
2ace ≤ e 1/c 2ace e -M 2ace ,

whereas if σ = 0 t M ≤ ∞ j=M +1 (αq j +β) -1 σ ≤ ∞ j=M +1 α j 2e + β -1 σ = 2e α -1 σ ∞ j=M +1 j + 2eβ α -1 σ .
The result follows by bounding the last sum by

∞ M x + 2eβ α -1 σ dx.
The bound t M obtained in Proposition 4 is exponential when σ = 0 and polynomial when σ = 0, but it is very conservative as already pointed out by [START_REF] Brix | Generalized gamma measures and shot-noise Cox processes[END_REF]. This finding is further highlighted in the table associated to Figure 6, where the bound t M is computed with appropriate constants derived from the proof. In contrast, the bound t M obtained by direct calculation of the quantiles q j (instead of resorting to a lower bound on them) is much sharper. Figure 6 displays the sharper bound t M . Inspection of the plot demonstrates a decrease pattern in this bound in probability which is reminiscent of the ones for the indices M and e M studied in the paper. This observation is a further indication that the Ferguson & Klass algorithm is a tool with well-behaved approximation error. q q q q q q q q q q q q q q q q q q q q q q q q q q SBP : a = 1, c = 1 σ = 0 σ = 0. Stable-beta process with parameters σ = 0 and σ = 0.5. Left: Bound in probability t M of the tail sum T M obtained by direct calculation of the quantiles q j with = 10 -2 as the truncation level M increases. Right: Bounds t M (provided in Proposition 4) and t M (obtained by direct calculation of the quantiles q j ) of the tail sum after M jumps with = 10 -2 .

Fig 1 :Fig 2 :

 12 Fig 1: Left panel: M as M varies; right panel: e M as M varies. Top row: generalized gamma process (GG) with varying γ and a = 1 fixed; middle row: inverse-Gaussian process (IG), γ = 0.5, with varying total mass a; bottom row: stablebeta process (SBP) with a = 1, c = 0.5 fixed and varying discount parameter σ. The points are connected by straight lines only for visual simplification.

Fig 3 :

 3 Fig 3: NGG posterior: density f U |Y with n = 10 observations, a = 1, γ = 0.5, and number of clusters k ∈ {1, . . . , 10}; k = 1 corresponds to the most peaked density and k = 10 to the flattest.

Fig 4 :

 4 Fig 4: Inverse-Gaussian process (γ = 0.5) with a = 1. Left: Moment-matching errors M as the number of jumps M varies. M corresponding to prior μ (continuous line) and posterior μ * under Y clustering scenarios (i) (dashed line), (ii) (dotted line), (iii) (dotted-dashed line). Right: Index of relative importance E k j=1 J * j /E μ * (X) for varying (n, k).

Fig 5 :

 5 Fig 5: Moment-matching errors M as the number of jumps M varies for the stable-beta process with c = 1, a = 1, and, respectively, σ = 0 (left panel) and σ = 0.5 (right panel). M corresponding to prior μ (continuous line) and the posterior μ * given with n = 5 (dashed line) and n = 10 (dotted line) and n = 20 (dashed-dotted line) observations.

  Fig 6:Stable-beta process with parameters σ = 0 and σ = 0.5. Left: Bound in probability t M of the tail sum T M obtained by direct calculation of the quantiles q j with = 10 -2 as the truncation level M increases. Right: Bounds t M (provided in Proposition 4) and t M (obtained by direct calculation of the quantiles q j ) of the tail sum after M jumps with = 10 -2 .

  and ξ 1 , ξ 2 , . . . are jump times of a standard Poisson process (PP) of unit rate i.e. ξ 1 , ξ 2 -ξ 1 , . . .

i.i.d.

  with Q the law of NRMI, and set Y = (Y 1 , . . . , Y n ). Due to the discreteness of NRMIs, ties will appear with positive probability in Y and, therefore, the sample information can be encoded by the K n = k distinct observations (Y * 1 , . . . , Y * k ) with frequencies (n 1 , . . . , n k ) such that k j=1 n j = n. Moreover, introduce the nonnegative random variable U such that the distribution of [U |Y] has density, w.r.t. the Lebesgue measure, given by (

	(17)	i | P ∼ Q, P i.i.d. ∼ P , for i = 1, . . . , n,

  . Let (Y n ) n≥1 be as in (17) where P is an NRMI defined in (16) with Lévy intensity as in (4). Then the posterior distribution of the unnormalized CRM μ, given a sample Y, is a mixture of the distribution of [μ|U, Y] with respect to the distribution of [U |Y]. The latter is identified by (18), whereas [μ|U = u, Y] is equal in distribution to a CRM with fixed points of discontinuity at the distinct observations Y * j ,

		k
	(19)	μ * +	J * j δ Y * j
		j=1
	such that:
	(a) μ * is a CRM characterized by the Lévy intensity
	(20)	ν

* (dv, dx) = e -uv ν(dv, dx), (b) the jump height J * j corresponding to Y * j has density, w.r.t. the Lebesgue measure, given by

Table 3

 3 Galaxy dataset. Kolmogorov-Smirnov distance dKS( F

	γ	eM = 0.1 eM = 0.05 eM = 0.01
	0	19.4	15.5	9.2
	0.25	31.3	23.7	15.1
	0.5	42.4	28.9	18.3
	0.75	64.8	41.0	23.2

Research supported by the European Research Council (ERC) through StG "N-BNP" 306406.

clearly away from zero (range from 9.2 to 34), Gaussian kernels, although having the whole real line as support, are typically employed in its analysis.

As far as the simulation algorithm is concerned, based on Sections 3.1 to 3.3, the following moment-matching Ferguson & Klass posterior sampling strategy is implemented: (1) evaluate the threshold M ( ) which validates trajectories of the CRM using Algorithm 1 on the prior distribution; (2) implement Algorithm 1 on the posterior distribution using the threshold M ( ). More elaborate and suitably tailored moment-matching strategies can be devised for specific models. However, to showcase the generality and simplicity of our proposal we do not pursue this here.

In particular, we set M = 0.01. We compare the output to the Ferguson & Klass algorithm with heuristic relative error e M criterion, which consists of step (2) only with truncation dictated by the relative error for which we set e M ∈ {0.1, 0.05, 0.01}. For both algorithms the Gibbs sampler is run for 20, 000 iterations with a burn-in of 4, 000, thinned by a factor of 5.

In order to compare the results, we compute the Kolmogorov-Smirnov distance d KS ( F M , Fe M ) between associated estimated cumulative distribution functions (cdf) F

M and Fe M under, respectively, the moment-match and the relative error criteria. The results are displayed in Table 3. The estimated cdf F M with M = 0.01 can be seen as a reference estimate since the truncation error is controlled uniformly across the different values of γ by the moment-match at the CRM level. First, one immediately notes that the smaller e M , the closer the two estimates become (in the d KS distance). Second, and more importantly, the numerical values of the distances heavily depend on the particular choice of the parameter γ for any given e M . In fact, F

M and Fe M are significantly further apart for large values of γ than for small ones. This clearly shows that the quality of approximation with the heuristic criterion of the relative index is highly variable in terms of a single parameter; in passing from γ = 0 to γ = 0.75 the distance increases by at least a factor of 2. This means that for comparing correctly CRM based models with different parameters one would need to pick different relative indices for each value of the parameter. However, there is no way to guess such thresholds without the guidance of an analytic criterion. And, this already happens by varying a single parameter, let alone when changing CRMs for which the same e M could imply drastically different truncation errors. This seems quite convincing evidence supporting the abandonment of heuristic criteria for determining the truncation threshold and the adoption of principled approaches such as the moment-matching criterion proposed in this paper.

APPENDIX A: PROOF OF PROPOSITION 1

For any measurable set A of X, the n-th moment of μ(A), if it exists, is given by m

A (0) denotes the n-th derivative of the Laplace transform L A in (2) evaluated at 0. The result is proved by applying Faà di Bruno's formula to (2) for obtaining the derivatives.

APPENDIX B: EVALUATION OF THE TAIL SUM OF THE STABLE-BETA PROCESS

Here we provide an evaluation of the tail sum (11) in the case of the stable-beta process. We start by stating a lemma useful for upper bounding the tail sum.