Irina Kovalenko

Daniel Hestroffer

A Doressoundiram

Nicolai Emelyanov

Radu S Stoica

Statistical inversion method for binary asteroids orbit determination

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

We focus on the study of binary asteroids, which are common in the Solar system from its inner to its outer regions. These objects provide fundamental physical parameters such as mass and density, and hence clues on the early Solar System, or other processes that are affecting asteroid over time. The present method of orbit computation for resolved binaries is based on Markov Chain Monte-Carlo statistical inversion technique. Particularly, we use the Metropolis-Hasting algorithm with Thiele-Innes equation for sampling the orbital elements and system mass through the sampling of observations. The method requires a minimum of four observations, made at the same tangent plane; it is of particular interest for orbit determination over short arcs or with sparse data. The observations are sampled within their observational errors with an assumed distribution. The sampling yields the whole region of possible orbits, including the one that is most probable. • 𝑝 𝝋 𝑿 = 𝑝 𝜀 (observational error p. d. f. )

Statistical inversion method

= exp[- 1 2 𝝋 -𝜓 𝑿 𝑇 Λ -1 𝝋 -𝜓 𝑿 ]
N observations 𝑺 = 𝜌 1 , 𝜃 1 ; 𝜌 2 , 𝜃 2 ; 𝜌 3 , 𝜃 3 ; 𝜌 4 Select a random set of 4 observations

The Metropolis-Hastings algorithm will be used for sampling the parameters 𝑿.

For each t iteration: This process is repeated until the stationary a posteriori density is reached.

𝑺 ′ -
z, km

x, km y, km

  for binary asteroids orbit determination I.Kovalenko, D. Hestroffer, A. Doressoundiram, N. Emelyanov, R. Stoica ikovalenko@imcce.fr (IMCCE -LESIA, Paris Observatory) Orbit determination References 1) R. Palacios 1958 AJ 63, 395 2) D. Oszkiewicz et al. 2013 ,SF2A 237 The astrometric observations are related to the theoretical positions through the observational equation: 𝝋 = 𝜓 𝑿 + 𝜀 • Observations: 𝝋 = 𝜌 1 , 𝜃 1 ; … ; 𝜌 N , 𝜃 N • Sky-plane position: 𝜓 𝑿 • Orbital elements + system's mass: 𝑿 = (𝑎, 𝑒, 𝑖 , Ω, 𝜔, 𝑀, 𝑚 𝑠𝑦𝑠 ) • Observational errors: 𝜀 = (𝜀 𝛼1 , 𝜀 𝛿1 ; … ; 𝜀 𝛼𝑁 , 𝜀 𝛿𝑁 ) 𝑝 𝑋 𝜑 = 𝑝 𝜑 𝑋 𝑝 𝑋 𝑝(𝜑) Markov Chain Monte-Carlo method 𝑝 𝑋 𝜑 ∝ 𝑝 𝜑 𝑋 𝑝 𝑋 Where • 𝑝(𝑋) ∝ 𝑑𝑒𝑡Λ -1

[ 2 ]

 2 𝐽 ′ and 𝐽 𝑡 -the determinants of Jacobians from coordinates to orbital parameters 𝐽 = det 𝜕𝑆 𝜕𝑋 𝑿 𝑡+1 = 𝑿′ 𝑿 𝒕+𝟏 = 𝑿 ′ , with probability 𝑎 or 𝑿 𝒕+𝟏 = 𝑿 𝒕 ,with probability 1 -𝑎 Thiele -Innes method [1]: -𝐸 𝑗 -sin (𝐸 𝑖 -𝐸 𝑗 )] 𝑖, 𝑗 = 1, 2, 3,4 Δ𝑖𝑗 = 𝜌 𝑖 𝜌 𝑗 sin (𝜃 𝑗 -𝜃 𝑖 ) => 𝑿 = (𝑎, 𝑒 𝑖 , 𝜔, Ω, 𝑀, 𝑇) 𝐼𝑓 𝑎 ≥ 1 𝐼𝑓 𝑎 < 1

  proposal set of observations 𝑺 𝑡 -last accepted set of observations

	p 𝑡 (𝑺 ′ ; 𝑆 𝑡 )-proposal density ′ ∝ 𝐺 𝛼 𝑖 𝛼 𝑖 𝑡 , 𝜎(𝛼 𝑖 ) , 𝛿 𝑖 ′ ∝ 𝐺 𝛿 𝑖 𝑡 , 𝜎(𝛿 𝑖 )
	Orbit 𝑿 ′ + monitor the fit to all
	observations
	𝐴𝑐𝑐𝑒𝑝𝑡𝑒𝑛𝑐𝑒 𝑐𝑟𝑖𝑡𝑒𝑟𝑖𝑎
	𝑎 =	𝑝 𝑋 ′ |𝜑 |𝐽 𝑡 | 𝑝 𝑋 t |𝜑 |𝐽 ′ |

Acknowledgements

This work is supported by the Labex ESEP (ANR N° 2011-LABX-030).

The algorithm is run for a large number of iterations until the entire possible orbital-element space is mapped


Alireza Roodaki

Julien Bect

Gilles Fleury

Relabeling and Summarizing Posterior Distributions in Signal Decomposition Problems when the Number of Components is Unknown

Keywords: Bayesian inference, Signal decomposition, Trans-dimensional MCMC, Label-switching, Stochastic EM

This paper addresses the problems of relabeling and summarizing posterior distributions that typically arise, in a Bayesian framework, when dealing with signal decomposition problems with an unknown number of components. Such posterior distributions are defined over union of subspaces of differing dimensionality and can be sampled from using modern Monte Carlo techniques, for instance the increasingly popular RJ-MCMC method. No generic approach is available, however, to summarize the resulting variable-dimensional samples and extract from them component-specific parameters.

We propose a novel approach, named Variable-dimensional Approximate Posterior for Relabeling and Summarizing (VAPoRS), to this problem, which consists in approximating the posterior distribution of interest by a "simple"-but still variabledimensional-parametric distribution. The distance between the two distributions is measured using the Kullback-Leibler divergence, and a Stochastic EM-type algorithm, driven by the RJ-MCMC sampler, is proposed to estimate the parameters. Two signal decomposition problems are considered, to show the capability of VAPoRS both for relabeling and for summarizing variable dimensional posterior distributions: the classical problem of detecting and estimating sinusoids in white Gaussian noise on the one hand, and a particle counting problem motivated by the Pierre Auger project in astrophysics on the other hand.

I. INTRODUCTION

Nowadays, owing to the advent of Markov Chain Monte Carlo (MCMC) sampling methods [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF][START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF][START_REF] Liu | Monte Carlo Strategies in Scientific Computing[END_REF][START_REF] Robert | Monte Carlo Statistical Methods[END_REF], Bayesian data analysis is considered as a conventional approach in machine learning, signal and image processing, and data mining problemsto name but a few. Nevertheless, in many applications, practical challenges remain in the process of extracting, from the generated samples, quantities of interest to summarize the posterior distribution.

Summarization consists, loosely speaking, in providing a few simple yet interpretable parameters and/or graphics to Copyright (c) 2014 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works.

Alireza Roodaki did this work during his Ph.D. in SUPELEC (1) and during his post-doctoral stay in LTCI (2) . The results presented here are also part of his Ph.D. dissertation [START_REF] Roodaki | Signal decompositions using trans-dimensional Bayesian methods[END_REF]. Email: al.roodaki@gmail.com.

Julien Bect is Associate Professor in SUPELEC (1) , Department of Signal Processing and Electronic System. Gilles Fleury is Head of Research in SUPELEC (1) . Email: firstname.lastname@supelec.fr.

(1) E3S-SUPELEC Systems Sciences, SUPELEC, Gif-sur-Yvette, France.

(2) LTCI, CNRS / Télécom ParisTech, Paris, France.

