Gaussian regression and power spectral density estimation with missing data: The MICROSCOPE space mission as a case study

Quentin Baghi, Gilles Métris, Joël Bergé, Bruno Christophe, Pierre Touboul, and Manuel Rodrigues

Phys. Rev. D 93, 122007 – Published 16 June 2016

Abstract

We present a Gaussian regression method for time series with missing data and stationary residuals of unknown power spectral density (PSD). The missing data are efficiently estimated by their conditional expectation as in universal Kriging based on the circulant approximation of the complete data covariance. After initialization with an autoregressive fit of the noise, a few iterations of estimation/reconstruction steps are performed until convergence of the regression and PSD estimates, in a way similar to the expectation-conditional-maximization algorithm. The estimation can be performed for an arbitrary PSD provided that it is sufficiently smooth. The algorithm is developed in the framework of the MICROSCOPE space mission whose goal is to test the weak equivalence principle (WEP) with a precision of $10^{- 15}$ . We show by numerical simulations that the developed method allows us to meet three major requirements: to maintain the targeted precision of the WEP test in spite of the loss of data, to calculate a reliable estimate of this precision and of the noise level, and finally to provide consistent and faithful reconstructed data to the scientific community.

Received 21 March 2016

DOI:https://doi.org/10.1103/PhysRevD.93.122007

Physics Subject Headings (PhySH)

Experimental studies of gravity Gravitational waves

Data analysis

Gravitation, Cosmology & Astrophysics

Authors & Affiliations

Quentin Baghi^1,*, Gilles Métris^2,†, Joël Bergé¹, Bruno Christophe¹, Pierre Touboul¹, and Manuel Rodrigues¹

¹ONERA—The French Aerospace Lab, 29 avenue de la Division Leclerc, 92320 Châtillon, France
²Université Nice Sophia Antipolis, CNRS, Observatoire de la Côte d’Azur, Géoazur UMR 7329, 250 rue Albert Einstein, Sophia Antipolis, 06560 Valbonne, France

^*quentin.baghi@onera.fr
^†gilles.metris@oca.eu

Article Text (Subscription Required)

Click to Expand

References (Subscription Required)

Click to Expand

Issue

Vol. 93, Iss. 12 — 15 June 2016

Reuse & Permissions

Access Options

Author publication services for translation and copyediting assistance advertisement

Images

Figure 1
Extract of the time series obtained for a 20 orbit spin session sampled at 4 Hz, for the periodic gaps (a) and the random gaps (b). For clarity, we plotted a sample which is longer for periodic gaps (lasting more than one orbital period $T_{orb}$ , or about six WEP periods $1 / f_{EP}$ ) than for random gaps (which are 260 times shorter and more scattered). Observed data are in black and reconstructed data from one Monte Carlo conditional draw are in blue. The missing data spans are indicated by gray areas.
Reuse & Permissions
Figure 2
Lomb-Scargle periodogram of the observed data (gray), periodogram of the reconstructed data (blue), and of the original data (black) for a 20 orbit spin session and for the periodic (a) and random (b) gap patterns. In order not to overload the figure, the values have been plotted on a subgrid of about 700 Fourier frequencies. The harmonic peaks of highest amplitudes visible around $2 f_{EP} = 1.87 \times 10^{- 3} Hz$ and common to the blue and black periodograms are not an artifact but are due to the gravitational gradient perturbation $T_{x x}$ and $T_{x z}$ (which are included in the model $A$ ). This is also the case for the peak at $f_{orb} = 1.70 \times 10^{- 4} Hz$ .
Reuse & Permissions
Figure 3
Average of the 400 Lomb-Scargle periodogram of the observed data (gray), of the periodogram of the reconstructed data (blue), and of the original data (black) for a 20 orbit spin session and for the periodic (a) and random (b) gap patterns. The same remarks made in Fig. 2 about the peaks of the gravitational gradient perturbation apply. Averaging reveals the peak at $f_{EP} = 9.35 \times 10^{- 4} Hz$ that is due to the simulated WEP violation. It also makes another faint gradient harmonic visible at $2 f_{orb} = 3.40 \times 10^{- 4} Hz$ , which is due to the influence of the second zonal harmonic $J_{2}$ .
Reuse & Permissions
Figure 4
Sample average of the PSD estimates (blue) along with its 99% confidence interval (light blue area), autoregressive estimate (dashed red), and true PSD (black) for a 20 orbit spin session and for the periodic (a) and random (b) gap patterns.
Reuse & Permissions
Figure 5
Sample average of the normalized autocovariance estimated from random-gapped data with the M-ECM algorithm (blue) along with its 99% confidence interval (dashed blue) and the true autocovariance (black). The autocovariance alternates between positive and negative values. The bar heights of the histogram correspond to the values of $R (τ) / R (0)$ at lags $τ$ located at the left-hand-side edge of the bars. To be comprehensive, the confidence interval should include both the upper and lower bounds for positive and negative values of the autocovariance. For the sake of clarity, only the upper bound of positive values and the lower bound of negative values are displayed.
Reuse & Permissions
Figure 6
Average of 400 conditional periodograms and PSD estimates obtained with the M-ECM algorithm applied to data samples generated from the toy model defined by Eq. (c2) with $N = 1000$ data points. The blue and red curves correspond to numbers of MC draws, respectively, equal to $N_{d} = 5$ and $N_{d} = 100$ . Light solid lines represent average periodograms while dashed lines represent PSDs. The gray curve is the average of the complete data periodogram, and the black dashed curve is the true PSD. Therefore, increasing $N_{d}$ lowers the bias of the PSD estimate.
Reuse & Permissions

Physical Review D

covering particles, fields, gravitation, and cosmology