{'fraction_non_alphanumeric': 0.032580197409006785, 'fraction_numerical': 0.011335595311536089, 'mean_word_length': 4.502121340687315, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Studying stellar binary systems with the Laser Interferometer Space Antenna using Delayed Rejection Markov chain Monte Carlo methods

Bayesian analysis of LISA data sets based on Markov chain Monte Carlo methods has been shown to be a challenging problem, in part due to the complicated structure of the likelihood function consisting of several isolated local maxima that dramatically reduces the efficiency of the sampling techniques. Here we introduce a new fully Markovian algorithm, a Delayed Rejection Metropolis-Hastings Markov chain Monte Carlo method, to efficiently explore these kind of structures and we demonstrate its performance on selected LISA data sets containing a known number of stellar-mass binary signals embedded in Gaussian stationary noise.


Introduction
Surveys of our galaxy in the gravitational-wave band with the Laser Interferometer Space Antenna, LISA [1] will yield the largest sample of close stellar-mass binary systems, in particular white dwarfs. LISA is in fact expected to individually resolve ∼ 10 4 stellar-mass binary systems that are sufficiently bright to stand above the confusion noise produced by the whole population of binaries, most of which are unresolved [2,3]. Gravitational waves from stellar-mass compact objects in the lowfrequency observational window are long lived, and in order to extract the maximum amount of information, one needs to be able to disentangle and accurately measure the parameters of a large number of sources overlapping in time and frequency space. This represents a challenging data analysis problem, whose successful solution is essential for the full science exploitation of the LISA data set.
Due to the large number of sources to be analysed at the same time and the fact that this number is unknown, Markov chain Monte Carlo (MCMC) techniques [4,5,6,7,8,9] and their extension to Reversible Jump Markov chain Monte Carlo (RJMCMCs) [10,11,12] are expected to be one of the most powerful search methods, although analysis approaches based on matched-filtering with a template bank have also been explored [13,14].
In MCMC applications to LISA data analysis, one of the key problems that one needs to tackle is the fact that the joint posterior probability density function (PDF) of the source parameters -the target distribution -presents a multimodal structure, often with secondary maxima well separated from the mode, see Figure 1 for an example. In such situations, maintaining a high acceptance ratio while sampling the full structure of the posterior PDF becomes very difficult, and may lead to a negligible efficiency of the algorithm. Several solutions have been proposed to increase the mixing and exploration ability of the chains, such as using parallel chains exploring the parameter space that, at a certain point, can be swapped (e.g. Metropoliscoupled MCMC algorithms [15]), or other methods like simulated annealing [16] or simulated tempering [17,18] that consist in adding a 'temperature' factor to the target distribution in order to make it smoother at the beginning and therefore to promote the movement of the chain. Examples of the use of some of these techniques to compute the posterior PDF for white dwarf binaries in LISA data sets are given in e.g. Refs. [4,5,6,11,12]. However, as the number of dimensions increases and the target distribution becomes more complicated, it is not clear which technique works best. For instance, running N parallel chains requires N times more computing time; moreover, it has been shown that in some cases, this approach may yield low efficiency [19]. Non-Markovian explorations during an initial "search phase" have been considered to quickly identify the neighbourhood of the main mode of the distribution; however these strategies suffer from the significant draw-back that they cannot be readily extended to the relevant case of an unknown number of dimensions (i.e. unknown number of signals) and as a consequence model selection. For all these reasons, it is useful to explore alternative approaches.
In this paper we explore a different strategy for a fully Markovian sampling of multimodal posterior PDFs: a Delayed Rejection (DR) Metropolis-Hastings Markov chain Monte Carlo method. Delayed Rejection was originally introduced by Tierney and Mira [20,21] for fixed dimension problems and then extended by Green and Mira [22] to transdimensional problems. The algorithm is now found in a wide variety of science applications, see e.g. Refs. [23,24,25,26,27,28]. In general, in all applications so far, either the number of steps was very small (2 or 3) or the assumption of symmetrical proposals is made. We have recently extended the DR algorithm to its most general scheme, using an arbitrary number of stages [29]. Here we demonstrate its power in applications to the LISA white dwarf problem, by applying it to the analysis of data sets in the restricted case in which the dimensionality of the problem (that is proportional to the number of white dwarf presents in the data set) is fixed and known.
We would like to emphasise that although the waveforms that we consider here are the simplest present in the LISA data set, this specific example tackles all the conceptual issues that will be present in other LISA data analysis problems involving more complex signals, such as those from spinning massive black hole binaries and extreme-mass ratio inspirals. The DR scheme will therefore be applicable to the development of MCMC analysis algorithms targeted at the study of those sources.

The general method
The Metropolis-Hastings algorithm is widely used in applications of Markov chain Monte Carlo methods to sample a target distribution π(x), in our case the posterior PDF of a set of parameters. Given a state of the Markov chain λ, one proposes a new state λ ′ drawn from a proposal distribution (or transition kernel) q( λ, λ ′ ), the probability of λ ′ given λ. This new state is accepted with probability and the chain remains at λ with probability 1 − α( λ, λ ′ ). In the previous equation, the notation a ∧ b (for any real numbers a and b) stands for the minimum between a and b, and as a consequence the acceptance probability is limited to unity in the case where the ratio is > 1.
An effective way to improve the efficiency of the MCMC algorithm has been proposed in the form of Delayed Rejection [20,21,22], and now we briefly summarise it. Suppose that the current state of the chain is λ. A candidate, β 1 , is generated from a proposal distribution q 1 ( λ, β 1 ) and accepted with probability as in a standard Metropolis-Hastings algorithm, see Equation (1). If β 1 is not accepted, instead of rejecting it, one can use this information and propose a new state β 2 from a (in principle different) proposal distribution q 2 ( λ, β 1 , β 2 ) which may use information about the rejected state β 1 . Moreover, in general the i−th stage of the DR chain is as follows [29]: if the state β i−1 is proposed and rejected, one can propose a new candidate β i from the proposal q i ( λ, β 1 , . . . , β i ) and accept it with probability Notice that with the notation we are using, q i ( λ, β 1 , . . . , β i−1 , β i ) represents the proposal probability of β i given { λ, β 1 , . . . , β i−1 } where the order of the parameters does matter. By imposing detailed balance at each stage and deriving the acceptance probability that preserves it, the resulting chain will be a reversible Markov chain with invariant distribution π. From a theoretical point of view, it can be demonstrated in general (see Sec. 2.2 of Ref. [29]) that the variance of an estimate made from a chain using Delayed Rejection is always smaller than that produced with a standard Markov chain. Moreover, besides increasing the efficiency of the sampling, this method allows us to propose a new candidate using the information of the past elements of the chain. We will profit from these two properties in increasing the efficiency of the computation of the marginalised posterior distributions generated by integrating multimodal likelihood functions.

Key rules on the choice of proposals
An essential point of the DR scheme is that Eq. (3) has been derived by imposing that the backward path, from β i to λ, follows the forward path, from λ to β i , in reverse order. As a consequence, this property forces one to choose proposal distributions that preserve, in some way, the reversibility of the chain in order to avoid negligible acceptance probabilities. Full details about the general implementation of the algorithm are provided in Ref. [29]; here we focus on the choice of proposals relevant to searches for white dwarf binary systems.
When we run a Metropolis-Hastings MCMC chain -likely exploring a secondary mode of the distribution -we can start a DR chain with a large jump that attempts to reach the mode of the distribution; such transition, likely to be rejected but hopefully closer to the mode, is followed by many small ones in order to explore the new region of the parameter space. In terms of the proposal PDFs that appear in Eq. (3), this translates into having q 1 ( λ, β 1 ) ≡ q a ( λ, β 1 ) mainly proposing large transitions in the parameter space, and q j ( λ, proposing small transitions. This choice inevitably adds some degree of non-reversibility that significantly degrade the acceptance probabilities, see Eq. (3). One could however achieve an efficient DR algorithm, by using proposals that allow for both small and big jumps, such as a mixture of three Gaussians (3-Gaussian, from now on) symmetrical functions: in the previous equation,x may be a function of the old elements of the DR chain, N a,b represents the probability of proposing a value from the central Gaussian and . Notice that an efficient exploration of multimodal structures could in principle be achieved by taking N a → 0 and N b → 1, but this in general would lead to negligible acceptance probabilities. In addition, the central value of the 3-Gaussian,x[ λ, . . . , β j−1 ] that is computed at each DR stage making use of the information from the old elements of the chain, should be independent of the particular sample of the chain. In our application we use the mean of all the past elements of the chain after the big jump, With these choices, all the proposals used during the DR are characterised by 6 parameters: {σ 1 , σ 2 , µ, N a , N b } and the maximum number of elements in a DR chain, n dr ; as a consequence, the contribution of the ratio of proposals to the final acceptance probability, Eq. (3), can be quantified (see Figs. 4 and 6 of [29]) in terms of these parameters. The choice of these parameters is the essential step in obtaining a DR algorithm that probes in an efficient way the posterior PDF. In practice, they can be chosen as follows: (i) {σ 1 , σ 2 , µ}, are chosen to have a proposal adapted to the multimodal structure of the target distribution, µ being the typical distance between two maxima of the target distribution and σ i their typical width; as consequence a qualitative understanding of the typical scales in parameter space of the likelihood function is important (we deal in detail with this point in the next Section); (ii) from {σ 1 , σ 2 , µ}, one can then identify the panels from Figs. 4 and 6 of [29] corresponding to the closest σ 1 /µ − σ 2 /µ values; first, one looks at the losses due to the central proposal evolution (CPE) (Figure 6 of [29]) to select N b high enough ‡ to have acceptable expected values of this ratio of proposals; as one normally works with n dr 100, the results depend very weakly on n dr ; (iii) one then needs to consider the losses due to an asymmetrical proposal (AP) (Figure 4 of [29]) in order to choose N a ; (iv) finally, as again the dependency on n dr is weak, its choice will be mainly driven by the constraints on the computational costs.
So, by following the prescriptions described here and having some previous knowledge about the multimodal structure of the target distributions that has to be sampled, one can implement a DR MCMC algorithm with an arbitrary numer of steps, n dr , capable of efficiently sampling these multimodal distributions. The algorithm is completely general, and can be applied to a wide range of problems. In the following sections we will focus on the application of the DR MCMC algorithm to the search of a known number of stellar-mass binary signals in Gaussian and stationary noise with LISA.

Application to LISA observations of white dwarfs
We consider a LISA data set -that we indicate with d A,E , to highlight the use of the two noise-orthogonal Time-Delay-Interferometry observables A and E, see e.g. Ref. [30] for a review -containing a known number, N , of signals from stellar-mass binary systems (h A,E ) buried in Gaussian stationary noise (n A,E ): The vector ζ = { ζ n , ζ s } contains all the M unknown parameters that describe the problem, where ζ n and ζ s are the noise and signal parameters, respectively. For the specific case considered here, we focus on monochromatic (in the source reference frame) white dwarfs and purely Gaussian and stationary (with zero mean and variance σ 2 n ) instrumental noise. The total number of parameters that describe the problemi.e. the dimensionality of ζ -is M = 1 + 7N where f 0 is the constant frequency of the signal (with respect to the Solar System Barycentre), A + and A × are the amplitudes of the two GW polarizations, the longitude φ and co-latitude θ are the two angles that define the source's sky location, ψ is the polarization angle of the gravitational wave and ϕ 0 is a constant that fixes the initial phase of the signal. We make use of the Delayed Rejection MCMC algorithm introduced in the previous Section (see [29] for more details) to compute the (marginalised) posterior PDFs of the noise and signal parameters. By applying the Bayes' theorem, we can compute the joint posterior density function, p( ζ|d), of ζ given the data sets, d, where L(d| ζ) is the likelihood function, p( ζ) the prior probability density function, and p(d) is the marginal likelihood. From a Markov chain Monte Carlo algorithm perspective, the target distribution that needs to be sampled is π( ζ) ≡ p( ζ) L(d| ζ) ∝ p( ζ|d). Since the resulting chain is a Markov chain with invariant distribution π, the average along a sample path of a function f is an asymptotically unbiased estimate of f (x)π(x)dx; in particular, the marginalised posterior PDF over a set of parameters, can be directly extracted from the output Markov chain of the algorithm.

Structure of a single source target distribution
The application of the DR algorithm requires some previous knowledge of the multimodal structure present in the target distribution of the problem at hand. This structure is intrinsic to the single source posterior PDF and, in particular, to its likelihood function, L(d| ζ). Local maxima [6] occur in the likelihood surface at approximately multiples of the modulation frequency (f m = 1 yr −1 ) around the frequency of the signal measured at the Solar System Barycentre. Although the main multimodal structure occurs in the frequency parameter, other parameters, like the sky position, are correlated to the frequency, and this is reflected in the structure of secondary maxima of the likelihood function.
In Fig. 1 we explore the structure of the likelihood function in the 3-dimensional space {f, θ, φ}, around the global maximum and the secondary mode, for three different frequency values. It can be observed how local modes of the maximised likelihood function appear around the actual frequency value at distances that are multiples of f m and, at the same time, the sky location is affected. As the signal frequency increases, a smaller change in the sky location can produce the same frequency shift due to the Doppler effect produced by the LISA motion, which explains the finer structures in the panels at the bottom of Fig. 1. In order to have a more systematic and quantitative understanding of the multimodal structure of the likelihood function, we have performed 500 Monte Carlo (MC) simulations around different frequencies, taking random values for all the 7 parameters that characterise a single white dwarf binary signal, and have computed the relative position (in parameter space) of the secondary maxima with respect to the main mode. The results of this Monte Carlo's are summarised in Fig. 2; from those we can make a suitable choice of the parameters {σ 1 , σ 2 , µ} that determine the proposals of the DR stage.
Some general conclusions, that are important for the tuning of the DR algorithm can be drawn. The frequency separation between the two modes appears to be independent of the actual frequency of the signal and the values are mostly concentrated between 0.75f m and 1.75f m . Moreover, we can estimate the angular distance on the two-sphere, ∆Ω, between the sky position corresponding to the main and the secondary maximum. From the central plot of Fig. 2 we observe that ∆Ω progressively shifts to lower values as the signal's frequency increases, and that (righthand plot of Fig. 2) the mean value of ∆Ω changes with the co-latitude, θ, producing secondary maxima that are more separated near the ecliptic plane than at the poles. A qualitative understanding of the dependency observed in Figure 2 can be simply derived from the expression of the Doppler shift,φ D, wherek is the unit vector from the source to the detector and v i is the velocity vector of the i-th LISA spacecraft. From here, we observe that the required changes in the sky location needed to obtain a certain frequency shift, will be inversely proportional to f 0 .

Parameters for the Delayed Rejection proposal distributions
The multimodal structure of the target distribution relevant for a search for stellarmass binaries in the LISA data is mainly present in the frequency of the signal, although there are correlations with the sky location angles. The solution that we have adopted here, as a first implementation of the DR algorithm, consists in applying a scheme where only these three parameters, {f, θ, φ}, are updated: the first one, f , using a mixture of three Gaussian distributions and the two angles defining the sky location drawn from a single Gaussian proposal centered at the previous element of the DR chain and width according with the results obtained in Fig. 2. With this, only the frequency proposals can contribute with small factors to the acceptance probability of Eq. (3), so we will now focus on these proposals.
By following the prescription given in Sec. 2.2, we select the first three parameters, {σ 1 , σ 2 , µ}, that characterise the 3-Gaussian distribution for the frequency. In light of the results shown in Figs. 1 and 2, our proposals can have the modes separated by µ = 1.25f m , choosing a σ 1 = f m /2 wide central Gaussian when exploring local modes and σ 2 = 0.2f m for the external Gaussians. This means that we propose from the external Gaussians, 95.4% of the values in the (0.85f m , 1.65f m ) interval.
Moreover, by looking at the corresponding panels of Figures 6 and 4 of [29] (in our case, σ 1 /µ = 0.4 and σ 2 /µ = 0.16), one can read out the values for N b and N a . Finally, the maximum number of iterations allowed in a Delayed Rejection is fixed mainly for computational reasons; from our experience we think that n dr = 1000 to 3000 are adequate values, depending on whether one prefers to attempt more frequent short DR chains, or longer chains less frequently. In summary, a DR algorithm for efficiently explore the joint posterior PDF in the frequency space can be defined by the following parameters:    Fig. 3 but here working with a 5 sources data set around 1.6mHz, with a separation of 13fm between each source. Only the third chain (blue in the on-line version) converged to a wrong sky location (the opposite to the actual one) because of degeneracies present in the likelihood surface and its proximity to one of the poles.

Results
We have implemented a DR MCMC algorithm to search for a known number of stellarmass binary signals in LISA data, assuming Gaussian and stationary noise. In this section we provide example results obtained by applying the analysis method on two different 1-yr long data sets: (i) one containing three sources at 1.1 mHz separated in frequency by 2f m and 40f m , see Fig. 3, and (ii) the other including 5 sources around 1.6 mHz equally separated in frequency by 13f m , see Fig. 4. From these figures we can see that after the burn-in period -we note that the DR part is switched on only after burn-in -that has been discarded, most of the chains are sampling either the actual mode of the signals or secondary maxima. At this point, the Delayed Rejection is switched on (with probability 1/3000) and from time to time a new region ≈ 1f m away from the current position of the chain, is explored searching for higher maxima. We observe that once the chain finds the right mode in the frequency space, all the other parameters also move their actual values.
With the DR, we are drawing n dr = 1000 proposals mainly targeted to a different mode of the posterior PDF with a probability of 1/3000, most of which are rejected. By keeping the same odds with a standard MCMC, we would have the same ability of exploring different maxima but this algorithm would yield an extremely correlated Markov chain, i.e. greater variances of the estimates once it has reached the stationary distribution. On the other hand, one could use a standard MCMC proposing big jumps in the parameter space with probability 1/3000, and the resulting Markov chain would have the same autocorrelation as the one obtained from a DR MCMC; however, the ability to explore different maxima of the target distribution would decrease by a factor n dr = 1000. We refer the reader to Ref. [29] for further discussions and examples.
The results of Figs. 3 and 4 were obtained from a single Markov chain of 5×10 5 ×N and 15 × 10 5 × N elements, respectively; where the first 25000 × N elements were discarded as a part of the burn-in period. The length of the chains is sufficient to demonstrate the ability of the DR MCMC algorithm to find multiple sources, but in a real problem, we would run longer chains in order to get accurate samplings of the posterior PDFs.
In addition of the DR parameters, see Sec. 3.2, many others for the MCMC algorithms need to be chosen carefully. Here we summarise the most relevant for the present discussion. First of all, we are not always updating the same number of sources. We set probabilities of updating 1, 2, . . . , N sources, always selecting the closest ones in frequency, since this is the main parameter that determines the degree of correlation between two sources: two sources with exactly the same parameters but separated by, say 1000f m , are perfectly resolvable. Our proposal density functions during the normal MCMC evolution are single Gaussian distributions centered at the actual element of the chain, with constant values for the amplitude of the proposals, ∆ i , throughout the evolution of the chain; for the angles, the value of ∆ i was set to one-ninth of the prior range (which covered the entire physical space). For the other parameters we use ∆ f = 1f m = 3.17nHz, ∆ A×/A+ = 2/3, ∆ log A+ = 0.08 and ∆ log σ = 0.01. We also consider updates of all the parameters at a time, but always working with diagonal variance-covariance matrices for the proposals. When the number of sources is greater than 1, the DR stage is attempted with all the selected sources at a time. We account for symmetries and degeneracies of the target distribution by proposing, with a given probability, special jumps based on our prior knowledge. Most of these 'quasi-equivalent' points in the parameter space are actual degeneracies of the antenna patterns that characterise the LISA response in the low frequency range [31], and that at higher frequencies appear as secondary maxima. When searching for more than one source, we also consider big jumps in the parameter space: single Gaussian proposals of width equal to one-third of the whole prior range for all the parameters. We found them necessary (see e.g. Fig. 4) because converged chains act like "walls" in the frequency space, decreasing the chances of other chains to cross over them.

Conclusions and future work
The motion of the LISA instrument during the observations of long lived GWs modulates the phase and amplitude of the detected signal, producing likelihood functions with a complicated structure consisting of several isolated local maxima. This is a common feature of all long lived LISA sources and, in turn, makes standard MCMC samplers inefficient. Here, we have presented a new strategy for the fully Markovian sampling of multimodal posterior PDFs based on a Delayed Rejection MCMC method and we have demonstrated its performance by applying it on the search for multiple stellar binary systems in Gaussian and stationary noise.
The Delayed Rejection technique introduced here, can be fruitfully applied in many other LISA data analysis areas involving more complex signals, such as those from spinning massive black hole binaries and extreme-mass ratio inspirals; furthermore, it needs to be extended to the case of trans-dimensional problems, relevant to the analysis of LISA data sets to search for the whole galactic populations of stellar-mass binary systems, where the number of GW sources in unknown. In fact, an additional benefit of having a fully Markovian algorithm capable of efficiently exploring the target distribution is the possibility of estimating the Bayes factor of different models directly from the output chains. We are currently exploring the possibility of combining the Delayed Rejection technique with a Reversible Jump MCMC algorithm in order to search for an unknown number of sources.