the end-user of a statistical method. For instance, in the case of a scalar parameter with a unimodal posterior distribution, measures of location and dispersion (e.g., the empirical mean and the standard deviation, or the median and the interquartile range) are typically provided in addition to a graphical summary of the distribution (e.g., a histogram or a kernel density estimate). In the case of multimodal distributions, summarization becomes more difficult but can be carried out using, for instance, the approximation of the posterior by a Gaussian Mixture Model (GMM) [START_REF] West | Approximating posterior distributions by mixture[END_REF]. Summarizing or approximating posterior distributions has also been used in designing proposal distributions of Metropolis-Hastings (MH) samplers in an adaptive MCMC framework; see, e.g., [START_REF] Haario | An adaptive metropolis algorithm[END_REF][START_REF] Bai | Divide and conquer: a mixturebased approach to regional adaptation for MCMC[END_REF][START_REF] Bardenet | An adaptive Metropolis algorithm with online relabeling[END_REF]. This paper addresses the problem of summarizing posterior distributions in the case of some trans-dimensional problems (i.e., "problems in which the number of things that we don't know is one of the things that we don't know" [START_REF] Green | Reversible jump MCMC computation and Bayesian model determination[END_REF][START_REF] Green | Trans-dimensional Markov chain Monte Carlo[END_REF]). More specifically, we concentrate on the problem of signal decomposition when the number of components is unknown, which is an important case of trans-dimensional problem. Examples of such problems include the detection and estimation of sinusoids in white Gaussian noise [START_REF] Andrieu | Joint Bayesian model selection and estimation of noisy sinusoids via reversible jump MCMC[END_REF] and the related problem of estimating directions of arrival in array processing [START_REF] Larocque | Reversible jump MCMC for joint detection and estimation of sources in coloured noise[END_REF], the detection of objects in images [START_REF] Rue | Bayesian object identification[END_REF][START_REF] Ortner | Building outline extraction from digital elevation models using marked point processes[END_REF], and the detection of physical particles (neutrons, muons, . . . ) using noisy data from various types of sensors, for instance in spectroscopy [START_REF] Andrieu | Bayesian deconvolution of noisy filtered point processes[END_REF] or astrophysics [START_REF]The Pierre Auger Project Design Report[END_REF][START_REF]Properties and performance of the prototype instrument for the Pierre Auger Observatory[END_REF].

Let y = (y 1 , y 2 , . . . , y N ) t be a vector of N observations, where the superscript t stands for vector transposition. As a generic description of a signal decomposition problem, we consider a countable family of models (M k ) k∈N , where it is assumed that, under model M k , the observed signal y is made of k components. A "component" might be a sinusoid (see Section I-B) or a decaying exponential (see Section III-B) in a one-dimensional signal processing problem, for instance, or an elementary geometric form in an image processing problem [START_REF] Ortner | Building outline extraction from digital elevation models using marked point processes[END_REF]; we simply assume that each component is completely described by a vector of parameters θ θ θ j ∈ Θ ⊆ R d , 1 ≤ j ≤ k. We denote by θ θ θ 1:k = (θ θ θ 1 , . . . , θ θ θ k ) ∈ Θ k the vector of all component-specific parameters, where Θ 0 = {∅}.

One feature that the problems we are considering have in common is the invariance of the likelihood p (y | k, θ θ θ 1:k ) with respect to permutations (relabeling) of the components, which is called the "label-switching" issue in the literature; see, e.g., [START_REF] Richardson | On Bayesian analysis of mixtures with an unknown number of components[END_REF][START_REF] Stephens | Dealing with label switching in mixture models[END_REF][START_REF] Jasra | Markov chain Monte Carlo methods and the label switching problem in Bayesian mixture modeling[END_REF][START_REF] Celeux | Computational and inferential difficulties with mixture posterior distributions[END_REF][START_REF] Frühwirth-Schnatter | Dealing with label switching under model uncertainty[END_REF]. We will discuss this issue further in Section I-A.

In a Bayesian framework, a joint posterior density f (k, θ θ θ 1:k ) p (k, θ θ θ 1:k | y) is obtained through Bayes' formula for the number k of components and the vector of componentspecific parameters, after assigning prior distributions on them:

f (k, θ θ θ 1:k ) ∝ p (y | k, θ θ θ 1:k ) p (θ θ θ 1:k | k) p (k) , (1)
where ∝ indicates proportionality. This density is defined over a variable-dimensional space Θ Θ Θ, which is a union of subspaces of differing dimensionality, i.e., Θ Θ Θ = ∪ k∈N {k} × Θ k . The posterior density [START_REF] Roodaki | Signal decompositions using trans-dimensional Bayesian methods[END_REF] completely describes the information (and the associated uncertainty) provided by the data y about the candidate models and the vector of unknown parameters. Since it is only known up to a normalizing constant in most cases, and potentially multimodal, Monte Carlo simulation methods, such as Reversible Jump MCMC (RJ-MCMC) [START_REF] Green | Reversible jump MCMC computation and Bayesian model determination[END_REF], have been widely used to approximate it.

A. The label-switching issue

One of the most challenging issues when attempting at summarizing posterior distributions, that even occurs in fixeddimensional situations, is the label-switching phenomenon (see, e.g., [START_REF] Richardson | On Bayesian analysis of mixtures with an unknown number of components[END_REF][START_REF] Stephens | Dealing with label switching in mixture models[END_REF][START_REF] Jasra | Markov chain Monte Carlo methods and the label switching problem in Bayesian mixture modeling[END_REF][START_REF] Celeux | Computational and inferential difficulties with mixture posterior distributions[END_REF][START_REF] Frühwirth-Schnatter | Dealing with label switching under model uncertainty[END_REF][START_REF] Papastamoulis | An artificial allocations based solution to the label switching problem in Bayesian analysis of mixtures of distributions[END_REF][START_REF] Rodríguez | Label switching in bayesian mixture models: Deterministic relabeling strategies[END_REF]), which is caused by the invariance of both the likelihood and the prior distribution under permutations of the components. As a consequence, the componentspecific marginal posterior distributions are all equal, and therefore useless for the purpose of summarizing the information contained in the posterior distribution about individual components. A symptom of this issue is the multimodality of marginal posterior distributions.

The simplest way of dealing with the label-switching issue is to introduce an Identifiability Constraint (IC), such as sorting the components with respect to one of their parameters; see [START_REF] Richardson | On Bayesian analysis of mixtures with an unknown number of components[END_REF] for more discussion concerning the use of ICs in the Bayesian analysis of GMMs. However, in most practical examples, choosing an appropriate IC manually is not feasible. Many relabeling algorithms have therefore been developed to "undo" the label-switching effect automatically-i.e., change sample labels to make the marginals as unimodal as possiblebut all of them are restricted to the case of fixed-dimensional posterior distributions; see [START_REF] Frühwirth-Schnatter | Dealing with label switching under model uncertainty[END_REF][START_REF] Rodríguez | Label switching in bayesian mixture models: Deterministic relabeling strategies[END_REF][START_REF] Sperrin | Probabilistic relabelling strategies for the label switching problem in bayesian mixture models[END_REF][START_REF] Yao | Model based labeling for mixture models[END_REF] for recent advances and references.

In variable-dimensional posterior distributions, there is an extra uncertainty about the "presence" of components, as will be clarified in the following illustrative example. This additional difficulty has hindered previous attempts to undo label-switching in the variable-dimensional scenario, where, according to [START_REF] Robert | Discussion of "On Bayesian analysis of mixtures with an unknown number of components[END_REF] "the meaning of individual components is vacuous".

B. Illustrative example: joint Bayesian detection and estimation of sinusoids in white Gaussian noise

In this example, it is assumed that under M k , the observed signal y is composed of k sinusoidal components observed in white Gaussian noise. That is, under M k , where a c, j and a s, j are the cosine and sine amplitudes of the j th sinusoidal component and ω j its radial frequency. Moreover, n is a white Gaussian noise of variance σ 2 .

y[i] = k ∑ j=1 (a c, j cos(ω j i) + a s, j sin(ω j i)) + n[i],
The unknown parameters are the number k of sinusoidal components, the vectors θ θ θ j = (a c, j , a s, j , ω j ) of component- specific parameters, 1 ≤ j ≤ k, and the noise variance σ 2 . Thus,

Θ = R 2 × (0, π) and Θ Θ Θ = ∪ k∈N {k} × Θ k ∪ R + .
We use the hierarchical model, prior distributions, and the RJ-MCMC sampler proposed in [START_REF] Andrieu | Joint Bayesian model selection and estimation of noisy sinusoids via reversible jump MCMC[END_REF] for this problem; the interested reader is thus referred to [START_REF] Green | Reversible jump MCMC computation and Bayesian model determination[END_REF][START_REF] Andrieu | Joint Bayesian model selection and estimation of noisy sinusoids via reversible jump MCMC[END_REF] for more details 1 .

