Bayesian Model Selection for LISA Pathfinder
Abstract
The main goal of the LISA Pathfinder (LPF) mission is to fully characterize the acceleration noise models and to test key technologies for future spacebased gravitationalwave observatories similar to the eLISA concept. The data analysis team has developed complex threedimensional models of the LISA Technology Package (LTP) experiment onboard LPF. These models are used for simulations, but more importantly, they will be used for parameter estimation purposes during flight operations. One of the tasks of the data analysis team is to identify the physical effects that contribute significantly to the properties of the instrument noise. A way of approaching this problem is to recover the essential parameters of a LTP model fitting the data. Thus, we want to define the simplest model that efficiently explains the observations. To do so, adopting a Bayesian framework, one has to estimate the socalled Bayes Factor between two competing models. In our analysis, we use three main different methods to estimate it: the Reversible Jump Markov Chain Monte Carlo method, the Schwarz criterion, and the Laplace approximation. They are applied to simulated LPF experiments where the most probable LTP model that explains the observations is recovered. The same type of analysis presented in this paper is expected to be followed during flight operations. Moreover, the correlation of the output of the aforementioned methods with the design of the experiment is explored.
I Introduction
LISA Pathfinder (LPF) Antonucci et al. (2012a); McNamara et al. (2013) is an European Space Agency mission that will serve as a technology demonstrator for a future spacebased gravitationalwave observatory like eLISA Seoane et al. (2013). The LPF mission will prove geodesic motion by monitoring the relative acceleration of two test masses in nominally freefall in the frequency band of 1 to 30 mHz. The main instrument onboard LPF is the LISA Technology Package (LTP), which is a suite of experiments with the aim of measuring and characterizing the different contributions to the differential acceleration noise between the two test masses. This characterization is the main task of the data analysis team. To that end, dedicated experiments are going to be performed in order to estimate the unknown parameters of the system. And for that purpose, a number of parameter estimation methods and models of the LTP have been implemented Nofrarias et al. (2010); Congedo et al. (2012); Antonucci et al. (2011). The main question that arises is about the suitability of the different models implemented, or in simpler terms, which model can describe better the observations of the experiment.
The motivation to implement an algorithm that will help us classify our models is due to the nature of the LTP system. There has been a lot of work trying to understand the instrument, and the models implemented are based on theoretical and experimental measurements from test campaigns Cervantes et al. (2013); Audley et al. (2011). Selecting the most probable model is crucial for the analysis for two main reasons: first, the most suitable LTP model that describes the observed physical effects is chosen, and secondly, overfitting situations are avoided, together with biased estimation of the parameters of the system.
One could use several criteria and algorithms that classify competing models, but in the end, working in a Bayesian framework, the main aim is to calculate the Bayes Factor, which is a comparison between the evidences of the models Kass and Raferty (1995); MacKay (2003). The evidence is defined as the probability of the data given the model, that is, the probability distribution on that quantifies the predictive capabilities of the given model.
Most of the approaches are based upon the likelihood evaluated at the maximum and a penalty for the number of parameters in the model consisting in multiplying by the socalled Occam Factor. The Occam’s razor is the principle that states that the simplest hypothesis that explains the observations is the most favorable. In our case we assume that we have to compare two different LTP models; model and a simpler one, , over a data set . If is a more complex model, it presents more predictive capabilities than , which translates to more disperse evidence over the data set . Thus, in the case where the data is compatible with both models, will turn out to be more probable than , without having to express any subjective dislike for complex models MacKay (2003). In other words, it is almost certain that the more parameters in a model the better the fit to the data. But taking it to the extreme, we can imagine a model with as many parameters as the data to fit. In that case we have overparameterized the model while the aim is to include only the parameters which substantially improve it.
In this paper we investigate several methods to compare competing LTP models giving emphasis to Reversible Jump Markov Chain Monte Carlo techniques. The plan of the paper is as follows: In Sec. II, we make a brief introduction to Bayesian methods and in Sec. III, the LTP experiments and models are thoroughly explained and we investigate several applications of our Bayesian techniques to LTP experiments. In Sec. IV, we compare the available methods and discuss their output in connection with the experiment design. We end with a summary and conclusions in Sec. V. The codes for the different algorithms are integrated into the LISA Technology Package Data Analysis (LTPDA) Matlab toolbox Hewitson et al. (2009); web (2005), that comes together with proper documentation and help for the user.
Ii Background on Bayesian Statistics
An algorithm that automatically penalizes higherdimensional models is the Reverse Jump Markov Chain Monte Carlo (RJMCMC) algorithm. The RJMCMC method Green (1995); Lopes and West (2004); Dellaportas et al. (2002); Stroeer and Veitch (2009); Umstätter et al. (2005) is widely used when dealing with nested models, meaning that we need to compare a set of models, where simpler models are a subset of a more complicated one. In fact, the RJMCMC algorithm is the generalized case of Markov Chain Monte Carlo (MCMC) methods that is capable of sampling the parameter space and at the same time jumping between models with different dimensionality.
At this point it would be convenient to describe the philosophy of the Bayesian framework and the Metropolis algorithm and then move to the general case of the RJMCMC algorithm. The Bayes rule can be summarized by the following equation:
(1) 
where is a model parameter vector, is the likelihood of the parameters over the dataset , and and are the posterior and the prior distributions of the parameters respectively. Note that the evidence , or marginal likelihood, is often neglected in parameter estimation algorithms, as it serves only as a normalisation constant:
(2) 
The Metropolis algorithm is one of the MCMCtype methods available that is widely used for parameter estimation purposes. It is based on sampling the parameter space by proposing new samples and evaluating the likelihood at each step. By sampling the posterior distribution, probability distributions to the parameters to be estimated are assigned. It is certain that given a large amount of steps in the parameter space, the Metropolis algorithm will converge to the set of parameters that maximise the likelihood.
ii.1 Calculating the Bayes Factor
The evidence of a hypothesis given the dataset , , states the support for this hypothesis, or in other words, how much the data favors a given model. In our case, the hypothesis is the model implemented to describe the LTP system. Consequently, given two different models X and Y, the Bayes Factor is a comparison between the evidences of model and model given by their ratio
(3) 
where
(4) 
This integral is extremely costly to evaluate, specially when the model becomes complicated with higher dimensionality. In the next sections, a selection of estimators of the Bayes Factor are described. Most of them evaluate directly the Bayes factor, while other techniques are used to estimate the evidence for each model. All of them are used to investigate our LPF models and mission experiments. In the end, we can draw conclusions about the competing models given the estimated value of the Bayes Factor. If , the evidence is negative and the observations support model Y. If , the evidence is positive and model X is more favourable than model Y.
ii.2 The Reversible Jump Markov Chain Monte Carlo Algorithm
The RJMCMC algorithm is a robust and efficient tool to estimate the Bayes Factor. It can be shown Green (1995); Dellaportas et al. (2002) that after a large number of iterations it will converge to the true value of the Bayes Factor. The only drawback is the computational cost of the algorithm. When more than three models are being compared, meaning that many transdimensional moves have to be performed, a considerable amount of time is required for convergence. The algorithm implemented in this work is a special case of the Metropolized Carlin and Chib method Lopes and West (2004). More specifically, let us suppose that we have a total number K of models to compare given a data set . Then, the recipe for our RJMCMC method can be summarised in the following steps:

