{'fraction_non_alphanumeric': 0.03026299987990873, 'fraction_numerical': 0.009150954725591449, 'mean_word_length': 4.491427064099183, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian modeling of source confusion in LISA data

One of the greatest data analysis challenges for the Laser Interferometer Space Antenna (LISA) is the need to account for a large number of gravitational wave signals from compact binary systems expected to be present in the data. We introduce the basis of a Bayesian method that we believe can address this challenge, and demonstrate its effectiveness on a simplified problem involving one hundred synthetic sinusoidal signals in noise. We use a reversible jump Markov chain Monte Carlo technique to infer simultaneously the number of signals present, the parameters of each identified signal, and the noise level. Our approach therefore tackles the detection and parameter estimation problems simultaneously, without the need to evaluate formal model selection criteria, such as the Akaike Information Criterion or explicit Bayes factors. The method does not require a stopping criterion to determine the number of signals, and produces results which compare very favorably with classical spectral techniques.


I. INTRODUCTION
The Laser Interferometer Space Antenna (LISA) is designed to detect gravitational radiation from astrophysical sources in the 10 −2 mHz to 100 mHz band [1]. The sensitivity of LISA is such that very many such sources should be detectable. Paradoxically, this sensitivity will also make signal identification somewhat problematic: parts of the band will likely be swamped with tens of thousands of signals. One of the most abundant classes of source will be close-by white dwarf binaries, producing signals from 0.1 mHz to 3 mHz. There will be source confusion below 1 mHz and resolvable sources above 5 mHz; the 1 mHz to 5 mHz band for LISA therefore presents a tremendous data analysis challenge, potentially containing up to 10 5 sources [2,3,4]. LISA's ability to detect and characterize other astrophysical sources will be greatly helped if the thousands of background signals from binary systems can be identified. For a detailed look at the population of binary systems that produce signals in LISA's operating band, and how they affect LISA's performance, we direct the reader to Barack and Cutler [2] and to Nelemans et al. [5].
We approach this problem from a new direction. Markov chain Monte Carlo (MCMC) techniques have been demonstrated to be especially suited to parameter estimation problems involving numerous parameters [6]. We have previously used the Metropolis-Hastings (MH) algorithm [7,8] in other gravitational radiation problems, such as estimating astrophysical parameters for gravitational wave signals from coalescing compact binary systems [9] or pulsars [10,11]. It is our belief that MCMC methods could provide an effective means for identifying and characterizing the thousands of background binary signals to be found within the LISA data.
The method that we present in this paper is not a source subtraction approach [12], but one that identifies and characterizes binary produced periodic signals in the data. Signals that are sufficiently large in amplitude will have their parameters estimated. Sources that are weak will contribute to the noise; our method also produces an estimate of the overall level of the noise. We show that the noise level estimate from our method depends on the inherent detector noise level, and also the presence of unidentified signals. MCMC methods are robust and dynamic, and we believe that ultimately it will be possible to use them with LISA data to estimate simultaneously the parameters associated with a wide range of source types occurring in the presence of many thousands of white dwarf binaries.
In this paper we present the results of a simulation study, comprising a data stream of m sinusoidal signals embedded in gaussian noise. Although our simulation study is simple, it does highlight a number of issues relevant to the real LISA data analysis problem. For example, the number of signals present in the data, m, will be unknown.
Our Bayesian approach does not need to fit each model with m signals, for m = 1, . . . M, and then select the best fitting model via the evaluation of Bayes factors. The evaluation of Bayes factors [13,14] requires computation of the marginal likelihoods and thus marginalizations over the parameter vectors of each model. This is a formidable computational problem when the dimension of the parameter space is large. A shortcut to the calculation of Bayes factors, the harmonic mean of the likelihood values [15], is known to be unstable because the inverse likelihood does not possess a finite variance. Other large sample approximations to the Bayes factors such as the Bayesian Information Criterion (BIC), also referred to as Schwarz Criterion, and the related penalized likelihood ratio model choice criterion, Akaike Information Criterion (AIC), have been shown to be inconsistent when the dimension of the parameter space goes to infinity [18]. The newly developed deviance information criterion (DIC) [19] is known to be controversial in mixture models. While MCMC algorithms like the Gibbs sampler and MH algorithm yield posterior distributions of the parameters, they do not provide marginal likelihoods. An indirect method of estimating marginal likelihoods from Gibbs sampling output has been developed [16] and extended to output from the MH algorithm [17]. These, however, are impractical when the number of candidate models is very large, as is the case for LISA data. Therefore, we use another strategy, the reversible jump algorithm [27], that samples over the model and parameter space in order to estimate posterior model probabilities/marginal likelihoods. We consider the number of sinusoids as an additional parameter and determine its marginal posterior distribution. Its modal value will give us the model with the most probable number of sinusoids. As illustrated in Section II, a Bayesian analysis naturally encompasses Occam's Razor [20] and a preference for a simpler (smaller m) model. We thus address and solve the detection and estimation of sinusoids simultaneously. In addition, our MCMC method is better than a classical periodogram at resolving signals that are very close in frequency, and we provide a detailed discussion of how to identify these signals. Finally, our method infers the noise level in the data together with the parameters of each of the m sinusoids.
The problem of identifying an unknown number of sinusoids is neither new nor simple [21,22]. Whereas previous studies have looked for a handful of unknown signals, here we show results for 100 signals. Another benefit of using MCMC methods is that computation time scales roughly linearly as the number of parameters increases, and does not show an exponential increase in time [6]. In the future we will make the character of the signals more realistic, taking into account the orbit of the LISA satellites and the nature of the inspiral of binaries. With our study of sinusoids we hope to inform LISA data analysis researchers of another possible avenue for characterizing a large number of background signals.
The rest of the paper is organized as follows: Section II illustrates that Bayesian inference automatically implements the principle of Occam's Razor. Our Bayesian method and posterior computational algorithms are described in Section III. In Section IV we present the results of this study using synthesized data. We believe that this method offers great hope for signal extraction in real LISA data, and this point is discussed in Section V.

II. OCCAM'S RAZOR
The notion that Bayesian inference naturally implements a quantitative version of Occam's Razor [20] is well known, but this property is sufficiently central to our approach to LISA analysis to warrant a very simple demonstration: We imagine two physical models, M 1 and M 2 constrained by a single datum, d. M 1 has one parameter, a 1 , to describe the datum. M 2 uses the sum of two parameters, s = a 1 + a 2 to describe the same datum. Clearly, all other things being equal, if the datum is equally consistent with both models, we would prefer M 1 over M 2 on the grounds that M 1 is 'simpler'. To see how this naturally occurs within the Bayesian framework we proceed by considering the ratio of probabilities (the odds ratio) of the two models: If we set the prior odds ratio to unity, reflecting no prior preference for either model, the right-hand side of the above equation is simply the ratio of marginal likelihoods (sometimes called the evidences) of the two models. This ratio is conventionally known as the Bayes factor for the comparison. We will take the priors for our parameters a 1 and a 2 to be each uniform in the range 0 → R. Under M 1 we therefore have p(a 1 |M 1 ) = 1/R. The likelihood under M 1 is assumed to be p(d|a 1 , M 1 ) = δ(d − a 1 ), where δ denotes the Dirac function, i.e. δ(x) = 1 if x = 0 and 0 otherwise. This yields Under M 2 with parameter s = a 1 + a 2 , the prior probability density function (pdf) for s is the convolution of the two uniform prior pdfs for a 1 and a 2 and therefore has the 'triangle' form so the evidence for M 2 is (see Fig. 1). If the datum lies in the range 0 < d < R, the Bayes factor is M 2 has more flexibility than is necessary to explain the datum and so penalizes itself by spreading its evidence more thinly. We can also see that the effect is not simply dependent on the number of parameters. If M 2 depended on just a single parameter but with a uniform prior on this parameter over 0 → 2R, then it too would be discriminated against (in the case 0 < d < R) for using more parameter space to do the same job as M 1 .
This quantitative ability to identify the least flexible model consistent with the data makes Bayesian inference the natural methodology in problems where the number of parameters is unknown. Indeed we do not need to consider explicitly Bayes factors at all if we consider the number of parameters itself to be a parameter, and determine its marginal posterior probability in the usual way. This is the approach we will apply to the the LISA white dwarf binary problem, as described in the next section.
where the noise terms ǫ j are assumed to be i.i.d. N(0, σ 2 m ) random variables. The signal of model M m has the form so that each sinusoid component is characterised by one frequency and two amplitudes.
Model M m is therefore characterized by a vector of 3m + 1 unknown parameters which we The objective is to find the model M m that best fits the data and to estimate its parameters. We use a Bayesian approach as in [25], but instead of fitting each model separately and selecting the best using Bayes factors, we treat the number m of unknown sinusoids as an additional unknown parameter and determine the mode of its marginal posterior distribution (together with the posterior distribution of all the other parameters). The joint probability that these data d arise from the parameter vector a m and model M m is given by For simplicity we take the prior distribution of the model dimension parameter m as uniform over {0, . . . , M}. The variances σ 2 m are given noninformative inverse gamma priors, discussed in Section III C.
Again, for simplicity our calculations use dimensionless frequencies with a Nyquist frequency of 0.5, and we use uniform priors for the component frequencies f Given m, σ 2 m and the frequency vector f (m) , our model (5) is a linear regression model which can be written in matrix form as   Note that these priors are not uninformative, and we express the parameter space in Cartesian coordinates simply because this is convenient for the implementation of the MCMC algorithms. A model expressed in polar coordinates with uniform priors on amplitude and phase would correspond to a different prior distribution for the transformed amplitudes in Cartesian coordinates. The effect of these different priors on the frequency estimates is generally negligible when the data is informative, and is discussed in [25]. The posterior pdf of the frequency reaches its maximum at the same value for both prior choices but with different curvatures at this maximum.
By applying Bayes' theorem, we obtain the posterior pdf for our model parameters of where p(d) = M i=0 p(m, a m )p(d|m, a m ) da m . The direct evaluation of the normalization constant p(d) is difficult due to the (3m + 1)-dimensional integration involved. Moreover, the computation of marginal posterior pdfs would require subsequent 3m-dimensional integration. To overcome this problem, we use sampling-based MCMC techniques to carry out posterior inference (see [6] for an introduction and overview.) These only require the unnormalized posterior p(m, a m |d) ∝ p(m, a m )p(d|m, a m ) to sample from Eq. (9) and to estimate the quantities of interest. However, in the present context the overall model does not have fixed dimension and classical MH techniques [7,8] cannot be used to propose transdimensional moves. We therefore use the Reversible Jump Markov Chain Monte Carlo (RJMCMC) algorithm for our model determination [27,28]. Additionally we use the delayed rejection (DR) method [29,30] for transitions within the same model. This allows better adaptation of the proposals in different parts of the state space by allowing the choice of the proposal distribution to depend on the proposed but rejected state as well as the current state.

A. The RJMCMC for model determination
To sample from the joint posterior p(m, a m |d) we construct a Markov chain simulation where m is the current number of signals. When a new model is proposed we attempt a step between state spaces of different dimensionality.
Suppose that at the nth iteration of the Markov chain we are in state (k, a k ). If model M k ′ with parameter vector a ′ k ′ is proposed, a reversible move has to be considered in order to preserve the detailed balance equations of the Markov chain. Therefore the dimensions of the models have to be matched by involving a random vector r sampled from a proposal distribution with pdf q(r), say, for proposing the new parameters a ′ k ′ = t(a k , r) where t is a suitable deterministic transformation function of the current state and r. Here we focus on transitions that either decrease or increase models by one signal, i.e. k ′ ∈ {k − 1, k + 1}.
We use equal probabilities p k →k ′ = p k ′ →k to either move up or down in dimensionality, and without loss of generality we consider the upward move k ′ = k + 1.
If the transformation t k →k ′ from (a k , r) to a ′ k ′ and its inverse t −1 k →k ′ = t k ′ →k are both differentiable, then reversibility is guaranteed if we define the acceptance probability for increasing a model by one signal as where is the Jacobian determinant of this transformation [27]. In this context, we suggest two types of transformations, 'split-and-merge' and 'birth-and-death'.

Split-and-merge transitions
For a 'split' transition we randomly choose one of our signals with parameter subvec- This signal is chosen by sampling i uniformly from {1, . . . , k}. The proposed parameter vector a ′ k ′ comprises all the other (k − 1) subvectors of a k and two additional 3-dimensional subvectors, say a ′ i 2 ) each with half the amplitude of a (i) , but same frequency, to replace a (i) . A three-dimensional Gaussian random vector (with mean zero), r = (r A , r B , r f ), changes the current state a (i) to the two resulting states a ′ (1i) , a ′ (2i) through a linear transformation The inverse transformation t −1 k →k ′ := t k ′ →k accounts for the merger of two signals and can be written as Note that the determinant of the Jacobian of the transformation t k →k ′ , (i.e. the determinant of the above 6-by-6 matrix) is |J k →k ′ | = 2, and that of its inverse is 1/2. Thus, the acceptance probability for increasing a model by one signal is given by By analogy, the acceptance probability for the reversible move of a fusion of two signals is given by where , for q(r). Suitable values for the variances of this multivariate normal distribution can be chosen by considering their effect on the acceptance probabilities. First, we consider the effect of a small proposal variance for a splitting transition. In this case the proposal has an insignificant effect on the likelihood since the two new signals in the model function are almost linearly dependent. On the other hand, the resulting large value for q(r) considerably decreases the acceptance probability when proposing a split, and increases it when a fusion of two signals is proposed. Now we consider the effect of a large proposal variance: the value of q(r), and therefore its influence on the acceptance probability, is moderate. However it causes the likelihood to change considerably, resulting in a small acceptance probability. The choice of σ 2 A , σ 2 B and σ 2 f is therefore an important consideration for improving mixing. In each iteration we set σ 2 A , σ 2 B equal to the noise level σ 2 m of the current model m. The posterior precision of the frequency in a single-frequency model depends on the signal to noise ratio γ = (A 2 + B 2 )/σ 2 and the number of samples N of the data set [31]. Using a Gaussian approximation to the posterior pdf of the frequency [31], its standard deviation is given by σ ′′ f = (2πγ) −1 48/N 3 . This yields a distance in frequency for which two sinusoids can still be identified as distinct, neglecting the interference of other sinusoids. We therefore use σ ′′ f as a frequency perturbation when splitting two sinusoids and use a normal distribution with this particular standard deviation when merging.