Here, we consider an experiment in which the observed signal y of length N = 64 consists of three sinusoids with energies A A A 1:k = (20, 6.32, 20) t , where A j = a 2 c, j + a 2 s, j , phases φ φ φ 1:k = (0, π/4, π/3) t , where φ j = -arctan(a s, j /a c, j ), and true radial frequencies ω ω ω 1:k = (0.63, 0.68, 0.73) t . The signal-to- noise ratio 2 is set to the moderate value of 7 dB. Figure 1 represents the posterior distributions of both the number k of components and the sorted radial frequencies ω ω ω 1:k = (ω 1 , . . . , ω k ) t given k obtained using 100 000 samples generated by the RJ-MCMC sampler. Note that, here, we used sorting to mitigate the effect of label-switching for visualization. Each row is dedicated to one value of k, for 2 ≤ k ≤ 4. Observe that other models have negligible posterior probabilities, since

p(2 ≤ k ≤ 4 | y) = 0.981.
Roughly speaking, two approaches co-exist in the literature for summarizing variable-dimensional posterior distributions: Bayesian Model Selection (BMS) and Bayesian Model Averaging (BMA). The BMS approach ranks models according to their posterior probabilities p(k | y), selects one model, denoted by k MAP here, where MAP stands for Maximum A Posteriori, and then summarizes the posterior distribution of the component-specific parameters under the (fixed-dimensional) 1 In fact, the "Birth-or-Death" moves' acceptance ratio provided in the seminal paper [START_REF] Andrieu | Joint Bayesian model selection and estimation of noisy sinusoids via reversible jump MCMC[END_REF] is erroneous. See [START_REF] Roodaki | Signal decompositions using trans-dimensional Bayesian methods[END_REF]Chapter 1] or [START_REF] Roodaki | Comments on "Joint Bayesian model selection and estimation of noisy sinusoids via reversible jump MCMC[END_REF] for justification and true expression of the acceptance ratio, which is used in this paper.

2 defined here as D a 1:k 2 / Nσ 2 , where D is the N × 2k "design matrix" of sines and cosines associated to ω ω ω 1:k and a

1:k = a c,1 , a s,1 , . . . , a c,k , a s,k t .
selected model. This is at the price of losing valuable information provided by the other (discarded) models. For instance, in the example of Figure 1, all information about the smalland therefore harder to detect-middle component is lost by selecting the most a posteriori probable model M 2 . On the other hand, the BMA approach consists in reporting results that are averaged over all possible models. Although BMA is suitable for signal reconstruction and prediction, it is not appropriate for studying component-specific parameters, the number of which changes in each model 3 . More information concerning these two approaches can be found in [START_REF] Green | Reversible jump MCMC computation and Bayesian model determination[END_REF][START_REF] Clyde | Model uncertainty[END_REF] and references therein.

To the best of our knowledge, no generic method is currently available that would allow to summarize the information that is so easily read on Figure 1 for this very simple example 4 : namely, that there seem to be three sinusoidal components in the observed noisy signal, the middle one having a smaller "probability of presence" than the others.

C. Outline of the paper

In this paper, we propose a novel approach, named Variabledimensional Approximate Posterior for Relabeling and Summarizing (VAPoRS), for relabeling and summarizing posterior distributions defined over variable-dimensional subspaces that typically arise in signal decomposition problems when the number of components is unknown. It consists in approximating the true posterior distribution with a parametric model (of varying-dimensionality), by minimization of the Kullback-Leibler (KL) divergence between the two distributions. A Stochastic Expectation Maximization (SEM)-type algorithm [START_REF] Celeux | The SEM algorithm : a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem[END_REF][START_REF] Celeux | A stochastic approximation type EM algorithm for the mixture problem[END_REF][START_REF] Nielsen | The stochastic EM algorithm: estimation and asymptotic results[END_REF], driven by the output of an RJ-MCMC sampler, is used to estimate the parameters of the approximate model.

VAPoRS shares some similarities with the relabeling algorithms proposed in [START_REF] Stephens | Dealing with label switching in mixture models[END_REF][START_REF] Sperrin | Probabilistic relabelling strategies for the label switching problem in bayesian mixture models[END_REF][START_REF] Yao | Model based labeling for mixture models[END_REF] to solve the label switching problem, and also with the EM-type algorithm used in [START_REF] Bai | Divide and conquer: a mixturebased approach to regional adaptation for MCMC[END_REF] in the context of adaptive MCMC algorithms (both in a fixeddimensional setting). The main contribution of this paper is the introduction of an original variable-dimensional parametric model, which allows to tackle directly the difficult problem of approximating a distribution defined over a union of subspaces of differing dimensionality, and thus provides a first solution to the "trans-dimensional label-switching" problem, so to speak.

Perhaps, the algorithm that we propose can be seen as a realization of the idea that M. Stephens had in mind when he wrote [34, page 94]:

"This raises the question of whether we might be able to obtain an alternative view of the [variable-dimensional] posterior by combining the results for all different k's, and grouping together components which are "similar", in that they have similar predictive density estimates. However, attempts to do this have failed to produce any easily interpretable results."

The paper is organized as follows. Section II introduces the proposed model and stochastic algorithm for relabeling and summarizing variable dimensional posterior distributions. Section III illustrates the performance of VAPoRS using two signal decomposition examples, namely, the problem of joint Bayesian detection and estimation of sinusoids in white Gaussian noise and the problem of joint Bayesian detection and estimation of particles in the Auger project (in astrophysics). Section IV confirms the performances of VAPoRS using a Monte Carlo experiment. Finally, Section V concludes the paper and gives directions for future work.

II. VAPORS

A. Introduction: observed components and t-components

We assume that the target posterior distribution, defined on the variable-dimensional space Θ Θ Θ = k∈N {k} × Θ k , admits a probability density function (pdf) f with respect to the kddimensional Lebesgue measure on each {k} × Θ k , k ∈ N * .

Our objective is to approximate the true posterior density f using a "simple" parametric model. This parametric model will also be defined on the variable-dimensional space Θ Θ Θ (i.e., it is not a fixed-dimensional approximation as in BMS).

We assume that a Monte Carlo sampling method-e.g., an RJ-MCMC sampler [START_REF] Green | Reversible jump MCMC computation and Bayesian model determination[END_REF]-is available to generate M samples from f , which we denote by θ θ θ (i) = k (i) , θ θ θ (i) 1:k (i) , for i = 1, . . . , M. These M variable-dimensional samples, which will be used to fit the approximate model to the true posterior, will play the role of input data for our method; therefore, an element θ θ θ

(i) j ∈ Θ (1 ≤ j ≤ k (i)
) of the sample θ θ θ (i) will be referred to as an "observed component" (remember that a "component" is a point in the set Θ).

As our main device to handle the variable dimensionality of Θ Θ Θ, our parametric model (described more precisely in the next section) introduces "transdimensional components", that we will refer to as t-components for the sake of brevity. Intuitively, in the example shown on Figure 1, these t-components will serve the purpose of aggregating bumps in parameter posterior probability across dimensions, thus identifying observed components occurring in different models as manifestations of the same true component in the signal.

B. Variable-dimensional parametric model

Instead of trying to describe the proposed parametric model directly, let us now adopt a generative point of view, i.e., let us describe how to sample an Θ Θ Θ-valued random variable θ θ θ = (k, θ θ θ 1:k ) from the corresponding probability distribution. We assume that a positive integer L is given, which represents the number of t-components in the model. Each t-component can be thought of as a "virtual component" that can generate zero or one (or several, in the case of the Poisson point process component that will be introduced later) observed component(s), according to some prescribed distribution on Θ.

To generate a sample θ θ θ ∈ Θ Θ Θ, we first generate, independently for each of the L t-components, a binary indicator variable ξ l ∈ {0, 1} drawn from the Bernoulli distribution Ber(π l ), where ξ l = 1 indicates that an observed component corresponding to the l th t-component is actually present in θ θ θ . The actual number k of components in the generated samples To generate a sample from the VAPoRS model, do:

1) For l = 1, . . . , L, generate ξ l ∼ Ber(π l ). 2) Set k = ∑ L l=1 ξ l . 3) For each l such that ξ l = 1, generate θ θ θ l ∼ N (µ µ µ l , Σ Σ Σ l ).

4) Randomly arrange the samples θ

θ θ l generated at step 3 in a vector θ θ θ 1:k = (θ θ θ 1 , . . . , θ θ θ k ). is thus k = ∑ L l=1 ξ l . The parameter π l ∈ (0, 1] will be called the probability of presence of the l th component.