Initialization: Choose an initial model and the corresponding parameters .

Apply the Metropolis algorithm for model . This step is also called the “in model step”.

Generate new from a multivariate Gaussian PDF and a random number from a uniform distribution. This is the step where we propose a new model .

Calculate the acceptance ratio :
(5) where is the proposal distribution from where the “dimension matching” parameters are drawn Godsill (2001), and is the Jacobian:
(6) 
If we accept the new model with parameters and set .

Iterate from step 2 until convergence.
The old set of parameters is connected to the new one by a well defined function (and of course ). We use independent proposals, so and , thus, the Jacobian term in equation (5) is unity. The algorithm spends most of the time iterating “inside” the model that best describes the data. The RJMCMC method autopenalizes high dimension models, also by taking into account the priors of each model . They serve as an Occam Factor integrated within the algorithm. Convergence is achieved if two main conditions are satisfied. First, the condition of reversibility, which is stated in a simple way: The proposal function must be invertible, meaning that we can jump from the proposed parameters back to the current parameters. And second, we must satisfy the dimension matching condition which in our case is always true since we use independent proposals in the acceptance ratio. After convergence has been achieved, a good approximation to the Bayes Factor is given by Littenberg and Cornish (2009); Bartolucci et al. (2006); Dellaportas et al. (2002); Cornish and Littenberg (2007)
(7) 
ii.3 Other approximations
While the RJMCMC method directly estimates the Bayes Factor, the other methods implemented in our work make an approximation to the evidence of each model. We implemented them as a crosscheck for the RJMCMC method, but they also provide certain freedom of choice, depending on the nature of the problem, as well as the data available.
The first approximations we consider are the Laplace approximations. The Laplace approximations perform a comparison between the volume of the models in the parameter space and the volume of the uncertainty ellipsoid of the parameters Kass and Raferty (1995). This comparison is feasible if we assume that we work within a high SignaltoNoise Ratio (SNR) and therefore the posterior PDF is Gaussian near the maximum likelihood parameters, . Another requirement is that a sufficient large number of samples must have been collected. Then, the evidence in Eq. (4) becomes:
(8) 
where is the number of dimensions of model X and is the Hessian matrix of the posterior . There are two main variations of the Laplace approximation. In the first one we make use of the Fisher Information Matrix , calculated at , as an approximation to the expected covariance matrix Cornish and Littenberg (2007); Vallisneri (2008). Then, the evidence of the model becomes:
(9) 
Of course, the main limitations of this method are associated with the confidence we have on the calculation of the Fisher Matrix. Furthermore, as expected, the results appear to be poorer in comparison with the other methods as we move towards lower SNR areas. We can follow the notation of Cornish and Littenberg (2007) and call this particular approximation the LaplaceFisher (LF) approximation. Another wellknown variation is the LaplaceMetropolis (LM) estimator of the marginal likelihood Kass and Raferty (1995). In this case, we use all necessary components for the calculation of the evidence from previous MCMC estimates. The parameters are extracted from the chains of a MCMC parameter estimation run for the particular model, while we use the weighted covariance matrix of the chains , using a Minimum Volume Ellipsoid (MVE) or a Minimum Covariance Determinant estimator (MCD) Rousseeuw and Leroy Annick (1987). The MVE method has been implemented and integrated in the LTPDA toolbox. In this case, the evidence of model X can be written as
(10) 
The LM method is considered to be a very reliable tool for the computation of the evidence of a model Kass and Raferty (1995). The third method we have used is the SchwarzBayes Information Criterion (SBIC) and is based on the following: After the assumption that the priors for each model follow a multivariate Gaussian PDF, the Schwarz criterion is defined as:
(11) 
where and are the dimensions of each model and is the number of samples in the data. It can be proven that if , then and the Schwarz criterion can be a good approximation to the Bayes factor. In fact, must be chosen carefully so that , where is the number of effective samples in the data that represent the growth of the Hessian matrix of the loglikelihood Kass and Raferty (1995).
ii.4 Implementation
If the system’s transfer function is and we consider known injection signals, , then the output of the system, , can be expressed as:
(12) 
where is the overall instrument noise and is the response or template of the system. The noise in the LTP experiment can be approximated with Gaussian noise since we work in a high SNR regime and we can safely consider that any dependence on the system parameters is negligible. Consequently, the likelihood for a model with parameters can be written as
(13) 
where the angular brackets denote the noise weighted inner product Barack and Cutler (2004)
(14) 
and where is the power spectral density of the noise. For the proposal distribution appearing in the acceptance ratio of the MetropolisHastings algorithm a multivariate Gaussian distribution is used. Although the choice of the proposal distribution should not affect the final estimated parameter PDFs, it greatly affects the convergence speed of the chains to the target posterior distribution. The covariance of the multivariate Gaussian is calculated beforehand or during the search phase using the Fisher Information Matrix (FIM) Vallisneri (2008). In the case of the RJMCMC method, the FIM is calculated once for each model before the execution of the algorithm. It is computed numerically or analytically depending on the model. The LTP models are coded either in StateSpace format Nofrarias et al. (2013); Hewitson et al. (2011); Grenadier and Weyrich (2008), or in pure analytical format in the acceleration domain Congedo et al. (2013).
The first step in the analysis is to apply a Fourier transformation to the data and estimate the power spectral density of the noise , i.e. . As described in Nofrarias et al. (2010), the nature of the LTP experiment allows us to combine all the available information from different investigations. This means that we can combine the posterior distribution for each data set for each model. Then, we can apply MCMC analysis strategies to sample the posteriors and make use of any of the criteria described in section II.3 or sample the joint parameter spaces with the RJMCMC method that was described in section II.2. Technically, since the RJMCMC algorithm can be seen as a generalised MCMC method, it performs a parameter search by simultaneously sampling the likelihood and consequently it maps the joint posterior distribution and assigns PDFs to the parameters. The extra point is that we can estimate directly the ratio of the evidences of the models and test the hypothesis made for the data.
Iii Modeling the LTP experiment
As a technology demonstrator of a spacebased gravitationalwave observatory,
the LPF mission will place two kg Test Masses ( and ) in free fall.
The goal is to estimate the residual differential acceleration between and Antonucci et al. (2011).
To minimize all external forces acting on along the xaxis, a
dragfree control loop Fichter et al. (2005) has been designed. The coordinates of are controlled via electrostatic
actuators and the Spacecraft (SC) is controlled by microNewton thrusters. More specifically, the main components of the LPF mission are:

The Gravitational Reference Sensor (GRS) Dolesi et al. (2003). It consists of the testmasses and the vacuum chamber around them. It is mounted with two identical electrostatic actuators to control the 6 degrees of freedom of each testmass, and also to apply forces and torques to keep the testmasses in free fall.

The Optical Metrology System Heinzel et al. (2005) consists of the optical bench, its subsystems and the processing computer. It performs the sensitive optical measurements of the positions of the testmasses along the xaxis. For the simplified onedimensional model version, we consider two interferometer inputs and two outputs of the system. The measured displacements are and , the distance of to the SC and the distance between and respectively.

The current design of the propulsion system for LPF is based on Cold Gas McNamara et al. (2013) thrusters. The main function of the thrusters is to maintain the reference testmasses in free fall conditions in the measured bandwidth.

Finally, the DragFree and Attitude Control System (DFACS) calculates all forces and torques acting on the SC and testmasses and computes the commanded forces/response of the system in order to maintain the nominal free fall of .
The LTP is a complicated instrument composed by a plethora of subsystems and coupled control loops. During the mission, there will be dedicated experiments with the aim of characterizing the different noise contributions and couplings. In the following section, the LTP models implemented for system identification experiments are described. We present simulations of the planned system identification experiments and use all techniques described to perform model selection.
iii.1 The LTP system dynamics
For the purposes of this work we consider a simplified version of the LTP operating in the socalled main science mode and we confine the system dynamics to the onedimensional case. There are two approaches in modeling the dynamics of the threebody system: The first one, shown in Fig. 1, is to fully describe the dynamics together with the controllers in a statespace format (SSM). In this case, it is very convenient to represent the model as a closed loop system, due to closedloop dynamics and controllers of the instrument. The white boxes represent standalone SSMs, all together assembled to form the final LTP experiment Nofrarias et al. (2013). Each submodule has injection and output ports and can be represented as:
(15) 
where are the states of the system, A is the state matrix, B is the input matrix, C is the output matrix, and D is the feedthrough matrix Nofrarias et al. (2013). This form of implementation has numerous advantages, like modularity and flexibility in modeling the LPF mission. For example, the user is able to combine any given SSM module and use a custom noise model for each particular subsystem.
The second approach is to represent the dynamics in an analytical way in the acceleration domain Congedo et al. (2013) which is described in section III.3. A more indepth investigation on the analysis in acceleration domain and its practical advantages is to appear soon.
For the purposes of this paper, we simulate the experiments as explained in Nofrarias et al. (2010); Congedo et al. (2012); Nofrarias et al. (2013, 2012); Vitale (2007). Then, the complete system is controlled by optical readouts that measure the positions of the Test Masses, and , where the axis is defined by the line joining them. In order to perform parameter estimation we inject sinusoidal signals to each interferometric channel alternately, or direct forces to the three bodies of the system (testmasses and SC), and measure the response of the system. For this type of experiments, we can define the following timeseries arranged as vectors:
(16) 
(17) 
where is the measured signals vector, is the interferometer injection signals vector, is the commanded forces signals vector, and is the overall noise for both channels. Here and are the commanded forces applied to and respectively while is the force applied to the SC. The noise contributions of each LTP subsystem are denoted in Fig. 1 as . It can be instrumental or readout noise, or it may originate from external sources (like solar radiation).
For the sake of simplicity, we use the onedimensional version of the models. The set of parameters that we recover from the system is : and stand for the electrostatic stiffness of the testmasses to the surroundings, and are time delays, and the capacitance and thruster gains respectively, and denotes the crosscoupling between the two interferometric channels.
As described in Vitale (2007); Nofrarias et al. (2012), during the planned mission experiments, sinusoidal signals of different frequencies are going to be injected to the control loop, simulating displacements of the test masses. Besides the LTP dynamics system, the data analysis team has modeled noise sources for the various subsystems of the LPF mission based on theoretical predictions or characterization of each element from onground test campaigns. The main noise contributions come from the thrusters, the interferometric readouts, and the capacitance actuators Antonucci et al. (2012b). We can simulate synthetic noise for any particular experiment and perform parameter estimation exercises in order to test the system identification algorithm’s readiness. A MCMC search of the parameter space for such a simulated run returns satisfactory results as shown in Table 1.
Parameter  Real value  Estimated  