Birth-and-death transitions
A 'birth' transformational step simply creates a new signal with parameter triple a ′ (i) independent of other existing signals in the current model M k . The one-to-one transformation in this case is very simply given by t k →k ′ (r) = r = a ′ (i) . The inverse ('death') transformation that annihilates signal i ′ , t −1 k →k ′ := t k ′ →k , has form t k ′ →k a ′ (i) = a ′ (i) = r. The Jacobian for both of these is 1. The acceptance probability for the creation process is therefore and that for the annihilation process is As for split and merge transitions, a bold proposal distribution q results in a small acceptance probability due to the strong effect on the likelihood, whereas timid proposals have minor effects on the likelihood but are often rejected due the higher values of the proposal distribution. An effective way to make a proposal for the frequencies is to base it on the Schuster periodogram of the data [32], given by where are the real and imaginary parts from the sums of the discrete Fourier transformation of the data. This technique has already been applied successfully by [22]. We use proposals from a normal distribution for the amplitudes with zero mean and a standard deviation derived from the average amplitude of all remaining amplitudes of the current state of the Markov chain.
Classical MCMC methods could be used for transitions within a particular model M m , however we use an adaptive MCMC technique here. The delayed rejection (DR) method has been introduced by [29,30,33] and several of us have successfully applied it to estimate the frequency and frequency derivative of potential gravitational radiation signals produced by a triaxial neutron star [11].
B. The delayed rejection method for parameter estimation As sampling progresses, suppose that at the nth iteration the state of the Markov chain is a = a m from model M m . We can choose a new state within the same model by first sampling a candidate state a ′ from a proposal distribution q 1 (a ′ |a) and then accepting or rejecting it with an MH-probability α 1 (a ′ |a) depending on the distribution of interest.
When a proposed MH move is rejected, a second candidate a ′′ can be sampled with a different proposal distribution q 2 (a ′′ |a ′ , a) that may depend on the previously rejected proposal. To preserve reversibility of the Markov chain and thus to comply with the detailed balance condition, the acceptance probabilities for both the first and the second stage are given by [30] and We therefore apply distinct types of DR transition for the amplitudes and the frequency of a sinusoid, and these are considered below. The transitions are performed randomly and with equal probability for a randomly chosen sinusoid i.