Second, given the vector of indicator variables ξ ξ ξ = (ξ 1 , . . . , ξ L ), a Θ-valued random vector is generated for each t-component that is present (i.e., for each l such that ξ l = 1). This random vector is generated according to some probability distribution associated to the t-component, that will be assumed, in this paper 5 , to be a d-dimensional Gaussian distribution with mean µ µ µ l and covariance matrix Σ Σ Σ l . In order to achieve the required invariance with respect to component relabeling, the generated vectors are randomly arranged in a vector θ θ θ 1:k = (θ θ θ 1 , . . . , θ θ θ k )-i.e., a permutation of the k components that are present is drawn uniformly in the set of all permutations. The above mentioned steps to generate samples from the proposed parametric model are summarized in Figure 2.

Example. Assume that there are L = 3 univariate Gaussian tcomponents in the model, with means µ µ µ = (0.63, 0.68, 0.73), variances s s s 2 = (0.01, 0.02, 0.01) and probabilities of presence π π π = (0.8, 0.3, 0.8). Figure 3 shows the pdf's of the three Gaussian t-components along with six random samples generated from the parametric model. Moreover, the kernel density estimates of 10 000 random samples generated from the parametric model are depicted in Figure 4 (a). It can be seen from both figures that the dimension of the generated samples varies from k = 0 to k = L = 3.

Let us look at the marginal posterior distributions of the sorted radial frequencies, depicted in the right panel of Figure 1. An important fraction of the posterior distribution is concentrated in compact bumps6 ; as will be seen later (see Section III-A), our Gaussian t-component will be effective at describing this fraction of the posterior. However, it can be observed on the plots related to the models with three and four sinusoidal components that the posterior also contains "diffuse" parts, meaning that a small fraction of the total probability mass is spread across a wide region (resulting, after sorting at fixed k, in heavy asymmetric tails for some marginal distributions). Another small fraction of the posterior distribution, not shown on Figure 1, is similarly scattered over the frequency axis for k = 4 and k = 5. It is clear that a model made of Gaussian t-components only will not be able to provide a parsimonious representation of this feature of the posterior f . Moreover, the observed components from the "diffuse" part, which would behave as outliers with respect to a model built with a small number of Gaussian tcomponents, could adversely influence the process of fitting the approximate posterior to the true posterior distribution of interest.

To overcome this issue, we propose to include in the model a "noise-like" Poisson Point Process (PPP; see, e.g., [START_REF] Karr | Point Processes and their Statistical Inference[END_REF]) to account for the presence of outliers in the observed samples. For simplicity 7 , we assume in this paper that Θ is bounded and that the PPP is homogeneous on Θ, with intensity λ /|Θ|. We denote by ξ L+1 ∈ N the number of components generated by the PPP, which follows a Poisson distribution with mean λ , and keep using the notation ξ ξ ξ for the extended vector (ξ 1 , . . . , ξ L , ξ L+1 ). Finally, ξ L+1 random samples are generated uniformly on Θ and inserted randomly among the samples drawn from the Gaussian components. The actual number k of components in the generated sample is thus, now, redefined as k = ∑ L+1 l=1 ξ l . Figure 5 provides the directed acyclic graph of the model.

Example (continued).

Setting Θ to the interval (0, π) (as is the case in the sinusoid detection problem) and λ = 0.5, Figure 4 (b) shows the intensities of generated samples from the toy example's parametric model equipped with the Poisson point process component. It can be observed that the model equipped with PPP is capable of generating diffuse samples and thus, provides a better approximation to the distribution of the observed samples in practice (see, for example, Figure 1). Another interesting point that can be seen in Figure 4 (b) is that the model with PPP is able to generate samples with dimensions greater than the given number L of Gaussian tcomponents. This allows the model to deal with vector of observed samples of dimension greater than L.

Notations. We denote by η η η l = (µ µ µ l , Σ Σ Σ l , π l ) the vector of parameters for the l th Gaussian t-component, and by η η η = (η η η 1 , . . . , η η η L , λ ) the full vector of parameters of the model. We denote by q η η η the pdf of the random variable θ θ θ = (k, θ θ θ 1:k ) generated by the above construction. As a convenient abuse of notation, we will also use q η η η in Sections II-C and II-D, to denote all joint, marginal and conditional distributions involving θ θ θ and the auxiliary variables used for its generation. θ θ θ

ξ l µ µ µ l Σ Σ Σ l π l ξ L+1 λ l = 1, 2, • • • , L

C. Distribution of the labeled samples

A random variable θ θ θ = (k, (θ θ θ 1 , . . . , θ θ θ k )) drawn from the density q η η η can be thought of as an "unlabeled sample", since the label l ∈ L {1, . . . , L + 1} of the component from which each θ θ θ j (1 ≤ j ≤ k) originates cannot be recovered from θ θ θ itself. Let us now introduce the (variable-dimensional) allocation vector

z = (k, (z 1 , . . . , z k )) ∈ k∈N {k} × L k ,
which provides the missing piece of information: z j = l indicates that θ θ θ j originates from the l th (Gaussian) t-component if l ≤ L, while z j = L + 1 indicates that θ θ θ j originates from the point process t-component. We will refer to the pair (θ θ θ , z) as a labeled sample. In the following, we will derive its joint

distribution q η η η (θ θ θ , z) = q η η η (θ θ θ | z)q η η η (z).
The distribution of the allocation vector z is

q η η η (z) = q η η η (z | ξ ξ ξ ) q η η η (ξ ξ ξ ), (2)
where q η η η (ξ ξ ξ ) is given by

q η η η (ξ ξ ξ ) = e -λ • λ ξ L+1 ξ L+1 ! L ∏ l=1 π ξ l l (1 -π l ) (1-ξ l ) . (3)
Note that ξ ξ ξ is a deterministic function of z:

ξ ξ ξ = n(z), with n l (z) = ∑ k j=1 1 z j =l , for 1 ≤ l ≤ L + 1.
To compute the first term of (2), remember that the points generated by the components of the parametric model are randomly arranged in θ θ θ 1:k . There- fore, for all ξ ξ ξ ∈ {0,

1} L × N such that ∑ L+1 l=1 ξ l = k, q η η η (z | ξ ξ ξ ) = ξ L+1 ! k! 1 ξ ξ ξ =n(z) , (4)
since two arrangements that differ only by the position of the points corresponding to the PPP give rise to the same allocation vector.

The conditional distribution q η η η (θ θ θ | z) reads

q η η η (θ θ θ | z) = k ∏ j=1 q η η η (θ θ θ j | z j ), (5)
where

q η η η (θ θ θ j | z j ) =    N θ θ θ j | µ µ µ z j , Σ Σ Σ z j if z j ≤ L, 1 |Θ| if z j = L + 1. (6)
Therefore, from Equations (3) to (6), we have

q η η η (θ θ θ , z) = e -λ k! λ |Θ| ξ L+1 ∏ 1≤ j≤k z j =L+1 N θ θ θ j | µ µ µ z j , Σ Σ Σ z j × L ∏ l=1 π ξ l l (1 -π l ) (1-ξ l ) 1 Z (z) , (7)
where (ξ 1 , . . . , ξ L+1 ) = n(z) and Z is the set of all allocation vectors (i.e., the set of all z ∈ ∪ k∈N {k} × L k such that ξ l = n l (z) ∈ {0, 1}, for 1 ≤ l ≤ L).

D. Estimating the model parameters

We propose to fit the parametric distribution q η η η to the posterior f of interest by minimizing a divergence measure 8 from f to q η η η . We use the KL divergence as a divergence measure in this paper, though other divergence measures can be used as well (see, e.g., [START_REF] Roodaki | Signal decompositions using trans-dimensional Bayesian methods[END_REF]Chapter 2]).

Denoting the KL divergence from f to q η η η by D KL ( f (θ θ θ ) q η η η (θ θ θ )), we define the criterion to be minimized as

J (η η η) D KL ( f (θ θ θ ) q η η η (θ θ θ )) = Θ Θ Θ f (θ θ θ ) log f (θ θ θ ) q η η η (θ θ θ ) dθ θ θ .
8 It would also be possible to assign prior distributions over the unknown parameters η η η and study their posterior distributions (for example, using an MCMC sampler with the latent variable z added to the state of the chain, in the spirit of the "data augmentation" algorithm [START_REF] Tanner | The calculation of posterior distributions by data augmentation[END_REF]). This would, however, leave the label-switching issue unsolved (because of the invariance of q η η η to permutations of its components).

At the (r + 1) th iteration, do:

(S-step) Draw allocation vectors z (i,r+1) , 1 ≤ i ≤ M, using an IMH step with target q η η η (r) (• • • | θ θ θ (i) ). (E-step) Construct the pseudo-completed log-likelihood