0.82  
1.08  
0.2  
0.2  
0.0004 
By using the LTPDA machinery, synthetic noise can be generated and the response of a known system can be simulated. For the following study cases, the LTP default parameters were always set to the same values given in Table 1. It should be noted that there is no reason to have any strong prior belief about the true values of the parameters. Therefore, it is reasonable to assume uniform priors for all the tests.
iii.2 Application to a simplified LTP model
In order to demonstrate a first application to the LTP, we can investigate a model selection case that was first encountered during a data analysis exercise Nofrarias et al. (2012). This type of exercises are organized by the DA team with the aim of testing the developed tools in more realistic scenarios. In this case, the “true” values of the parameters of the system to identify, were totally unknown. In that situation, the data analysis team was able to perform parameter estimation, but it was noted that the fit was improved when two extra parameters were introduced. These parameters are the guidance delays, and . They are simply time delays to the application of the guidance signals (see Fig. 1), due to operations of the Data Management Unit (DMU), which is the onboard computer controlling the LTP experiment. The final robustness of the fit indicates that the new parameters substantially improve the model.
This problem can be better studied by reducing it to a model selection problem and can be tested with synthetic datasets. The simulated data source is a LTP system with the default characteristics given in Table 1 with the exception of the delays, which were set to: . With this configuration, we assume a true system where the application signals are applied instantaneously. For our first tryout we injected “fake” interferometric displacements, while for the second investigation we used the same structure of injections but with much lower Signal to Noise Ratio (SNR ).
In order to determine the importance of the two extra parameters, we have to verify which model describes better the particular dataset: The seven parameter model X with parameters , or the five parameter model Y with parameters . While both models, X and Y, are capable of explaining the observations, we expect the simpler model Y to be more favorable from a RJMCMC output, since the extra parameters and are not significant for the data. The evolution of the Bayes factor for such an investigation for two different levels of SNR is shown in Fig. 2.
Method  