Proposing new amplitudes
Our amplitude proposal distribution is multivariate Normal, with a covariance matrix that depends on the basis functions of the sinusoids. The closer two sinusoids are in frequency the more correlation there is in their recovered amplitudes.
In [22] the amplitudes were regarded as nuisance parameters and integrated out by treating them as parameters of a multiple regression model with a conditional posterior that is normally distributed with known mean and covariance matrix according to the g-prior. If deemed necessary the amplitudes can be updated by Gibbs sampling. However, the computation of the covariance matrix involves determining the inverse of D (m) (f ) T D (m) (f ) for each iteration, which is impractical considering the number of sinusoids expected for LISA.
Mindful of this we consider the interference between the sinusoids by proposing those pairs of sinusoids that have the strongest correlation, namely sinusoid pairs which are closest in frequency. Therefore for the randomly chosen sinusoid i we also involve the sinusoid i c which is the closest in frequency. The matrix for this subset of basis vectors reduces to which results in a covariance matrix merely of type 4 × 4 for computing proposals for a sinusoid pair (i, i c ). Now suppose that a (i,ic) = A T is a vector containing the amplitudes of the sinusoid pair (but not their frequencies). Then is the covariance matrix of the proposal distribution q 1 (a ′ (i,ic) |a (i,ic) ) = N(a (i,ic) , Σ (i,ic) ), where r is a factor controlling the step size of the proposal.
If the first proposal is rejected, another attempt is made by proposing the amplitudes of a single sinusoid i without considering correlations.
The variance of the proposal is assessed by the mean amplitude of all signals which is given bȳ is the subvector of a containing the amplitudes of sinusoid i.