J M (η η η) = -∑ M i=1 log q η η η (θ θ θ (i) , z (i,r+1) ) .
(M-step) Estimate η η η (r+1) such that η η η (r+1) = argmin η η η J M (η η η). Using samples θ θ θ (i) , i = 1, . . . , M, generated by the RJ-MCMC sampler, this criterion can be approximated as

J (η η η) ≃ - 1 M M ∑ i=1 log q η η η (θ θ θ (i) ) + C, ( 8
)
where C = f (θ θ θ ) log f (θ θ θ ) dθ θ θ is a constant that does not depend on η η η. One should note that minimizing the right-hand side of (8) amounts to choosing

η η η = argmax η η η M ∑ i=1 log q η η η (θ θ θ (i) ) . (9)
To estimate the model parameters η η η, one of the extensively used algorithms for Maximum Likelihood (ML) parameter estimation in latent variable models is the EM algorithm proposed by [START_REF] Dempster | Maximum Likelihood from Incomplete Data via the EM Algorithm[END_REF]. However, it turns out that the EM algorithm, which has been used in similar works [START_REF] Bai | Divide and conquer: a mixturebased approach to regional adaptation for MCMC[END_REF][START_REF] Stephens | Dealing with label switching in mixture models[END_REF][START_REF] Sperrin | Probabilistic relabelling strategies for the label switching problem in bayesian mixture models[END_REF], is not appropriate for solving this problem, as computing the expectation in the E-step is intricate. More explicitly, in our problem the computational burden of the summation in the E-step over the set of all possible allocation vectors z increases very rapidly with both L and k. In fact, even for moderate values of L and k, say, L = 15 and k = 10, the summation is far too expensive to compute as it involves ∑ k m=0 L! (L-k+m)! ≈ 1.3 10 10 terms. In this paper, we propose to use the SEM algorithm [START_REF] Celeux | The SEM algorithm : a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem[END_REF][START_REF] Celeux | A stochastic approximation type EM algorithm for the mixture problem[END_REF][START_REF] Nielsen | The stochastic EM algorithm: estimation and asymptotic results[END_REF], a variation of the EM algorithm in which the E-step is substituted with stochastic simulation of the latent variables from their conditional posterior distributions given the previous estimates of the unknown parameters. In other words, at the iteration r + 1 of the SEM algorithm, denoting the estimated parameters at iteration r by η η η (r) , for i = 1, . . . , M, the allocation vectors z (i) are drawn from

q η η η (r) (• • • | θ θ θ (i)
). This step is called the Stochastic (S)-step. Then, these random samples are used to construct the so-called pseudo-completed log-likelihood.

Exact sampling from q η η η (r) (• • • | θ θ θ (i) ), as required by the S-step of the SEM-type algorithm, is unfortunately not feasible-not even using the accept-reject algorithm, due to the heavily combinatorial expression of the normalizing constant q η η η (r) (θ θ θ (i) ). Instead, since

q η η η (r) (z (i) | θ θ θ (i) ) ∝ q η η η (r) (θ θ θ (i) , z (i) )
can be computed up to a normalizing constant, we choose to use an Independent Metropolis-Hasting (IMH) step with q η η η (r) (z (i) | θ θ θ (i) ) as its stationary distribution; see [1, Algo- rithm 2.2] for more details.

The proposed SEM-type algorithm is summarized in Figure 6. Explicit expressions for the M-step are easily obtained from Equation [START_REF] Haario | An adaptive metropolis algorithm[END_REF] (see [1, page 76]). The computational cost of one iteration is of the order O ML 2 for a given d, the most expensive part being the IMH algorithm used in the S-step.

Remark 1. Convergence results of the SEM algorithm in the general form are provided by [START_REF] Nielsen | The stochastic EM algorithm: estimation and asymptotic results[END_REF] and, in the particular example of mixture analysis problems, by [START_REF] Diebolt | Asymptotic properties of a stochastic EM algorithm for estimating mixing proportions[END_REF]. Unfortunately, the assumptions in [START_REF] Nielsen | The stochastic EM algorithm: estimation and asymptotic results[END_REF][START_REF] Diebolt | Asymptotic properties of a stochastic EM algorithm for estimating mixing proportions[END_REF] do not hold in the problem we are dealing with as, 1) the observed samples θ θ θ (i) are correlated, owing to the fact that they are generated from the true posterior distribution using some MCMC methods, e.g., the RJ-MCMC sampler; 2) an IMH sampler is used to draw z (i) from the conditional posterior distribution. Nevertheless, empirical evidence of the "good" convergence properties of the SEMtype algorithm we proposed will be provided in the next two sections.

E. Robustified algorithm

Preliminary experiments with the SEM-type algorithm described in Figure 6 were not satisfactory, because the sample mean and (co)variance estimates in the M-step, obtained from minimizing the KL divergence from the posterior distribution f to the parametric model q η η η , still suffer from sensitivity to the outliers in the observed samples, even after including the Poisson point process component. Similar robustness concerns are widespread in the clustering literature; see, e.g., [START_REF] Davé | Robust clustering methods: a unified view[END_REF] and references therein.

As a workaround, we propose to use robust estimates [START_REF] Huber | Robust statistics[END_REF] of the and (co)variances of Gaussian distributions instead of the empirical means and (co)variances in the M-step 9 . For example, in the case of univariate Gaussian distributions, we use the median and 1.349 times the interquartile range as robust estimators of the mean and standard deviation, respectively. See [1, Section 2.5] for more discussion of this robustness issue, including an alternative solution using the (robust) α-divergence of [START_REF] Basu | Robust and efficient estimation by minimising a density power divergence[END_REF] instead of the KL divergence.

Remark 2. Note that the robustification introduced in this section does not render the PPP component useless. Indeed, removing it would result (experiments not shown) in the introduction of one or several large-variance small-probability Gaussian components in the model, to account for the diffuse part of the posterior discussed in Section II-B, and would also impose to choose L at least as large as the largest value of k (i) in the MCMC samples.

III. ILLUSTRATIVE EXAMPLES

In this section, we will investigate the capability of VAPoRS for summarizing variable-dimensional posterior distributions using two signal decomposition examples; 1) joint Bayesian detection and estimation of sinusoids in white Gaussian noise [START_REF] Andrieu | Joint Bayesian model selection and estimation of noisy sinusoids via reversible jump MCMC[END_REF] and 2) joint Bayesian detection and estimation of astrophysical particles in the Auger project [START_REF]The Pierre Auger Project Design Report[END_REF][START_REF]Properties and performance of the prototype instrument for the Pierre Auger Observatory[END_REF]; see [START_REF] Roodaki | Signal decompositions using trans-dimensional Bayesian methods[END_REF]Chapters 3 and 4] for more results and discussion. We emphasize again that the output of the trans-dimensional Monte Carlo sampler, e.g, the. RJ-MCMC sampler in this paper, is considered as the observed data for VAPoRS.

A. Joint Bayesian detection and estimation of sinusoids in white Gaussian noise

Let us consider the problem of detection and estimation of sinusoidal components introduced in Section I-B where the unknown parameters are the number k of components, the component-specific parameters (a c, j , a s, j , ω j ), 1 ≤ j ≤ k, and the noise variance σ 2 . Since the amplitudes and the noise variance can be analytically integrated out, we focus on summarizing the joint posterior distribution p(k, ω ω ω 1:k | y) of the form illustrated in Figure 1. Therefore, we assume that the proposed parametric model introduced in Section II-B consists of univariate Gaussian components, with means µ l , variances s 2 l , and probabilities of presence π l , 1 ≤ l ≤ L, to be estimated. Moreover, the space of component-specific parameters is Θ = (0, π) ⊂ R.

Before launching VAPoRS, we need first to initialize the parametric model. It is natural to deduce the number L of Gaussian components from the posterior distribution of k.

Here, we set it to the 90 th percentile of p(k | y) to keep all the probable models in the play. To initialize the Gaussian components' parameters, i.e., µ l and s 2 l , 1 ≤ l ≤ L, we used the robust estimates of the means and variances of the marginal posterior distributions of the sorted radial frequencies given k = L. Finally, we set π l = 0.9, for 1 ≤ l ≤ L, and λ = 0.1.