RJMCMC  0.309 
LF  0.124 
LM  0.078 
SBIC  0.768 
The results we obtain for all the approximations verify that the simpler model Y is much more probable than the more complicated model X. For the particular experiments proposed in Vitale (2007) where the SNR is high, the Bayes factor computed tends to zero. In fact, for the case of the RJMCMC method, there is no single iteration “inside” the more complex model. This changes dramatically depending on the nature of the problem and, of course, as we show in section IV.1, on the SNR. In Table 2 we present only the lowSNR experiment, for the sake of comparison between the methods. Each method seems to favor the simpler model but they are not in total agreement between them. This is to be expected, as the SNR of this investigation is very low and the approximations of the evidence become more sensitive. For this particular case of the injections, the models are not competitive and therefore, the resulting estimated Bayes Factor is extremely small. For a direct comparison, the RJMCMC algorithm requires more than  iterations.
iii.3 More realistic applications
During the last years a series of data analysis exercises have been performed. These so called Operational Exercises are organized by the LPF Data Analysis team with the aim of training the scientists and engineers for the upcoming flight operations. Data generated by the ESA simulator are used and reallife shifts and duties are rehearsed. The ESA simulator for the LPF is an industry developed software that is considered the most accurate simulator to date. The scientific and organizational results are valuable since they reproduce situations and data analysis challenges that are expected to occur during flight operations.
An outcome of these data challenges was the realization of the necessity to change or manipulate the model of the dynamics of the system shortly after receiving telemetry from the satellite, thus facing once again model selection problems. Expressing the dynamics in the acceleration domain Congedo et al. (2013), one could write a simplified model for the differential acceleration, like the one below:
(18) 
where is the differential acceleration and and denote the applied forces on the first and second testmasses respectively. The real motion of the testmasses can be approximated by the delayed interferometer readouts:
(19) 
The parameters appearing in Eq. (18) and (19) are the stiffnesses of the two testmasses and , the interferometer readout delay , and the gains of the capacitance actuators (here represented as for simplicity). Note that for this first approximation we have assumed identical actuators for both testmasses. The commanded forces and are available as telemetry, but the real applied forces on the three bodies are to be determined by the measurements and the analysis itself. For example, in the real datastream there might be additional delays or even filtering of the applied forces coming from the controllers, so in reality a gain might be proven to be frequency dependent: . This situation will appear in the received telemetry and we need the means to disentangle those two physical effects in a quantitative way.
For the particular simulated dataset, the model was not able to remove all the injected signals and even for the simple case of Eq. (18), this led to a biased estimate for the parameters. Apart from a simple time delay on the commanded signals, there might be another process that causes a difference between the commanded and and the actual applied forces and . For a first approximation of such a process, we assume a single real pole filter, filtering the timeseries of the applied force on the second testmass. For this investigation we can propose two models where and respectively. Here denotes the actuators time delay. In the end we can apply the model selection methods to these two models: the first one, X, where the applied forces are time delayed, and the second one, Y, where the forces are frequency dependent (filtered by a single real pole filter). The calculated Bayes Factor between those two models is:
(20) 
clearly indicating that the most probable process on the forces is the one described by model Y. The real pole was estimated to be Hz and was confirmed for the complete set of system identification experiments. The estimated acceleration residuals for both models can be seen in Fig. 4.
The same analysis strategy applies of course to more complicated versions of the analytical models, where the number of parameter increases with the more terms of the equation are added, as happens in the following more realistic case. As already described in section III, the two test masses are controlled with identical actuators. In reality, a small misbalance between the electrostatic actuators surrounding each test mass might be present. If we introduce this “asymmetry” to the system, immediately for the simple case of Eq. (18), we can increase the dimensionality of the model by four parameters:
(21) 
These parameters are gains, delays and filters that are different for the actuator of each test mass. Theoretically the highest in dimensions model of Eq.(21) can describe the observations, but the problem to solve appears to be overparametrized. A solution is to generate a set of nested models under the highest in dimensions of Eq. (21) and apply the RJMCMC algorithm. The result of such a run is shown in Fig. 4 and it reveals the most favorable model and consequently the underlying procedure that describe best the physical system. For the particular simulation we can verify that the fivedimensional model is the best, concluding that no “asymmetry” in the hardware of the LTP is present. This changed in the following Data Challenge, where the Bayes factor between the simple model of Eq. (18) and a sixparameter asymmetric one, is greater than , clearly supporting the correct higher dimension model where .
Iv Using the Bayes factor for experiment optimization
In this last section we explore the capability of the implemented model selection framework, not only as a set of tools for data analysis purposes, but also as a way to evaluate the efficiency of the experiments we are planning to run in the satellite. As we are going to show, by comparing different models under different input signal conditions, we can safely determine the best range of parameters that define our experiments, or verify the injection frequencies that maximize the information extracted from the system.
iv.1 The Bayes Factor as function of the SNR
It has been shown Cornish and Littenberg (2007); Veitch and Vecchio (2008) that there is a dependence of the Bayes Factor output on the SNR regime of the investigations. This, of course, holds true in the case of the LTP as it can be seen in Fig. 6. This figure was created by simulating LTP experiments for each value of the SNR, while the injection signals were single frequency () sinusoidal inputs into the system. The LTP models under comparison were quite similar with the exception of a different realization of the response model of the thrusters. It is clear that above the critical value of the the results obtained with the different techniques are consistent and in good agreement. Below that value of the SNR we cannot make clear decisions about the competing models, as the wrong model is showing preference, or we poorly approximate the Bayes Factor.
Although this SNR limit varies, as expected, depending on the type of investigation and model, the current result is already providing an estimation of the required power of the injection signals that we need to consider in the LTP experiment. This information and the method used here will be of particular interest for the design of inflight experiments for the LPF mission.
iv.2 The Bayes Factor as function of the injection frequencies
Furthermore, for system identification experiments, as in the case of the LTP, the computed evidence of a model depends on the design of the experiment itself. The information obtained from the system differs depending on the injection frequencies. An interesting study is to explore in detail this relation. A four (X) and a fivedimensional (Y) models are examined, given different injection frequencies. More precisely, since the difference between the models is the crosscoupling as shown in Fig. 1, which describes the signal leakage from the first to the differential interferometric channel, we examine the Bayes Factor given different injection frequencies to the first channel, while keeping constant the injection to the differential channel ( Hz). The SNR of this experiment is kept at the “low” value of 28.
The data generation model is mounted with a “perfect” interferometer () and model X is the same as the one used to produce the data, while model Y is the one with the extra parameter . The expected outcome of this exploration is that if the system is more sensitive to the parameter at some particular frequencies, we must detect an increase in the Bayes Factor which underlines a more clear decision towards the correct model.
In Fig. 6 we can see the corresponding Bayes Factor versus the injected frequencies to the first channel. Given the low SNR of the investigation, while model X should be more favorable, a preference for the more complex model Y is shown for a certain set of frequencies. This result is a clear indicator of the set of preferred frequencies that can be injected to the system for its characterization given the current configuration and SNR. Indeed, injections around Hz and Hz promote the identification of the correct model, while injections at both the high and low frequency limit, together with Hz, may induce the analysis into an error. This frequency dependence must be associated to the sensitivity of the experiment to a given parameter, the parameter in this particular case. This observed dependency, when considering a more realistic model, will be of particular interest in the selection of injection signals for the experiments to be run inflight.
V Conclusions
We have implemented three different methods to compare competing models of the LTP experiment onboard the LPF mission: The RJMCMC algorithm, the Laplace approximations, and the Bayes Information Criterion. The results from each method seem to be in agreement, but the output strongly depends on the expected SNR and of course on the models under investigation. Considering the LPF mission planned experiments, the SNR is high enough to safely use any of all the available techniques, but probably the most computationally demanding methods will be used for offline analysis to confirm our first computations.
The RJMCMC algorithm (together with the Laplace methods and the Bayes Information Criterion) employed in this work has been integrated in the LTPDA toolbox as part of the LPF data analysis software. The RJMCMC algorithm is by far the most computationally costly, but at the same time it is the more suitable one when we compare more than two nested models or we work with inputs with low SNR. The LaplaceMetropolis and the LaplaceFisher methods are reliable when we work in the high SNR regime, but they also require significant computing time, specially when one has to use outlier detection methods to estimate the weighted covariance matrix. On the other hand, the LaplaceFisher approximation is limited by the use of the Fisher Information Matrix, which for the case of LTP state space models is computed numerically.
Moreover, an attempt to associate the output of the aforementioned methods with the actual system identification experiment has been made. We have used different experiment setups to demonstrate that the Bayes Factor depends not only on the SNR, but also on the injection frequencies to the system.
The developed algorithms were successfully applied to model selection problems for the LPF data analysis for the first time.
Two different cases of LTP model selection problems have been investigated over datasets that were produced by both the
LTPDA and the ESA simulator. For the first case, we have considered an easy case of five and sevendimensional statespace models, where the importance of the
extra two parameters was examined. These two extra parameters are time delays caused by the LPF hardware and
they can be characterized as essential parameters of the model. For the second case, we explored the most suitable dimensionality of analytic models. There, the simplest model that described efficiently the observations
was recovered, excluding the more complicated ones that caused overfitting issues. This type of analysis is expected to be performed
during operations due to the broad spectrum of possible applications, like identifying external disturbances that result into forces applied to the threebody
system.
Acknowledgements.
We acknowledge support from contracts AYA201015709 (MICINN) and 2009SGR935 (AGAUR). NK acknowledges support from the FI PhD program of the AGAUR (Generalitat de Catalunya). MN acknowledges a JAEdoc grant from CSIC and support from the EU Marie Curie CIG 322288. GC acknowledges support from the Beecroft Institute for Particle Astrophysics and Cosmology.References
 Antonucci et al. (2012a) F. Antonucci, M. Armano, H. Audley, G. Auger, M. Benedetti, et al., Class.Quant.Grav. 29, 124014 (2012a).
 McNamara et al. (2013) P. McNamara, F. Antonucci, M. Armano, H. Audley, G. Auger, M. Benedetti, P. Binetruy, J. Bogenstahl, D. Bortoluzzi, N. Brandt, et al., in Astronomical Society of the Pacific Conference Series, edited by G. Auger, P. Binétruy, and E. Plagnol (2013), vol. 467 of Astronomical Society of the Pacific Conference Series, p. 5.
 Seoane et al. (2013) P. A. Seoane et al. (eLISA Collaboration) (2013), eprint 1305.5720.
 Nofrarias et al. (2010) M. Nofrarias, C. Röver, M. Hewitson, a. Monsky, G. Heinzel, K. Danzmann, L. Ferraioli, M. Hueller, and S. Vitale, Physical Review D 82, 1 (2010), ISSN 15507998, URL http://link.aps.org/doi/10.1103/PhysRevD.82.122002.
 Congedo et al. (2012) G. Congedo, L. Ferraioli, M. Hueller, F. De Marchi, S. Vitale, et al., Phys.Rev. D85, 122004 (2012), eprint 1108.0862.
 Antonucci et al. (2011) F. Antonucci, M. Armano, H. Audley, G. Auger, M. Benedetti, P. Binetruy, C. Boatella, J. Bogenstahl, D. Bortoluzzi, P. Bosetti, et al., Classical and Quantum Gravity 28, 094006 (2011), ISSN 02649381, URL http://stacks.iop.org/02649381/28/i=9/a=094006?key=crossref.41fc94b94ce025596175f5fe463e5a86.
 Cervantes et al. (2013) F. G. Cervantes, R. Flatscher, D. Gerardi, J. Burkhardt, R. Gerndt, M. Nofrarias, J. Reiche, G. Heinzel, K. Danzmann, L. G. Boté, et al., in Astronomical Society of the Pacific Conference Series, edited by G. Auger, P. Binétruy, and E. Plagnol (2013), vol. 467 of Astronomical Society of the Pacific Conference Series, p. 141.
 Audley et al. (2011) H. Audley, K. Danzmann, A. Garcia Marin, G. Heinzel, A. Monsky, et al., Class.Quant.Grav. 28, 094003 (2011).
 Kass and Raferty (1995) R. Kass and A. Raferty, Journal of American Statistical Assosiation 90 (1995).
 MacKay (2003) D. J. C. MacKay, Information Theory, Inference, and Learning Algorithms (Cambridge University Press, 2003).
 Hewitson et al. (2009) M. Hewitson, M. Armano, M. Benedetti, J. Bogenstahl, D. Bortoluzzi, et al., Class.Quant.Grav. 26, 094003 (2009).
 web (2005) LTPDA Matlab toolbox (2005), URL http://www.lisa.aeihannover.de/ltpda/.
 Green (1995) P. J. Green, Biometrica pp. 711–32 (1995).
 Lopes and West (2004) H. F. Lopes and M. West, Statistica Sinica 14, 41 (2004).
 Dellaportas et al. (2002) P. Dellaportas, J. J. Forster, and I. Ntzoufras, Statistics and Computing pp. 27–36 (2002).
 Stroeer and Veitch (2009) A. Stroeer and J. Veitch, Phys. Rev. D 80, 064032 (2009), URL http://link.aps.org/doi/10.1103/PhysRevD.80.064032.
 Umstätter et al. (2005) R. Umstätter, N. Christensen, M. Hendry, R. Meyer, V. Simha, J. Veitch, S. Vigeland, and G. Woan, Phys. Rev. D 72, 022001 (2005), URL http://link.aps.org/doi/10.1103/PhysRevD.72.022001.
 Godsill (2001) S. J. Godsill, JOURNAL OF COMPUTATIONAL AND GRAPHICAL STATISTICS 10, 230 (2001).
 Littenberg and Cornish (2009) T. B. Littenberg and N. J. Cornish, Phys. Rev. D 80, 063007 (2009), URL http://link.aps.org/doi/10.1103/PhysRevD.80.063007.
 Bartolucci et al. (2006) F. Bartolucci, L. Scaccia, and A. Mira, Biometrika 93, 41 (2006), ISSN 00063444.
 Cornish and Littenberg (2007) N. Cornish and T. Littenberg, Physical Review D 76, 1 (2007), ISSN 15507998, URL http://link.aps.org/doi/10.1103/PhysRevD.76.083006.
 Vallisneri (2008) M. Vallisneri, Physical Review D 77, 1 (2008), ISSN 15507998, URL http://link.aps.org/doi/10.1103/PhysRevD.77.042001.
 Rousseeuw and Leroy Annick (1987) P. J. Rousseeuw and M. Leroy Annick, Robust Regression and outlier detection (john Wiley and Sons, 1987).
 Barack and Cutler (2004) L. Barack and C. Cutler, Phys.Rev. D69, 082005 (2004), eprint grqc/0310125.
 Nofrarias et al. (2013) M. Nofrarias, F. Antonucci, M. Armano, H. Audley, G. Auger, M. Benedetti, P. Binetruy, J. Bogenstahl, D. Bortoluzzi, N. Brandt, et al., in Astronomical Society of the Pacific Conference Series, edited by G. Auger, P. Binétruy, and E. Plagnol (2013), vol. 467 of Astronomical Society of the Pacific Conference Series, p. 161.
 Hewitson et al. (2011) M. Hewitson, M. DiazAguilo, and A. Grynagier, Tech. Rep. S2AEITN3069 (2011).
 Grenadier and Weyrich (2008) A. Grenadier and M. Weyrich, Ph.D. thesis (2008).
 Congedo et al. (2013) G. Congedo, R. Dolesi, M. Hueller, S. Vitale, and W. J. Weber, Phys. Rev. D 88, 082003 (2013), URL http://link.aps.org/doi/10.1103/PhysRevD.88.082003.
 Vitale (2007) S. Vitale, Tech. Rep. S2UTNTN3045 (2007).
 Fichter et al. (2005) W. Fichter, P. Gath, S. Vitale, and D. Bortoluzzi, Classical and Quantum Gravity 22, S139 (2005), ISSN 02649381, URL http://stacks.iop.org/02649381/22/i=10/a=002?key=crossref.87240f1d344f89a51d67135cf12a4c33.
 Dolesi et al. (2003) R. Dolesi, D. Bortoluzzi, P. Bosetti, L. Carbone, A. Cavalleri, et al., Class.Quant.Grav. 20, S99 (2003).
 Heinzel et al. (2005) G. Heinzel, C. Braxmaier, M. Caldwell, K. Danzmann, F. Draaisma, et al., Class.Quant.Grav. 22, S149 (2005).
 Nofrarias et al. (2012) M. Nofrarias, L. Ferraioli, G. Congedo, M. Hueller, M. Armano, et al., J.Phys.Conf.Ser. 363, 012053 (2012), eprint 1111.4916.
 Antonucci et al. (2012b) F. Antonucci, M. Armano, H. Audley, G. Auger, M. Benedetti, P. Binetruy, J. Bogenstahl, D. Bortoluzzi, P. Bosetti, N. Brandt, et al., Classical and Quantum Gravity 29, 124014 (2012b), URL http://stacks.iop.org/02649381/29/i=12/a=124014.
 Veitch and Vecchio (2008) J. Veitch and A. Vecchio, pp. 1–10 (2008), eprint arXiv:0801.4313v3.