Proposing a new frequency
A new frequency is proposed as follows: In the first stage a new frequency for sinusoid i is sampled from a proposal density that is proportional to the periodogram q 1 (a ′ (i) ) ∝ (0, 0, C(f )) T and independent of the actual stage. This is similar to the sinusoid proposal scheme for RJMCMC. The main objective of this stage is to coarsely scan the whole parameter space for frequencies.
The event of a rejection suggests to sample from the local frequency mode and a proposal is made conditional on the actual state of the frequency by slightly perturbing the state. In the same manner as in the split-and-merge transition the perturbation is oriented on the possible achievable accuracy σ f (m) frequency by [25]. This yields the proposal q 2 (a ′′ (i) |a (i) ) = N a (i) , diag 0, 0, σ f (m) i and aims to draw representative samples from the local mode in the second stage.

C. Updating the noise parameter
The sum of the squared residuals between the model and the data, taken from the likelihood in Eq. 7, is Using this, we choose a vague prior for the noise parameter σ 2 m , defined by IG(α, β) = IG(N p /2, N p · S 2 p /2) with a shape parameter α = N p /2 = 0.001 and a scale parameter β = N p S 2 p /2 = 0.001. This yields N p = 0.002 and S 2 p = 1 for the parameters of the vague prior. The full conditional distribution is used for drawing samples for σ 2 m in a Gibbs update.