We ran the "robustified" stochastic algorithm introduced in Section II on the specific example shown in Figure 1, for 100 iterations, with L = 3 Gaussian components (note that the posterior probability of {k ≤ 3} is approximately 90.3%). To assess the convergence of VAPoRS, Figure 7 illustrates the evolution of the model parameters η η η together with the criterion J . Two substantial facts showing the convergence of VAPoRS can be deduced from this figure: first, the decreasing behavior of the criterion J M , which is almost constant after the 10 th iteration; second, the convergence of the parameters of the parametric model, particularly the means µ l and proba- bilities of presence π l , 1 ≤ l ≤ L, even though we used a naive initialization procedure. Indeed after the 40 th iteration there is no significant move in the parameter estimates.

As discussed in Section I, one of the main objectives of the algorithm we proposed is to solve the label-switching issue in a trans-dimensional setting. Figures 8 shows the histograms of the labeled samples, i.e., (θ θ θ (i) , z (i) ), with i = 1, . . . , M, along with the pdf's of the estimated Gaussian components (black solid line). Moreover, the summaries provided by VAPoRS for each component are presented in its corresponding panel. We used the average of the last 50 SEM iterations as parameter estimates, as recommended in the SEM literature; see, for example, [START_REF] Celeux | A stochastic approximation type EM algorithm for the mixture problem[END_REF][START_REF] Nielsen | The stochastic EM algorithm: estimation and asymptotic results[END_REF]. Comparing the unimodal distributions of the labeled samples with the ones of the posterior distributions of the sorted radial frequencies given k = 3 shown in Figure 1, which are highly multimodal, reveals the capability of VAPoRS in solving label-switching in a variable-dimensional setting.

Looking at the bottom right panel of Figure 8, the role of the point process component in capturing the outliers in the observed samples, which cannot be described by the Gaussian components, becomes clearer 10 . Note that, without the point process component, these outliers would be allocated to the Gaussian components and would, consequently, induce a significant deterioration of the parameter estimates.

Table I presents the summaries provided using VAPoRS along with the ones obtained using the BMS approach. Contrary to the BMS approach, VAPoRS has enabled us to benefit from the information of all probable models to give summaries about the middle harder to detect component. Turning to the results of VAPoRS, it can be seen that the estimated means are compatible with the true radial frequencies. Furthermore, the estimated probabilities of presence are consistent with uncertainty of them in the variable-dimensional posterior shown in Figure 1.

To observe better the "goodness-of-fit" of the estimated Gaussian components, the bottom panel of Figure 9 depicts 10 The presence of two peaks on the left indicates that additional Gaussian t-components could be added here to get a more accurate approximation of the posterior. See Section III-B for a discussion of the interpretation of the estimated intensity as a residual. their normalized densities 11 , under the posterior distributions of the sorted radial frequencies given k. This figure can be used to validate the coherency of the estimated summaries with the information in the variable-dimensional posterior distribution. It can be seen from the figures that the shape of the pdf's of the estimated Gaussian components are coherent in both the location and dispersion with the ones of the posterior of the sorted radial frequencies.

It is also useful, for the validation of the estimated parametric model, to compare the intensity [see, e.g., 35] intensity (which averages over MCMC samples from all model sizes, and is therefore an example of using the BMA approach; see [1, Chapter 2] for more information). Figure 10 shows such a figure for the specific example of this section where the solid black line indicates the intensity of the estimated parametric model. These figures also indicate the "goodness-of-fit" of the fitted approximate posterior and the true one.

h = L ∑ l=1 πl • N (• • • | μ µ µ l , Σ Σ Σ l ), (10)
Finally, to validate both the estimated probabilities of presence of the Gaussian components and the mean parameter λ of the Poisson point process component, Figure 11 illustrates the posterior distribution of the number k of components together with its approximated versions using VAPoRS. It can be seen from the figure that, in this example, VAPoRS successfully captured the information provided in the true posterior of the number k of components.

Similar results have been obtained for a wide range of configurations (with varying number of sinusoids, radial frequencies, signal-to-noise ratios, . . . ; not shown in the paper). In particular, an example with k = 30 sinusoids is discussed in [START_REF] Roodaki | Signal decompositions using trans-dimensional Bayesian methods[END_REF]Section 3.3.4]; the computation time of our Matlab/C implementation, on this example where L = 30 and M = 20 000, is approximately 15 seconds per SEM iteration (on a laptop with an Intel Core i5 M540 running at 2.53 GHz and 4 GB of RAM).

B. Joint Bayesian detection and estimation of astrophysical particles in the Auger project

As the second illustrative example, we show results on a signal decomposition problem encountered in the international astrophysics collaboration called Auger [START_REF]The Pierre Auger Project Design Report[END_REF][START_REF]Properties and performance of the prototype instrument for the Pierre Auger Observatory[END_REF]. The Auger project is aimed at studying ultra-high energy cosmic rays, with energies in order of 10 19 eV, the most energetic particles found so far in the universe. The long-term objective of this project is to study the nature of those ultra-high energy particles and determine their origin in the universe. These particles are not observed directly. When they collide the earth's atmosphere, a host of secondary particles are generated, some of which, mostly muons, finally reach the ground. To detect these muons, the Pierre Auger Cosmic Ray Observatory was built, which consists of two independent detectors: an array of Surface Detectors (SD) and a number of Fluorescence Detectors (FD). There are in total 1600 SD tanks, each separated from its neighbors by 1.5 kilometers, covering a surface of about 3000 km 2 .

The number of muons and their arrival times can be used as indications of both the chemical composition and the origin of the primary particles (see [START_REF]The Pierre Auger Project Design Report[END_REF][START_REF]Properties and performance of the prototype instrument for the Pierre Auger Observatory[END_REF] for more information). Here, we concentrate on the signal decomposition problem, where the goal is to count the number of muons and estimate their individual parameters from the signals observed by SD detectors; while noting that, to investigate the characteristics of the primary particles, one needs to use information obtained from a few tens of SD's.

This problem has been addressed by [START_REF] Bardenet | An adaptive Metropolis algorithm with online relabeling[END_REF][START_REF] Kégl | Bayesian estimation, the Metropolis-Hastings algorithm, and a simple example[END_REF][START_REF] Bardenet | Single muon response: The signal model[END_REF] in a Bayesian framework, in which they developed an RJ-MCMC sampler to jointly count the muons and estimate their parameters. They run thousands of iterations of the RJ-MCMC sampler on the signals captured by each individual SD tank and, then, aggregate all the samples in a secondary analysis step to make inference on the characteristics of the primary particle [START_REF] Bardenet | Towards adaptive learning and inference : Applications to hyperparameter tuning and astroparticle physics[END_REF]Chapter 7]. To make the whole process more efficient and easier to interpret, it would be preferable to "digest" the MCMC samples associated to each SD tank, and to convey only the resulting summaries to the next inference level. This would also make it much cheaper to store the data required for the second stage of inference, e.g., for traceability or future studies. Therefore, an algorithm that faithfully summarizes the RJ-MCMC samples is needed. In this section, we first briefly describe the problem and then use VAPoRS to relabel and summarize variable-dimensional output samples of the RJ-MCMC sampler developed by [START_REF] Bardenet | An adaptive Metropolis algorithm with online relabeling[END_REF][START_REF] Kégl | Bayesian estimation, the Metropolis-Hastings algorithm, and a simple example[END_REF][START_REF] Bardenet | Single muon response: The signal model[END_REF].

When a muon crosses a SD tank, it generates photoelectrons (PE's) along its track that are, then, captured by detectors and create a discrete observed signal. We denote the vector of observed signal by n = (n 1 , . . . , n N ) ∈ N N , where the element n i indicates the number of PE's deposited by the muons in the time interval

[t i-1 , t i ) [t 0 + (i -1)t ∆ , t 0 + it ∆ ),
where t 0 is the absolute starting time of the signal and t ∆ = 25 ns is the signal resolution (length of one bin). Each muon has two component-specific parameters, namely, the arrival time t µ and the signal amplitude a µ . The absorption process of the photons generated by a muon is modeled by a non-homogeneous Poisson point process with intensity [43, Section 2.2] where p τ,t d (t) is the time response distribution, t d is the risetime and τ is the exponential decay (both measured in ns); see Figure 12 (bottom) for such exponential shape intensities. Then, the expected number of PE's in the bin i is obtained by integrating the intensity [START_REF] Green | Trans-dimensional Markov chain Monte Carlo[END_REF] in the corresponding bin:

h(t | a µ ,t µ ) = a µ p τ,t d (t -t µ ), (11)
ni (a µ , t µ ) = a µ t i t i-1 p τ,t d (t -t µ )dt. (12)
Conditioning on the number k of muons and the vector of parameters t t t µ = (t µ,1 , . . . ,t µ,k ) and a a a µ = (a µ,1 , . . . , a µ,k ), and assuming that the number of PE's in each bin are independent, the likelihood is written as

p(n | k,t t t µ , a a a µ ) = N ∏ i=1 p(n i | ni (k, a a a µ , t t t µ )), (13)
where

We will now illustrate the performance of VAPoRS on a simulated PE counting signal (see [START_REF] Roodaki | Signal decompositions using trans-dimensional Bayesian methods[END_REF]Chapter 4] for results on two other simulated experiments). The observed signal of the illustrative example considered here consists of five muons located at t t t µ = (105, 169, 267, 268, 498) (see Figure 12). The posterior distributions of the number k of muons and sorted arrival times are shown in Figure 13. Note that, in this example, there are two muons with almost equal arrival times, i.e., the third and fourth muons.

Using the BMS approach, the model with four muons would be selected (p(k = 4 | n) = 0.4), although M 5 has an almost identical posterior probability of 0.38. Moreover, observe that the marginal posterior of the arrival time of the third component is bimodal under both M 4 and, more significantly so, M 5 . We ran VAPoRS with L = 6 Gaussian components on the RJ-MCMC output samples shown in Figure 13 (note that p(k ≤ 6 | n) = 0.94). Figure 14 shows the histogram of the labeled samples and the estimated parameters of the components. From the figure, it can be seen that the bimodality effects caused by labelswitching exhibited in Figure 13 is removed completely and the estimated Gaussian components are unimodal and enjoy reasonable variances. In the presented summary, there are four muons with high probabilities of presence corresponding to the ones shown in the bottom row of Figure 13. There are also two other muons with comparatively low probabilities of presence.

In fact, the samples allocated to the point process component shown the bottom row of Figure 14 can be regarded as the residuals of the fitted model, that is, the observed samples which the L Gaussian components in q η η η have not been able to describe. These residuals can be used, as usual in statistics, as a tool for goodness-of-fit diagnostics and model choice.

Figure 15 illustrates the histograms of the residuals of the fitted model for different values of L ∈ {3, 4, 6, 8}. It can be seen from the top left panel of Figure 15 that the distribution of the residuals corresponding to the case where L = 3 contains a few "significant" peaks. The peaks are gradually removed by adding Gaussian components. When L = 4, a component is added at t µ = 261 that captures samples distributed around the most significant peak of the top left panel of Figure 15. However, there still exist a few peaks, particularly around t µ = 173 which are captured when L ≥ 6. However, the distribution of residuals for the case of L = 6 and L = 8 do not differ significantly. Note the decrease of value of λ by increasing L.

Figure 16 compares the normalized intensities of the estimated Gaussian components for 6 ≤ L ≤ 9. When moving from L = 6 to L = 9, the six Gaussian components that are estimated in the case with L = 6 always exist, but additional Gaussian components with very low probabilities of presence are added to the summary, which improve the fit but does not change much the final inference.

IV. MONTE CARLO EXPERIMENT

The examples of Section III have illustrated the capability of VAPoRS to relabel and summarize variable-dimensional posterior distributions encountered in two signal decomposition problems. In order to confirm these findings, we will now investigate more systematically, by means of a Monte Carlo simulation experiment, how faithfully the approximate posterior distribution preserves certain features of the true posterior distribution.

One hundred realizations of the sinusoid detection experiment described in Section I-B (see Figure 1) were simulated and analyzed using the same RJ-MCMC sampler as before. The number of RJ-MCMC iterations was set to 100 000 and the first 20 000 samples were discarded as the burn-in period. Then, the samples were thinned to one every fifth. To initialize the parametric model q η η η in a systematic fashion, we set L to the largest k such that its posterior probability is not less than 0.05. Then, during the process of the SEM-type algorithms, if sufficient number of samples, say, 10, is not allocated to a Gaussian component (or, equivalently, its probability of presence fades to zero), we will remove it from the parametric model and decrease L by one. Using this approach generally results in approximate posterior distributions which are "richer" than those provided by the BMS approach 12 , in the sense that L ≥ k MAP , where k MAP = argmax k p(k|y). To initialize the Gaussian components' parameters, i.e., the means µ l and variances s 2 l , we used as previously robust estimates of the mean and variances of the posterior distributions of sorted radial frequencies given k = L.

Figure 17 compares various features of the fitted approximate posterior distribution q η η , obtained using 100 iterations of VAPoRS, with the corresponding features of the true variable-dimensional posterior distribution. These features are described in the rest of this section.

The scatter plots shown in panels (a), (b), and (c) compare the posterior distribution of the number k of components, i.e., p(k|y), with its approximated version, denoted here by p(k|y), in 100 runs. We only show the posterior probabilities of k = 2 and k = 3 in this comparison, as the other probabilities were close to zero. The digits situated on the right of the points in the panel (a) indicate the number of occurrence of the corresponding event in 100 runs and kMAP = argmax k p(k|y). It can be seen from these three panels that the information in p(k|y) was well preserved by the approximated posterior distributions.

Next we compare the performance of VAPoRS with the one of the "direct" BMA approach (i.e., using the RJ-MCMC samples directly, and not the VAPoRS posterior) in reconstructing the noiseless signal y 0 = D a 1:k . To this end, the estimated reconstructed noiseless signal is defined as

ŷ0 = E(y 0 | y) = ∑ k∈N Θ k E(y 0 | k, θ θ θ 1:k , y) p(k, θ θ θ 1:k | y) dθ θ θ 1:k . (15)
In the direct BMA approach, using the samples generated with the RJ-MCMC sampler, the above integral is approximated by

ŷBMA 0 = 1 M M ∑ i=1 D (i) â a a (i) 1:k (i) ,
where D (i) is the design matrix of the i th vector of the sampled radial frequencies ω ω ω (i)

1:k (i) and â a a (i) 1:k (i)
is the posterior mean of the amplitudes given ω ω ω (i) 1:k (i) and its hyperparameters. To reconstruct the noiseless signal from the approximate posterior q η η η using VAPoRS, one can generate R pairs of samples (k (r) , ω ω ω Panel (d) compares the normalized reconstruction errors when using VAPoRS with the ones of the direct BMA approach in dB, defined as 10 log 10 ŷ0y 0 2

y 0 2 , ( 16
)
where • is the L 2 -norm and we set ŷ0 = ŷBMA 0 and ŷ0 = ŷVAPoRS 0 , when using the BMA approach and VAPoRS, respectively. It can be seen from the figure that the normalized errors of the reconstructed noiseless signals using the compact summary obtained by VAPoRS are quite comparable with the ones obtained using the BMA approach.

Finally, the scatter plots in the last two panels compare the expected number of components in the intervals (0, π/4) and (π/4, π/2) using VAPoRS with, again, the ones obtained using the direct BMA approach. For the BMA approach, the expected number of components in an interval T ⊂ (0; π) is

E(N(T ) | y) = ∑ k∈N E(N(T ) | k, y) p(k | y) ≈ 1 M M ∑ i=1 N (i) (T ) ,
where N (i) (T ) is the number of radial frequencies observed in T on the i th sample. On the other hand, from the summary provided by VAPoRS, the number of components in interval T is

E η η η (N(T ) | y) = L ∑ l=1 πl N (T ; η η η l ) + λ |T | |Θ| ,
where N (T ; η η η l ) denotes the probability of T under the Gaus- sian distribution with parameters η η η l . The figures confirm that the expected number of components in the chosen intervals computed using both approaches are very similar.

The results shown in this section confirmed that the approximate posterior distribution q η η η obtained using VAPoRS preserves faithfully several important features of the true posterior distribution; see [START_REF] Roodaki | Signal decompositions using trans-dimensional Bayesian methods[END_REF]Section 3.4] for more results in this vein, including a numerical investigation comparison of the properties of estimators derived from VAPoRS.

V. CONCLUSION

In this paper, we have proposed a novel algorithm to relabel and summarize variable dimensional posterior distributions encountered in signal decomposition problems when the number of component is unknown. For this purpose, a variable-dimensional parametric model has been designed to approximate the posterior of interest. The parameters of the approximate model have been estimated by means of an SEM-type algorithm, using samples from the true posterior distribution f generated by a trans-dimensional (e.g., RJ-MCMC) Monte Carlo sampler. Modifications of our initial SEM-type algorithm have been proposed, in order to cope with the lack of robustness of maximum likelihood-type estimates.