D. Initial values
The initial values of a Markov chain are crucial for the length of the burn-in period needed to converge to the real posterior distribution. We could start with and empty model, M 0 , but it is obvious that it would then take the sampler many steps to find all the signals. Instead, we perform a Fast Fourier Transformation (FFT) of the data and use this to generate our initial values. We use the frequencies corresponding to the local maxima in the periodogram, i , respectively. Theoretically we could use all the local maxima as initial values, but as most of them are due to noise we select only those that exceed a certain threshold. We set this threshold low, as it is easier to delete non-relevant sinusoids than create good ones.
The frequency resolution depends on sample size N, suggesting that the convergence of the Markov chain is also dependent on N. As we will see, spectral estimates based on the FFT are significantly worse than those we obtain by our MCMC method, but they are sufficient to serve as initial values.

E. Identifying the sinusoids
Although the RJMCMC method enables us to select the most probable model, we still encounter the label-switching problem. This is a general problem caused by the invariance of the likelihood under relabelling of the sinusoidal components, and has been extensively discussed in the context of mixture models [34]. During the MCMC simulation, parameter triples are constantly changing their affiliation to individual sinusoids, either due to the creation and annihilation of sinusoids or following transitions within the same model. We therefore need an additional step in our analysis if we are to break this symmetry and talk meaningfully about individual sinusoid components. This step involves associating the samples in the final Markov chain with particular sinusoids, which we know neither by number nor by location.
The sinusoids that are contained in the modelM, with the highest posterior probability of m, are permutations ofm coexistent sinusoids from a list of unknown size. We will assume therefore that this list is the same size as the upper limit m rmmax of the marginal posterior of m. However, in practice the coexistent sinusoids inM are usually clear, and it is rather unlikely that others present in higher order models could replace them. The parameter vector sampled in each iteration of the Markov chain (corresponding to model M) is a permutation ofm parameter triples determiningm out of m max sinusoids. The problem is to determine which parameter triple belongs to which sinusoid. We find that the parameter that contributes significantly to identifying a sinusoid is its frequency. Thus in order to obtain the dominant sinusoids that characterize modelM we calculate the marginal posterior of the frequency and obtain them strongest peaks together with their frequency ranges by finding the threshold that separates those peaks.
Since we have to deal with a vast number of output samples grouped within very small regions we cannot apply kernel density estimates or histograms, as the required fixed bin size for a histogram would be too small to be feasible. Instead, we use a variable bin size and calculate densities for a fixed number of samples per bin. We initially sort all individual parameter triples for all MCMC samples of the considered modelM by frequency.
If we have generated n MCMC samples of modelM during a run then we have to deal withmn parameter triples and hence f 1 < . . . < f nm ordered frequency samples. The density can be assessed by calculating the frequency range spanned by a fixed number of sorted frequencies. Themn parameter triples are fromm sinusoids and therefore we can expect n frequency samples per peak since we assume the number of peaks to be similar to the number of sinusoids. Hence the fixed number of sorted frequencies that spans the frequency ranges must be some fraction, r, of the number of MCMC samples n, and is the counterpart to the required bin width in a histogram or the bandwidth h of a kernel density estimate f(f ) = (2hnm) −1 nm i=1 I |f −f i |<h with uniform kernel, where I is the indicator function. The advantage of this approach is the automatic adaptation of the bandwidth to the situation by involving the information of the expected parameter triples per peak. We found r = 0.05 (5%) to be a good value, and hence rn serves as an estimate of the number of samples needed for assessing the spans of the frequency ranges. In analogy to the kernel density estimate with uniform kernel, where the number of samples are counted that fall into a range of length 2h, the density value ρ j corresponds to each frequency sample j and its rn − 1 subsequent samples that fall into a frequency range of length (f j+rn−1 − f j ).
Hence the density is given by where j = 1, · · · , n(m − r). Since each ρ j comprises the samples j, · · · , j + rn − 1 this has to be considered later when deriving spans for the peaks. We find the smallest density threshold l that separatesm distinct peaks with respect to the values of ρ j . The frequency range for each peak k ∈ {1, · · · ,m} is [f j k,start , f j k,end +rn−1 ] where j k,start and j k,end are the indices of the first and the last member of the set of ρ j 's in peak k, respectively. Due to the fact that we always focus on frequency ranges that contain a fixed amount of frequency samples we efficiently deal with large frequency ranges of low density. This technique is fast and requires a minimum of memory. The greatest computational cost is in sorting the frequencies, although this can be carried out fairly quickly using a heap sort.
It is still possible that individual peaks contain more than one sinusoid or even none. This can be assessed by the histogram of the number of sinusoids in the MCMC samples for the restricted frequency range under consideration. To separate more than one present sinusoid, we then consider the two amplitudes and apply an agglomerative hierarchical cluster analysis that involves all three parameters. We use a modified Ward technique [35] that minimizes the within-cluster variance using a normalized Euclidean distance between the parameters by adjusting the frequency range to the much larger range of the amplitudes. We use the software package R [36] for this task. The Ward technique starts with each parameter triple belonging to a singleton cluster. Iteratively, cluster pairs are joined that produce the smallest possible increase in within-cluster sum of squares. In this particular case we have no difficulty in detecting when to stop the agglomeration, which is an important issue in cluster analysis, as we know the expected number of sinusoids in a peak and hence the number of clusters. However, due to the vast number of parameter triples involved here it is not possible to carry out a cluster analysis simultaneously on all of the data. Instead we have to divide the set of samples into randomly chosen subsets of equal size and perform separate cluster analyses (with R) for each of these subsets. Finally, we perform a cluster analysis on the median points of the subset clusters to allocate each of those clusters to a super-cluster. Each single parameter triple is then allocated to a super-cluster and hence to a presumed sinusoid. Those parameters which fall into small peaks below the threshold can be allocated to an additional noise contribution.