The relevance of the proposed algorithm, both for summarizing and for relabeling variable-dimensional posterior distributions, has been illustrated on two signal decomposition examples, namely, the problem of detection and estimation of sinusoids in Gaussian white noise and a particle counting problem motivated by the astrophysics project Auger. Most notably, VAPoRS has been shown to be the first approach in the literature capable of solving the label-switching issue in trans-dimensional problems. We have shown that the proposed parametric model provides a good approximation for the posteriors encountered in both applications. Moreover, VAPoRS can provide the user with more insight concerning not only the component-specific parameters but also the uncertainties about their presence.

We believe that this algorithm can be useful in the vast domain of signal decomposition and mixture model analysis to enhance inference in trans-dimensional problems. In particular, it is useful in large-scale applications, such as the Auger project, where storing all samples can be problematic. Theoretical investigations are required in order to extend available existing convergence results to the SEM-type algorithm used in this paper (with correlated input data and Metropolis-Hastings updates). Future work will focus on using VAPoRS to design more efficient adaptive trans-dimensional MCMC methods, as a continuation of the ideas presented in [START_REF] Bai | Divide and conquer: a mixturebased approach to regional adaptation for MCMC[END_REF][START_REF] Bardenet | An adaptive Metropolis algorithm with online relabeling[END_REF].

Figure 1 :

 1 Figure 1: Posterior distributions of k (left) and sorted radial frequencies, ω ω ω 1:k , given k (right) from 100 000 output RJ-MCMC samples. The true number of components is three. The vertical dashed lines in the right figure locate the true radial frequencies. Not shown on the figure: p(k ≤ 1 | y) ≈ 0% and p(k ≥ 5 | y) ≈ 1.9%.

Figure 2 :

 2 Figure 2: Steps to generate samples from the proposed parametric model.

Figure 3 :Figure 4 :

 34 Figure3: Generated samples from an example of the proposed variable-dimensional parametric model. There are L = 3 Gaussian t-components in the model with the means µ µ µ = (0.63, 0.68, 0.73) t , the variances s s s 2 = (0.01, 0.02, 0.01) and the probabilities of presence π π π = (0.8, 0.3, 0.8). The × × × signs indicate the location of the generated random samples. (a) ξ ξ ξ = (1, 0, 1) and θ θ θ = (2, (0.63, 0.72)), (b) ξ ξ ξ = (1, 1, 0) and θ θ θ = (2, (0.63, 0.66)), (c) ξ ξ ξ = (0, 0, 1) and θ θ θ = (1, (0.71)), (d) ξ ξ ξ = (1, 0, 0) and θ θ θ = (1, (0.62)), (e) ξ ξ ξ = (1, 1, 1) and θ θ θ = (3, (0.62, 0.70, 0.73)), (f) ξ ξ ξ = (0, 0, 0) and θ θ θ = (0, ()).

Figure 5 :

 5 Figure5: Proposed variable-dimensional parametric model in a generative viewpoint. Square nodes show fixed variables while circle nodes denote random variables. Note that η η η = (η η η 1 , . . . , η η η L , λ ) and η η η l = (µ µ µ l , Σ Σ Σ l , π l ).

Figure 6 :

 6 Figure 6: Pseudo-code for the proposed SEM-type algorithm

7 :

 7 Evolution of the model parameters along with the criterion J M defined in (8) using 100 iterations of VAPoRS with L = 3 on the RJ-MCMC output samples shown in Figure1.

µFigure 8 :Figure 9 :

 89 Figure 8: Normalized histogram of the labeled samples, that is, the samples allocated to the Gaussian and Poisson point process components, versus the pdf's of estimated Gaussian components in the model (black solid line) using VAPoRS on the sinusoid detection example. The estimated parameters of each component are presented in the corresponding panel.

Figure 10 :Figure 11 :

 1011 Figure 10: Comparison of the intensity of the fitted parametric model obtained with VAPoRS (black line) with an histogram estimator of the intensity (light gray).

Figure 12 :

 12 Figure 12: (top) Simulated observed signal n. (bottom) Intensity of the model h(t | a a a µ , t t t µ ) defined in[START_REF] Green | Trans-dimensional Markov chain Monte Carlo[END_REF]. There are k = 5 muons in the signal with the true arrival times, i.e., t t t µ = (105, 169, 267, 268, 498), indicated by vertical dashed lines.

Figure 13 :

 13 Figure 13: Posterior distributions of the number k of muons (left) and the sorted arrival times, t t t µ , given k (right) constructed using 60 000 RJ-MCMC output samples after discarding the burn-in period. The true number of components is five. The vertical dashed lines in the right figure locate the arrival times.

µ 14 :Figure 15 :

 1415 Figure 15: Normalized histograms of the residuals of the fitted model using VAPoRS with different values of L = {3, 4, 6, 8}.

Figure 16 :

 16 Figure 16: Normalized pdf's of the fitted Gaussian components using VAPoRS with different values of 6 ≤ ≤ 9.

  E(N(π/4, π/2))

Figure 17 :

 17 Figure 17: Comparison of (some features of) the true posterior distribution with its VAPoRS approximation.

Table I :

 I Summaries of the variable-dimensional posterior distribution shown in Figure1; VAPoRS vs. the BMS approach.

	Comp. 1	µ 0.62 0.017 s	π 1	µ BMS 0.62	s BMS 0.016
	2	0.68 0.021 0.22	-	-
	3	0.73 0.011 0.97	0.73	0.012

  p(n i | ni (k, a a a µ , t t t µ )) is a Poisson distribution with the mean ni (k, a a a µ , t t t µ ). Then, assuming independence of the muons, the expected number of PE's in the i th bin, i.e., ni (k, a a a µ , t t t µ ), given k, t t t µ , and a a a µ becomes ni (k, a a a µ , t t t µ ) =

k ∑ j=1 ni (a µ, j , t µ, j ).

See, however, the intensity plot provided in Section III (Figure10) as an example of a BMA summary related to a component-specific parameter.

Reporting additional models would improve the situation for BMS, in this example, but would not directly provide a "probability of presence" for each component, as our approach does.

Any parametric family of d-variate distributions could be used at this point. As often in the literature[START_REF] Bai | Divide and conquer: a mixturebased approach to regional adaptation for MCMC[END_REF][START_REF] Stephens | Dealing with label switching in mixture models[END_REF][START_REF] Sperrin | Probabilistic relabelling strategies for the label switching problem in bayesian mixture models[END_REF][START_REF] Yao | Model based labeling for mixture models[END_REF], the Gaussian distribution is chosen as a convenient mean of describing a "compact" and unimodal d-dimensional distribution, nothing more. Note that, because the Gaussian distribution is supported by R d , our parametric model is actually defined on a variable-dimensional space bigger than Θ Θ Θ if Θ is a strict subset of R d .

which are multimodal, because sorting the frequencies at fixed k does not properly resolve the label switching problem

Homogeneity of the PPP component has been assumed for the sake of simplicity, but more elaborate (non-homogeneous) models are easily accommodated by our approach, if needed. A non homogeneous PPP is needed, in particular, if |Θ| = +∞; indeed, an homogeneous PPP on a such a set would generate an almost-surely infinite number of components. As an example: if one of the component-specific parameters is positive-an amplitude parameter in a positive mixture, say-then a lognormal distribution, a gamma distribution or uniform distribution on a bounded segment can be used to build a (bounded) non-homogeneous intensity.

As an alternative to the simple plug-in method proposed here, robustification could be achieved by using heavy-tailed distributions (e.g., Student distributions) instead of Gaussian distributions in the parametric model. Note that this approach, more elegant perhaps, makes the M-step more complicated since (at least a few steps of) an optimization algorithm must be used to find the new value of the parameters.

To obtain the normalized densities, first, we normalized the estimated pdf's to have their maximum equal to one. Then, we multiplied the estimated probability of presence of each Gaussian component to its corresponding normalized estimated pdf. Thus, the maximum of each normalized density is equal to the corresponding estimated probability of presence.

Later, in a post-processing step, since each Gaussian component has been endowed with a probability of presence π l , with 1 ≤ l ≤ L, one can decide to discard the ones with π l smaller than a certain threshold; see[START_REF] Roodaki | Signal decompositions using trans-dimensional Bayesian methods[END_REF] Section 3.4.3] for more discussion about this idea.

ACKNOWLEDGMENT

The authors would like to express their gratitude to B. Kégl for his collaboration and providing data for the Auger example.