IV. SIMULATIONS
We created an artificial data set of 1 000 samples containing 100 random sinusoids in gaussian noise. The amplitude coefficients were chosen randomly in the range [−1 . . . 1] and the noise standard deviation was σ = 1, making our maximum signal to noise ratio,   Fig. 2(b). Note that σ m decreases with higher model order m since a model comprising more sinusoids accounts for more of the available power. p(f i |d) as described in [31]. An advantage of using a Bayesian spectral analysis is the subsequent ability to calculate confidence areas for the spectrum by grouping our MCMC samples and calculating posterior confidence intervals for frequency bins. An appropriate width for the bins can be estimated using the frequency resolution of the spectrum, σ f = (2πγ) −1 (48/N 3 ) 1/2 , given by [25]. In this example the choice of 20 000 bins is sufficient to resolve sinusoids of γ ≤ 1.4. Fig. 3 compares the known real sinusoids from which the data set had been created, the Bayesian spectral density estimate as described above, and the classical periodogram derived from Eq. 15. We note here that Fig. 3(b) and Fig. 3(c) are spectral densities, whereas Fig. 3(a) is a plot of the component energies, defined as N(A 2 i + B 2 i )/2. The theoretical spectral density of these components would consist of delta functions of formally zero width and hence infinite density. Furthermore the accuracy of the Bayesian spectral estimate is better than the one obtained from the periodogram, which also also leads to different values on the ordinate. The Bayesian spectrum estimate usually reveals much higher values in the density since the same energy contribution from a sinusoid is distributed over a smaller frequency range. However the locations of the peaks can be compared directly, which is the main purpose of Fig. 3.
We have picked out three different frequency ranges for special attention. Each contain two close sinusoid pairs and we will refer to these as Region In Region-1 we see two sinusoids with a frequency gap of approximately 1/N = 0.001 which, for a discrete Fourier transform, is the size of a single frequency bin. Region-2 and Region-3 contain sinusoids with smaller frequency gaps of 0.000 413 and 0.000 517, respectively. The posteriors for Region-1 are strongly peaked at the true values. The same is true for Region-2 and Region-3, but we observe a larger uncertainty in their frequencies and even more in their amplitudes especially when the frequencies of the sinusoids are very close.
The reason for this is clear: when the frequencies of the sinusoids are very close together they slowly beat against each other, and for these particular signals the beat period is greater than the observing time so that their total energy is difficult to assess.
The spectral estimates for these regions are shown in greater detail in Fig. 5. Rather than scale the frequency posterior by the expectation of the energy, as in Fig. 3, here we scale by the 2.5%, 50% and 97.5% quantiles of the energy for each frequency bin to indicate the uncertainty in the value of the spectral density.
The gap between the two sinusoids in Region-1 of Fig. 5 is about 1/N, and the sinusoid with the lower frequency has a significantly larger energy than its neighbor, so its frequency is inferred with greater accuracy. Region-2 and Region-3 however comprise pairs of considerably closer sinusoids, making it considerably more difficult to infer their separate frequencies. show energy densities (inferred energy per unit frequency). The peaks generated by the Bayesian method are higher than those of the periodogram, reflecting its greater ability to constrain the frequency.
The corresponding uncertainties in their amplitudes can be seen in both Fig. 4 and Fig. 5, which shows large energy density values, especially in the area where the frequencies overlap.
Although spectral estimation is an important topic, for LISA analysis a more pertinent issue is the identification of orbital parameters for the compact binary systems generating the signals. In our example, we can therefore consider the joint posterior of the frequency f i and the energy E = N(A 2 i + B 2 i )/2 for the six sinusoids in our three regions (Fig. 6). Unlike in Fig. 5, the ordinate now displays an energy rather than an energy density, and the contour is at the 95% credible level.
To better point out the strength of this Bayesian approach, in Fig. 7 we display 95% confidence areas for a wider frequency range, including more sinusoids. As stated in [25], "if the signal one is analyzing is a simple harmonic frequency plus noise, then the maximum of the periodogram will be the best estimate of the frequency that we can make in the absence of additional prior information about it". This is demonstrated in Fig. 7 showing the concordance in frequency estimates for the well-separated sinusoids from both the joint posterior and periodogram plots. However, the periodogram peaks are significantly wider than the 95% areas, and are clearly sub-optimal for closely spaced sinusoids. In fact the posterior probability density for frequency is the exponent of the ratio of the periodogram C(f ) and the noise variance σ 2 [25]. The Bayesian approach therefore takes account of the noise variance in the estimation process, which explains why the confidence regions are significantly narrower than the periodogram peaks.
This Bayesian method therefore shows great power when tackling strong signals closely spaced in parameter space (in this case, frequency). In addition it delivers confidence intervals for the parameters (frequency and amplitude) and can take account of relevant prior information when applied to LISA data.  It is obvious that stronger sinusoids are easier to detect and separate, so model 2 is favored with much smaller frequency gaps in the left-hand column of Fig. 8 than in the right-hand column. Furthermore, in both the left and right-hand columns the separation ability is worse for small phase shifts and generally improves as the phase shift increases, but this again depends on the noise level. The most striking effect can be seen in the right column at a phase shift of 3π/4. For a small frequency gap the sinusoids can be detected but not separated, then from 14% of a (1/N)-step model 0 starts to be favored. One might expect an increase in the probability of of model 2 at this stage, but in fact their frequency and phase separation mean the pair almost annihilate each other over the time period considered, so that the resultant cannot be seen against the noisy background. With further separation the effect is reduced, and for larger frequency gaps (30% of a (1/N)-step) model 2 does indeed become dominant. This phenomenon can also be seen in the periodogram revealing a single peak slowly decreasing in height from 1% to about 19% of a (1/N)-step. Then the signal has the same height as the surrounding noise and is slowly increasing beyond 20% of a (1/N)-step. Note that for a different observing period, say t ∈ {t start , . . . , t start + N − 1}, the results would have been different. However even in the worst case the detection can be achieved within 40% of a (1/N)-step for the weaker sinusoids and within 26% for the stronger sinusoids.

V. DISCUSSION
In this paper we have presented a Bayesian approach to identifying a large number of unknown periodic signals in a set of noisy data. Our reversible jump Markov chain Monte Carlo method can be used to estimate the number of signals present in the data, their parameters, and the noise level. This method compares favorably with classical spectral techniques. Our approach allows for simultaneous detection and parameter estimation, and does not require a stopping criterion for determining the number of signals.
Although the parameters of even strong components are not well determined when they are sufficiently close together in frequency, we still obtain useful confidence intervals. Importantly, the noise level is itself a parameter in the overall fit so that the energy present in the data is automatically allocated to either signal or noise. The large uncertainties in amplitudes we obtained for the closely-spaced components in our example of 100 sinusoids will, for a practical problem, be lessened by a sensible choice of priors. A prior similar to the g-prior, using a diagonal matrix and a hyper-parameter g as a multiplicative factor, might be a good choice. For the real LISA data we would use priors for the signal amplitudes that would reflect our astrophysical knowledge.
The motivation for our research is to address the difficulty that LISA will ultimately encounter in what is loosely called the confusion problem. LISA may see as many as 100 000 signals from binary systems in the 1 mHz to 5 mHz band. We therefore view our work as a powerful new technique for identifying and characterizing these signals in the LISA data stream.
The work presented here is of course a highly simplified toy problem: we are neglecting signal modulation due to LISA's orbit and beam pattern and are not considering an appropriate data model for LISA. In addition, the signals from compact binary systems differ from simple sinusoids. However, we believe that it demonstrates the applicability of the approach to LISA data analysis, and the next step is to deal with these more complicated signals and to develop a realistic strategy for applying our MCMC methods. We do not claim that this will be a trivial extension; in fact, we acknowledge the complexity of the situation. However, we believe that MCMC methods, like those presented here, do give a realistic strategy for identifying and characterizing the large number of signals, of all types, that will exist in LISA data.
Besides LISA, the methods discussed here are likely to be useful in other fields of study where the data contain an unknown number of periodic signals. One timely area of astrophysics is in the field of asteroseismology, where attempts are made to measure vibrational modes of a star; see for example Handler et al. [23]. We would expect that our method could be directly applied to asteroseismic data on stellar oscillations. Similarly, spectroscopic nuclear magnetic resonance studies are often concerned with spectral decomposition into an unknown number of components, dependent on the composition of the sample [24], and we believe the parameter estimation technique we present could be usefully applied here too.