Christian Röver

Renate Meyer

Nelson Christensen

Modelling coloured residual noise in gravitational-wave signal processing

Keywords: numbers: 04.80.Nn, 02.30.Zz, 05.45.Tp

We introduce a signal processing model for signals in non-white noise, where the exact noise spectrum is a priori unknown. The model is based on a Student's t distribution and constitutes a natural generalization of the widely used normal (Gaussian) model. This way, it allows for uncertainty in the noise spectrum, or more generally is also able to accommodate outliers (heavy-tailed noise) in the data. Examples are given pertaining to data from gravitational wave detectors.

Introduction

The measurement of gravitational radiation holds great promise for exciting astronomical observations [START_REF] Thorne | Gravitational radiation[END_REF][START_REF] Schutz | Gravitational wave astronomy[END_REF]. Around the world, efforts are under way to construct and improve detectors for gravitational waves; among these are the LIGO detectors in the US [START_REF] Sigg | Status of the LIGO detectors[END_REF], GEO [START_REF] Grothe | The GEO 600 status[END_REF] and Virgo [START_REF] Accadia | Commissioning status of the Virgo interferometer[END_REF] in Europe, and TAMA in Japan [START_REF] Arai | Status of Japanese gravitational wave detectors[END_REF], some of which have already started taking data. Plans for future detectors include advanced LIGO [START_REF] Harry | Advanced LIGO: the next generation of gravitational wave detectors[END_REF], LCGT in Japan [START_REF] Kuroda | Status of LCGT[END_REF], the Laser Interferometer Space Antenna (LISA) [START_REF] Danzmann | LISA technology-concept, status, prospects[END_REF] and the Einstein Telescope (ET) [START_REF] Punturo | The third generation of gravitational wave observatories and their science reach[END_REF]. As existing instruments are becoming increasingly sensitive and future instruments are approaching completion, sophisticated signal processing methods are required in order to detect and accurately interpret these often weak signals within a noisy environment. For example, signal detection is usually implemented via matched filtering [START_REF] Finn | Detection, measurement, and gravitational radiation[END_REF][START_REF] Turin | An introduction to matched filters[END_REF][START_REF] Wainstein | Extraction of signals from noise[END_REF], while parameter estimation is commonly done using Bayesian methods [START_REF] Röver | Coherent Bayesian inference on compact binary inspirals using a network of interferometric gravitational wave detectors[END_REF][START_REF] Searle | Bayesian detection of unmodeled bursts of gravitational waves[END_REF]. Many of the data analysis procedures employed to date are based on the assumption of the noise being Gaussian with a known power spectral density [START_REF] Finn | Detection, measurement, and gravitational radiation[END_REF][START_REF] Jaranowski | Gravitational-wave data analysis. Formalism and sample applications: The Gaussian case[END_REF]. While these methods have proven to be computationally efficient and very powerful at discriminating rare, weak signals within noise, some more flexibility or robustness is sometimes desired. One example is the analysis of data to be expected from LISA, where the measurement noise originates partly from instrumental as well as astrophysical sources, and the noise spectrum is only vaguely known a priori [START_REF] Barack | Confusion noise from LISA capture sources[END_REF].

We introduce an approach that was developed in the context of the latter case, where, along with the signal parameters to be inferred, the noise's power spectrum needed to be incorporated as an unknown into the model [START_REF] Röver | Bayesian inference on astrophysical binary inspirals based on gravitationalwave measurements[END_REF][START_REF] Röver | Inference on inspiral signals using LISA MLDC data[END_REF]. The model developed here turns out to be computationally convenient, as the additional noise parameters Confidential: not for distribution. Submitted to IOP Publishing for peer review 9 November 2010 may be analytically integrated out, leading to a Student-t likelihood expression instead of the original normal formulation. We expect the same approach to be useful in many related signal processing contexts, as it allows to specify the prior information about the power spectrum, i.e., its expected magnitude and the associated certainty, in a straightforward manner. In this way, it should e.g. also be able to account for nonstationarities in the data, in the sense that the power spectrum does not necessarily need to be assumed to be exactly the same as at an earlier measurement; other fluctuations, like outliers, could also be accommodated in a similar manner. In fact, the resulting model actually falls into a class of generalizations to the standard Gaussian that have been proposed for their robustness properties: models accommodating heavier-tailed noise, or outliers, by either ('directly') implementing a non-Gaussian noise model or ('indirectly') substituting the corresponding least-squares procedures by less outlier-sensitive methods. Instead of the derivation based on an unknown power spectrum, the same model may be motivated by assuming the power spectrum itself to be random, so that the resulting noise mixture distribution exhibits a greater variability. In the limiting case of decreasing spectrum variability, the model again simplifies to the Gaussian. So the model will also be applicable as an ad-hoc alternative tunable robust model, with clearly interpretable "tuning parameters".

The organization of the paper is as follows. The time series setup is introduced in Section 2.1. In the subsequent sections, the probabilistic modelling is described, including its time-domain counterpart, the likelihood, prior distributions, posterior distribution, marginal likelihood and some implications. Section 3 describes the approach using two illustrative examples of simulated time series with and without a signal. The discussion section 4 puts the new approach into context. An appendix explicating the Discrete Fourier Transform conventions used in this paper is attached.

The time series model

The setup

Consider a time series x 1 , . . . , x N of N real-valued observations sampled at constant time intervals ∆ t , so that each observation x i corresponds to time t i = i∆ t . This set of N observations can equivalently be expressed in terms of sinusoids of the Fourier frequencies:

x i = 1 √ N ∆ t N/2 j=0 a j cos(2πf j t i ) + b j sin(2πf j t i ) ( 1 )
where the variables a j and b j each correspond to Fourier frequencies f j = j∆ f = j N ∆t . The summation in (1) runs from j = 0 to j = N/2 , with N/2 denoting the largest integer less than or equal to N/2. By definition, b 0 is always zero, and b N/2 is zero if N is even; going over from x i 's to a j 's and b j 's again yields the same number (N ) of non-zero figures. The set of (N ) frequency domain coefficients a j and b j and the time domain observations x i are related to each other through a discrete Fourier transform (and appropriate scaling; see appendix Appendix A.2). The set of trigonometric functions in (1) constitutes an orthonormal basis of the sample space, so that there is a unique one-to-one mapping of the observations in time and in frequency domain. Instead of the two amplitudes a j and b j , the definition in (1) may equivalently be expressed in terms of a single amplitude and a phase parameter at each frequency:

x i = 1 √ N ∆ t N/2 j=0 λ j sin(2πf j t i + ϕ j ), (2)
where λ j = a 2 j + b 2 j and

ϕ j =    arctan bj aj if a j > 0 arctan bj aj ± π if a j < 0. (3)
For each Fourier frequency f j , let κ j be the number of Fourier coefficients not being zero by definition, i.e.:

κ j = 1 if (j = 0) or (N is even and j = N/2) 2 otherwise (4)
so that N/2 j=0 κ j = N . Note that generally one may often simplify to κ j = 2 ∀j without introducing a noticeable discrepancy, but in the following we will stick to the accurate notation of κ depending on the index j.

For a given time series (either in terms of x i or a j and b j ), we define the functions of the Fourier frequencies

p 1 (f j ) = a 2 j + b 2 j κ j and (5)
p 2 (f j ) = a 2 j + b 2 j κ 2 j = p 1 (f j ) κ j = ∆t N |x(f j )| 2 (6)
for j = 0, . . . , N/2 , where x denotes the discretely Fourier-transformed time series x, as defined in appendix Appendix A.2. These are the empirical, discrete analogues of the one-sided and two-sided spectral power. The set of p 1 (f j ) is also known as the periodogram of the time series [START_REF] Jenkins | General considerations in the analysis of spectra[END_REF]. Now consider the case where the observations x i , and consequently the a j and b j , correspond to random variables X i , A j and B j , respectively. This may mean that these are realizations of a random process, where the random variables' probability distributions describe the randomness in the observations, or that they are merely unknown, where the probability distributions describe a state of information, or it may also be a mélange of both [START_REF] Jaynes | Probability theory: The logic of science[END_REF].

The time series has a zero mean if and only if the expectation of all frequency domain coefficients is zero as well:

E[X i ] = 0 ∀i ⇔ E[A j ] = E[B j ] = 0 ∀j. (7)
For the probabilistic time series, the spectral power (p 1 (f j ) or p 2 (f j )) consequently also is a random variable (P 1 (f j ) or P 2 (f j ), respectively). We denote its expectation value by S * (f j ); if the mean is zero as in [START_REF] Harry | Advanced LIGO: the next generation of gravitational wave detectors[END_REF], it is given by

S * 1 (f j ) = κ j S * 2 (f j ) = E[P 1 (f j )] = E A 2 j +B 2 j κj (7)

=

Var(Aj )+Var(Bj ) κj

.

It is important to note that while the above expectation is closely related to a time series' power spectral density, it is yet quite different. Up to here we have considered the discretely Fourier transformed data and some of its basic properties. The figure S * (f j ) refers to data sampled at a particular resolution and sample size and is defined only at a discrete set of Fourier frequencies f j . One might want to refer to S * as the discretized power spectrum. If in fact the data are a realization from a stationary random process with (continuous) power spectral density S(f ), then the expectation S * (f j ) is related to S(f ) through a convolution, depending on the sample size N . Only in the limit of an infinite sample size both of them are equal [START_REF] Jenkins | General considerations in the analysis of spectra[END_REF][START_REF] Gregory | Bayesian logical data analysis for the physical sciences[END_REF]. While both continuous and discretized spectrum are commonly used as an approximation to or in place of one another, it is crucial to realize that in the following inference will be done with reference to the case of a finite sample size and the discretized, convolved spectrum S * (f j ).

A related useful figure in this context is the integrated spectrum with respect to some frequency range [f 1 , f 2 ] [START_REF] Jenkins | Spectral analysis and its applications[END_REF], which here we define as

I [f1,f2] = ∆ f j2 j=j1 κ j 2 S * 1 (f j ) ( 9 )
where the summation is done over the corresponding range of Fourier frequencies (j 1 = min{k :

f k > f 1 }, j 2 = max{k : f k ≤ f 2 }, where f 2 -f 1 ≥ ∆ f ).
The integrated spectrum allows to investigate or compare (discrete) spectra independent of the particular sample size N .

The normal model

Suppose the moments as defined in [START_REF] Harry | Advanced LIGO: the next generation of gravitational wave detectors[END_REF] and ( 8) are given. We may set up a corresponding normal time series model by assuming all the (frequency domain) observables to be stochastically independent, have zero mean and

Var(A j ) = σ 2 j , (10)
Var(B j ) = (κ j -1) σ 2 j for j = 0, . . . , N/2 ,

where

σ 2 j = S * 1 (f j ) = κ j S * 2 (f j ) for j = 0, . . . , N/2 . ( 12
)
While other choices of variance settings matching the assumptions would be possible, the assumption of equal variances for A j and B j is the only one that makes the joint density of (A j , B j ) a function of overall amplitude λ j alone, independent of the phase ϕ j , and with that leaves the model invariant with respect to time shifts. The same model may also be motivated via the maximum entropy principle, as the (zero mean, equal variance, independent) normal distribution maximizes the entropy given the constraints on the moments given in ( 7) and ( 8) above [START_REF] Bretthorst | The near-irrelevance of sampling frequency distributions[END_REF]. In that way, the normal model constitutes a most convervative model setup under the given assumptions [START_REF] Jaynes | Probability theory: The logic of science[END_REF][START_REF] Gregory | Bayesian logical data analysis for the physical sciences[END_REF][START_REF] Bretthorst | The near-irrelevance of sampling frequency distributions[END_REF]. The joint normal distribution derived here has actually been commonly used before (see e.g. [START_REF] Fisher | Tests of significance in harmonic analysis[END_REF]); the normality assumption may not only appear as a natural choice, but will also turn out computationally convenient in the following. In the context of Fourier domain data in particular, the normality assumption may also be motivated via asymptotic arguments, by considering the limit of an infinite observation time [START_REF] Brillinger | Time series: Data analysis and theory[END_REF][START_REF] Champeney | A handbook of Fourier theorems[END_REF][START_REF] Kawata | On the Fourier series of a stationary stochastic process[END_REF]. The resulting normal model also is exactly the same as the one underlying the so-called Whittle likelihood [START_REF] Whittle | Curve and periodogram smoothing[END_REF][START_REF] Choudhuri | Contiguity of the Whittle measure for a Gaussian time series[END_REF], or the one at the basis of matched filtering [START_REF] Finn | Detection, measurement, and gravitational radiation[END_REF][START_REF] Turin | An introduction to matched filters[END_REF][START_REF] Wainstein | Extraction of signals from noise[END_REF]. While intuitively in other contexts it may often be sufficient to point out the asymptotic equality S * (f j ) ≈ S(f j ) for N → ∞, here it is crucial to appreciate what exactly the σ 2 j = S * 1 (f j ) stand for in the case of a finite sample size N .

The difference between the exact and approximate ("Whittle") model is explicated in more detail in [START_REF] Choudhuri | Contiguity of the Whittle measure for a Gaussian time series[END_REF].

The corresponding time-domain model

An immediate consequence of the normality assumption for the frequency-domain coefficients (A j , B j ) is that the time-domain variables (X i ), being linear combinations of the frequency-domain variables, also follow a normal distribution. The exact joint distribution of the X i is completely determined by their variance/covariance structure, which may be expressed in terms of the autocovariance function. The covariance for any pair of time-domain observations X m and X n (with m, n ∈ {1, . . . , N}, and corresponding to times t m and t n ) is given by:

Cov(X m , X n ) = 1 N ∆ t N/2 j=0 S * 1 (f j ) κj 2 cos 2πf j (t n -t m ) . ( 13
)
Note that the autocovariance Cov(X m , X n ) = γ(t n -t m ) depends on t m and t n only via their time lag t n -t m . Remarkably (though not surprisingly), following either the invariance or the maximum entropy argument to motivate the model structure in section 2.2, the result is a strictly stationary model for the data [START_REF] Taqqu | Weak stationarity[END_REF]. Note also that

γ(0) = Var(X i ) = I [0,f N/2 ] .
The autocovariance function allows one to express the distribution of the X i in terms of their (N × N ) covariance matrix Σ X , which is the symmetric, square Toeplitz matrix

Σ X =        γ(0) γ(∆ t ) γ(2∆ t ) • • • γ(∆ t ) γ(∆ t ) γ(0) γ(∆ t ) • • • γ(2∆ t ) γ(2∆ t ) γ(∆ t ) γ(0) • • • γ(3∆ t ) . . . . . . . . . . . . . . . γ(∆ t ) γ(2∆ t ) γ(3∆ t ) • • • γ(0)        , ( 14
)
since γ is periodic such that γ(i∆ t ) = γ((N -i)∆ t ). The periodic formulation in (1) makes the first and last observations X 1 and X N "neighbours", just as X 1 and X 2 are. This may seem odd (depending on the context, of course), but is not an unusual problem, as it arises for any conventional spectral analysis of discretely sampled data in the form of spectral leakage [START_REF] Harris | On the use of windows for harmonic analysis with the discrete Fourier transform[END_REF]; it may just not always be as obvious in its consequences. The above time domain expression sheds some light on the exact shortcomings of the Whittle likelihood approximation; you can see for example that it will be a poor approximation in case of predominantly low-frequency noise and a small number of observations. The problem may be tackled via windowing of the data, or it will also lessen with an increasing sample size; again, the features of this approximation are discussed in more detail in [START_REF] Choudhuri | Contiguity of the Whittle measure for a Gaussian time series[END_REF].

Incorporating an unknown spectrum 2.4.1. The likelihood function

Now suppose the spectrum is unknown and to be inferred from the data x 1 , . . . , x N . An unknown spectrum here is equivalent to the variance parameters σ 2 0 , . . . , σ 2 N/2 being unknown. Assuming normality of A j and B j , the likelihood function (as a function of the parameters σ 2 0 , . . . , σ 2 N/2 ) is:

p(x 1 , . . . , x N | σ 2 0 , . . . , σ 2 N/2 ) = N/2 j=0 1 √ 2π σ j exp - a 2 j 2σ 2 j × 1 √ 2π σ j κj -1 exp - b 2 j 2σ 2 j (15) = exp -N 2 log(2π) - N/2 j=0 κ j log(σ j ) + a 2 j + b 2 j 2σ 2 j ( 16
) ∝ exp - N/2 j=0 κ j 2 log S * 2 (f j ) + ∆t N |x j | 2 S * 2 (f j ) . ( 17
)
The term xj here denotes the jth element of the Fourier-transformed data vector, as defined in appendix Appendix A.1 and Appendix A.2. The above likelihood ( 17) is exactly equivalent to the so-called Whittle likelihood, which is an approximate expression for stationary time series [START_REF] Whittle | Curve and periodogram smoothing[END_REF][START_REF] Choudhuri | Contiguity of the Whittle measure for a Gaussian time series[END_REF].

If the noise spectrum is a priori known, the "log S * 2 (f j ) " term may also be irrelevant (see [START_REF] Taqqu | Weak stationarity[END_REF]), and the model again reduces to the normal model described e.g. in [START_REF] Finn | Detection, measurement, and gravitational radiation[END_REF] that is also the basis for matched filtering [START_REF] Turin | An introduction to matched filters[END_REF][START_REF] Wainstein | Extraction of signals from noise[END_REF]. In that case, the logarithmic likelihood may be conveniently expressed in terms of an inner product of vectors, allowing for a geometric interpretation that is expecially useful in the case of linear signal models [START_REF] Finn | Detection, measurement, and gravitational radiation[END_REF][START_REF] Jaranowski | Gravitational-wave data analysis. Formalism and sample applications: The Gaussian case[END_REF][START_REF] Mood | Introduction to the theory of statistics[END_REF].

The conjugate prior distribution

For the model defined above, the conjugate prior distribution for each of the σ 2 0 , . . . , σ 2 N/2 is the scaled inverse χ 2 -distribution with scale parameter s 2 j and degrees-of-freedom parameter ν j :

σ 2 j ∼ Inv-χ 2 (ν j , s 2 j ) ( 1 8 )
with density function

f νj ,s 2 j σ 2 j = νj s 2 j 2 νj /2 Γ(ν j /2) σ 2 j -(1+νj /2) exp - ν j s 2 j 2 σ 2 j ( 19
) [START_REF] Gelman | Bayesian data analysis[END_REF]. The degrees-of-freedom here denote the precision in the prior distribution, while the scale determines its order of magnitude. For increasing ν j the distribution's variance goes towards zero, and for ν j → 0 the density converges toward the noninformative (and improper) distribution with density f (σ

2 j ) = 1 σ 2 j
, that is uniform on log(σ 2 j ), and which also constitutes the corresponding Jeffreys prior for this problem [START_REF] Jeffreys | An invariant form for the prior probability in estimation problems[END_REF]. If σ 2 j follows an Inv-χ 2 (ν j , s 2 j ) distribution, then its expectation and variance are given by

E[σ 2 j ] = νj νj -2 s 2 j and Var(σ 2 j ) = 2ν 2 j (νj -2) 2 (νj -4) s 4 j , (20)
and the mean and variance are finite for ν j > 2 and ν j > 4, respectively [START_REF] Gelman | Bayesian data analysis[END_REF].

In addition to the Jeffreys prior (with ν = 0) already mentioned above, other improper prior distributions may be implemented as special cases of an Inv-χ 2 (ν, s 2 ) distribution. A uniform prior distribution on σ corresponds to (ν = -1, s 2 = 0), and a uniform prior on σ 2 corresponds to (ν = -2, s 2 = 0). Generally, a prior with density f (σ 2 ) ∝ 1 (σ 2 ) k (for k ≥ 0) corresponds to an Inv-χ 2 (ν = 2(k -1), s 2 = 0) distribution [START_REF] Bolstad | Introduction to Bayesian Statistics[END_REF]. As usual, care must be taken when using such improper priors, as the resulting posterior may not be a proper probability distribution [START_REF] Gelman | Bayesian data analysis[END_REF].

The posterior distribution

Due to the use of a conjugate prior distribution, the posterior distribution of the σ 2 j for given data again is of the same family. The posterior density is defined by the product of prior [START_REF] Röver | Inference on inspiral signals using LISA MLDC data[END_REF] and likelihood [START_REF] Jaranowski | Gravitational-wave data analysis. Formalism and sample applications: The Gaussian case[END_REF]:

p σ 2 j {x 1 , . . . , x N } ∝ σ 2 j -(1+νj /2) exp - ν j s 2 j 2 σ 2 j × σ 2 j -κj /2 exp - a 2 j + b 2 j 2 σ 2 j (21) = σ 2 j -(1+ ν j +κ j 2 ) exp - ν j s 2 j + a 2 j + b 2 j 2 σ 2 j ( 22
)
which again can be recognized as a scaled inverse χ 2 -density (cp. ( 19)):

σ 2 j {x 1 , . . . , x N } ∼ Inv-χ 2 ν j + κ j , νj s 2 j + a 2 j +b 2 j νj + κj , (23)
where all the different parameters corresponding to different frequencies are mutually independent. Comparing prior and posterior parameters (( 18), ( 23)), and the way the prior and likelihood are combined [START_REF] Jaynes | Probability theory: The logic of science[END_REF], one can see that the prior distribution might be thought of as providing the information equivalent to ν j observations (of coefficients a j or b j ) with average squared deviation s 2 j [START_REF] Gelman | Bayesian data analysis[END_REF]. The prior is essentially of the same functional form as the likelihood [START_REF] Jaynes | Probability theory: The logic of science[END_REF], modulo a 1 σ 2 term, which resembles an overall (uninformative) Jeffreys prior prefactor [START_REF] Jeffreys | An invariant form for the prior probability in estimation problems[END_REF]. The use of the conjugate prior distribution hence is not only a computationally convenient choice, but may also appear as a "natural" way of expressing prior information in this context.

The use of the conjugate prior distribution leads to a convenient expression for the posterior distribution of the variance parameters σ 2 j , and with that of the complete discrete spectrum. Also, if σ 2 j ∼ Inv-χ 2 (ν j , s 2 j ), then, since s 2 j is a scale parameter, it follows that the distribution of the two-sided spectrum S * 2 (f j ) = σ 2 j /κ j simply is Inv-χ 2 (ν j , s 2 j /κ j ). There is a deterministic relationship between given variance parameters σ 2 j and the implied autocorrelation function γ(t) (see ( 13)), and so for random σ 2 j the distribution of the corresponding γ(t) may be numerically explored via Monte Carlo sampling from the distribution of the σ 2 j . The expectation and variance of γ(t) may also be derived analytically: the expected autocovariance [START_REF] Wainstein | Extraction of signals from noise[END_REF], with respect to the distribution of the σ 2 j , is

E[γ(t)] = 1 N ∆ t N/2 j=0 E[σ 2 j ] κj 2 cos(2πf j t) , (24)
which is finite as long as E[σ 2 j ] is finite for all j. Similarly, the variance of the autocorrelation is

Var(γ(t)) = 1 N 2 ∆ 2 t N/2 j=0 Var(σ 2 j ) κ 2 j 4 cos(2πf j t) 2 , ( 25
)
which again is finite as long as Var(σ 2 j ) is finite for all j.

The marginal likelihood

The use of the conjugate Inv-χ 2 prior (with ν j > 0 degrees of freedom) for the variance parameters σ 2 j allows to integrate the unknown noise spectrum out of the likelihood expression [START_REF] Barack | Confusion noise from LISA capture sources[END_REF]. The marginal likelihood function then is

p(x 1 , . . . , x N | ν 0 , . . . , ν N/2 , s 2 0 , . . . , s 2 N/2 ) = N/2 j=0 ∞ 0 p(a j , b j |σ 2 j ) p(σ 2 j |ν j , s 2 j ) dσ 2 j (26) = N/2 j=0 (2 π) - κ j 2 νj s 2 j 2 ν j 2 νj s 2 j +(a 2 j +b 2 j ) 2 ν j +κ j 2 Γ νj +κj 2 Γ νj 2 (27)
∝ N/2 j=0 1 + κ 2 j ∆t N |x j | 2 ν j s 2 j - ν j +κ j 2 (28)
= exp -

N/2 j=0 νj +κj 2 log 1 + κ 2 j ∆t N |x j | 2 ν j s 2 j , (29)
which constitutes a product of Student-t densities with (ν j +κ j -1) degrees of freedom [START_REF] Gelman | Bayesian data analysis[END_REF]. Using the uninformative, improper Jeffreys prior (with ν j = 0 degrees of freedom and density p(σ 2 j ) ∝ 1

σ 2 j
) yields the marginal likelihood

p(x 1 , . . . , x N ) ∝ N/2 j=0 κ 2 j ∆t N |x j | 2 -κ j 2 . ( 30
)
One may then also get a mixture of the above expressions (( 28) and ( 30)) in case of a prior setting of partly zero and non-zero prior degrees of freedom. Note that the corresponding analogue likelihood expression in case of an a priori known spectrum (as e.g. in [START_REF] Whittle | Curve and periodogram smoothing[END_REF][START_REF] Finn | Detection, measurement, and gravitational radiation[END_REF][START_REF] Choudhuri | Contiguity of the Whittle measure for a Gaussian time series[END_REF], see also [START_REF] Barack | Confusion noise from LISA capture sources[END_REF]) was

p(x 1 , . . . , x N |σ 2 0 , . . . , σ 2 N/2 ) ∝ exp - N/2 j=0 κ 2 j ∆t N |x j | 2 2 σ 2 j , (31)
so that the (normal) sum-of-squares expression [START_REF] Taqqu | Weak stationarity[END_REF] in case of a known spectrum generalizes to the Student-t expression (29) once one takes uncertainty about the spectrum's scale into consideration. The normal model in turn is the limiting case for increasing degrees of freedom.

Robust inference via the Student-t model

The Student-t marginal likelihood expression [START_REF] Kawata | On the Fourier series of a stationary stochastic process[END_REF] suggests that considering noise as normal but with unknown spectrum is technically equivalent to viewing the noise itself as being t-distributed. In that way, the Student-t model may also be useful as an alternative for robust modelling, as the wider family of t-distributions includes the normal and Cauchy distributions as special or limiting cases. In other contexts the t-distribution is commonly used as a generalization of the normal model in order to accommodate heavier-tailed errors (see e.g. [START_REF] Lange | Robust statistical modeling using the t distribution[END_REF][START_REF] Geweke | Bayesian treatment of the independent Student-t linear model[END_REF][START_REF] Divgi | Robust estimation using Student's t distribution[END_REF][START_REF] Mcdonald | Partially adaptive estimation of regression models via the generalized t distribution[END_REF]). In contrast to the above derivation of the t-distributed noise based on prior/posterior distributions, the noise may also be considered as a scale mixture of normal distributions; that is, normal noise with a randomly varying scale (variance) [START_REF] Geweke | Bayesian treatment of the independent Student-t linear model[END_REF]. Both relations may actually be used to motivate the use of the Student-t distribution for noise exhibiting certain nonstationarities (e.g., a noise spectrum that is slightly fluctuating over time), be it to model the "variability" or the "uncertainty" in the spectrum, both of which are mathematically the same here (except that uncertainty is reduced through learning, while the random variation is not). The use of heavier-tailed noise models already has repeatedly been advocated in the gravitational wave detection context, in order to gain robustness against outliers and similar deviations from the commonly assumed normality. Creighton [START_REF] Creighton | Data analysis strategies for the detection of gravitational waves in non-Gaussian noise[END_REF] for example suggests the use of a mixture of Gaussian noise and a (uniform) burst component in order to account for wide-band noise burst events. Allen et al. [START_REF] Allen | Robust statistics for deterministic and stochastic gravitational waves in non-Gaussian noise: Frequentist analyses[END_REF] also illustrate the use of a two-component Gaussian mixture in order to reduce the influence of outliers in the data, but more generally they advocate an approach also known as M-estimation, namely the explicit downweighting or ignorance of outlier observations falling far into the tail of the noise distribution [START_REF] Hampel | Robust statistics: The approach based on influence functions[END_REF][START_REF] Huber | Robust statistics[END_REF]. Application of the above Student-t model may also be considered as a special case of robust Mestimation, based on a clearly interpretable noise model and associated parameters [START_REF] Divgi | Robust estimation using Student's t distribution[END_REF][START_REF] Mcdonald | Partially adaptive estimation of regression models via the generalized t distribution[END_REF].

In case of an a priori known power spectrum, maximization of the normal likelihood [START_REF] Taqqu | Weak stationarity[END_REF] is the basis for the matched filtering approach commonly used in signal detection problems, when looking for signals of parameterized shape in noise [START_REF] Turin | An introduction to matched filters[END_REF][START_REF] Wainstein | Extraction of signals from noise[END_REF][START_REF] Jaranowski | Gravitational-wave data analysis. Formalism and sample applications: The Gaussian case[END_REF]. In Gaussian noise, the likelihood maximization then is equivalent to a least-squares approach. A filter based on the Student-t likelihood (28) may therefore be useful for cases of an unknown noise spectrum or non-Gaussian noise.

Note that since the likelihood (28) implies independent (as opposed to merely uncorrelated) errors, likelihood maximization here is different from least-squares estimation [START_REF] Kelejian | Independent or uncorrelated disturbances in linear regression: An illustration of the difference[END_REF][START_REF] Breusch | The emperor's new clothes: a critique of the multivariate t regression model[END_REF][START_REF] Geweke | Bayesian treatment of the independent Student-t linear model[END_REF], and the likelihood might in fact exhibit multiple modes [START_REF] Mäkeläinen | On the existence and uniqueness of the maximum likelihood estimate of a vector-valued parameter in fixed-size samples[END_REF]. The (independent) Student-t model not only leads to a less drastic fall-off of the likelihood for extreme values, and hence reduced leverage of outliers, but it also implies non-spherical density contours for the joint distribution of the noise, effectively allowing for a fraction of scrambled individual noise residuals. A similar effect was pointed out by Creighton [START_REF] Creighton | Data analysis strategies for the detection of gravitational waves in non-Gaussian noise[END_REF] when implementing a robust Gaussian/uniform mixture noise model: the approach would allow for excess noise in individual interferometers, so that a noise burst would essentially be automatically "vetoed" if it is only measured in one of several interferometers. Similarly, the non-spherical density contours of the Student-t model make it robust against odd data values at individual frequencies.

Defining the prior distribution's parameters

Depending on the particular application and context in mind, there may be different ways to sensibly specifying prior parameters. Firstly, there are the supposedly uninformative priors; the Jeffreys prior [START_REF] Jeffreys | An invariant form for the prior probability in estimation problems[END_REF] with ν j = 0 degrees of freedom, and priors that are uniform on σ j or on σ 2 j were already mentioned above (see Sec. 2.4.2). Care needs to be taken here though, as these priors (with ν j ≤ 0) are improper distributions. These may lead to improper marginal likelihoods, as in the case of [START_REF] Choudhuri | Contiguity of the Whittle measure for a Gaussian time series[END_REF], which does not correspond to a normalizable probability distribution for the noise. The resulting posterior distribution of signal or noise parameters then may or may not be a proper probability distribution [START_REF] Gelman | Bayesian data analysis[END_REF].

If the prior information on the spectral parameters is in the form of measurements (or samples) of the spectrum (of same size and resolution), then this may be expressed in terms of an equivalent sample size and corresponding prior degrees-of-freedom, as suggested in Sec. 2.4.3. Choosing e.g. the ν j to be twice the number of measurements of the spectrum taken (i.e., equal to the number of observed Fourier coefficients) would then be equivalent to initially assuming an uninformative Jeffreys prior and then using the posterior based on the measurements as the prior for the actual (signal) analysis. This will of course only make sense if one assumes the spectrum unknown but stationary.

When using the Student-t model as a robust model accommodating for heavytailed noise, i.e., when the noise itself is assumed to actually follow a t-distribution that one can sample from, then one can use sample estimates for the Student-t parameters. Moment estimators for s 2 j and ν j (based on sample variance and kurtosis) are given e.g. in [START_REF] Sutradhar | Estimation of the parameters of a regression model with a multivariate t error variable[END_REF][START_REF] Singh | Estimation of error variance in linear regression models with errors having multivariate Student-t distribution with unknown degrees of freedom[END_REF].

When specifying prior parameters that are supposed to reflect information and/or variability, it may be helpful to consider the implied moments or quantiles, for individual frequency bins or frequency bands. The expressions for the prior's moments in [START_REF] Jenkins | General considerations in the analysis of spectra[END_REF] may be inverted to

ν j = 4 + 2 E[σ 2 j ] 2 Var(σ 2 j )
and

s 2 j = ν j -2 ν j E[σ 2 j ], (32)
which allows one to specify the scale s 2 j and degrees-of-freedom ν j based on pre-defined prior expectation and variance of σ 2 j , respectively. Note that the degrees-of-freedom ν j then are simply a function of the prior variation coefficient Var(σ 2 j ) E[σ 2 j ]. A specification of s 2 j independent of j means a priori white noise, and specifying individual ν j for different j indicates varying prior certainty across the spectrum. A sensible definition of the prior certainties for the individual spectrum parameters may be complicated by the fact that the exact meaning of this discrete set of parameters depends on the sample size N . In that case it may be helpful to instead consider the integrated spectrum (9) and its a priori properties. The (prior) moments of the integrated spectrum are given by

E[I [f1,f2] ] = ∆ f j2 j=j1 κj 2 E[σ 2 j ] = ∆ f j2 j=j1 κ j 2 ν j ν j -2 s 2 j ( 33
)
(if all ν j > 2), and

Var(I [f1,f2] ) = ∆ 2 f j2 j=j1 κ 2 j 4 2ν 2 j (ν j -2) 2 (ν j -4) s 4 j ( 34
)
(if all ν j > 4). The (prior) variation coefficient for the power within any frequency range then is:

Var(I [f1,f2] ) E[I [f1,f2] ] = j2 j=j1 κ 2 j 4 2ν 2 j (νj -2) 2 (νj -4) s 4 j j2 j=j1 κj 2 νj νj -2 s 2 j ( 35
)
which simplifies in case all d.f. parameters are taken to be equal (ν j ≡ ν):

Var(I [f1,f2] ) E[I [f1,f2] ] = 2 ν -4 j2 j=j1 κ 2 j 4 s 4 j j2 j=j1 κj 2 s 2 j (36)
and simplifies further in case all scale parameters are taken to be equal (s 2 j ≡ s 2 ):

Var(I [f1,f2] ) E[I [f1,f2] ] = 2 ν -4 j2 j=j1 κ 2 j 4 j2 j=j1 κj 2 . ( 37
)
In the general case of (j 1 > 0, j 2 < N/2 ) this is:

Var(I [f1,f2] ) E[I [f1,f2] ] = 2 ν -4 1 √ j 2 -j 1 + 1 ( 38
)
and similarly, in case (j 1 = 0, and j 2 = N/2, N even) it is:

Var(I [f0,f N/2 ] ) E[I [f0,f N/2 ] ] = 2 ν -4 N -1 2 N 2 . ( 39
)
This may be useful if one wants to specify piecewise constant prior settings with given constraints on the overall power per frequency range; this would then lead to the d.f. settings

ν = 4 + 2 j2-j1+1 E[I [f 1 ,f 2 ] ] √ Var(I [f 1 ,f 2 ] ) 2 (40) or ν = 4 + 2 N -1 N 2 E[I [f 0 ,f N/2 ] ] Var(I [f 0 ,f N/2 ] ) 2 (41)
respectively, for the above two cases (( 38), (39)).

For example, if one wants the spectrum to be a priori white with some scale s 2 , such that the marginal prior mean and variation coefficient of the integrated power

I [f0,f N/2 ] = Var(X i ) = γ(0) are given by E I [f0,f N/2 ] = ς 2 and
Var

I [f0,f N/2 ] E I [f0,f N/2 ] = c, (42)
then a setting of

ν j = 4 + N -1 N 2 2 c 2 and s 2 j = 2 ∆ t ν-2 ν ς 2 (43)
independent of j will yield an a priori white spectrum with constant expectation and variation in the overall variance I [f0,f N/2 ] for any sample size N . Piecewise constant settings across different frequency bands may be implemented analogously. In case all ν j > 4, the integrated power, being a sum of random variables with finite mean and variance, will be asymptotically normally distributed. Instead of using moments, the prior may also be specified in terms of quantiles, for example by first deciding on the degrees-of-freedom settings and then aiming the prior median at a certain value. The quantiles of an Inv-χ 2 (ν, s 2 ) distribution may be derived based on the quantiles of a χ 2 -distribution; the p-quantile is then given by ν s 2 χ 2 ν;1-p [START_REF] Huber | Robust statistics[END_REF] where χ 2 ν;1-p is the (1 -p)-quantile of a χ 2 -distribution with ν degrees of freedom. Finally, in the M-estimation context, the t-distribution's d.f. parameter may also be set based on the shape of the corresponding influence function [START_REF] Hampel | Robust statistics: The approach based on influence functions[END_REF][START_REF] Huber | Robust statistics[END_REF].

Examples

Noise only

Consider a time series n(t) of N = 100 data points sampled at times t i = i 100 . As a simple example of non-white noise, the data are generated from an autoregressive process n(t i ) = 3 4 n(t i-1 ) + x(t i ), [START_REF] Kelejian | Independent or uncorrelated disturbances in linear regression: An illustration of the difference[END_REF] where the innovations x(t i ) are drawn independently from a uniform distribution across the interval [-

√ 3, √ 3
], so that they have zero mean and unit variance. The overall variance of the process defined in this way is Var n(t i ) = 2.29; due to the positive correlation of subsequent samples it has a higher power at low frequencies and less at high frequencies. Figure 1 shows a noise sample generated using the above prescription, together with its theoretical power spectral density.

We will now apply the noise model introduced above and derive the posterior distribution of the noise parameters. Assuming one has a rough idea of the noise variance, we set the prior scale parameters s 2 j so that the noise is a priori white with a prior expectation of E Var n(t i ) = 2.50. We use ν j = 3 degrees of freedom, so that the noise parameters' (and with that, the overall power's) prior expectations are finite, while the variances are not. Prior scale and prior expectation then are s 2 j = 0.0166, and E[σ 2 j ] = 0.05. Figure 2 illustrates the resulting posterior distribution for all 51 noise parameters [START_REF] Jenkins | Spectral analysis and its applications[END_REF] in comparison to the case of using the uninformative (and improper) Jeffreys prior. The Jeffreys prior (with ν j = 0 degrees of freedom for each frequency bin) does not depend on the scale parameters s 2 j , and the resulting posterior with ≤ 2 degrees of freedom at each frequency does not have finite expectation values. Note also that the posterior distributions corresponding to the first and last frequency bin (zero and Nyquist frequencies, σ 2 0 and σ 2 50 ) are wider than the others in both cases, as they have one less degree of freedom.

Figure 3 shows the posterior distributions of autocovariance and variance, which are functions of the individual spectrum parameters σ 2 j . The variance may either be considered as the zero-lag autocovariance γ(0) (see [START_REF] Wainstein | Extraction of signals from noise[END_REF]), or as the integrated power I [0,f N/2 ] across the whole frequency range (see [START_REF] Danzmann | LISA technology-concept, status, prospects[END_REF]). Posterior distributions of autocovariance and variance; the distributions shown here were derived via Monte Carlo integration.

A signal with additive noise

In the following example we will consider a time series where the primary interest is in a signal component, while the additive noise component will be modelled using the approach introduced above. As in the previous section, we will again consider N = 100 data points y(t 1 ), . . . , y(t 100 ) that are modelled as

y(t i ) = g f, ḟ,a,φ (t i ) + n σ (t i ), (46)
where n σ (t i ) is non-white noise of unknown spectrum, and g f, ḟ,a,φ (t) is a "chirping" signal waveform of increasing frequency:

g f, ḟ ,a,φ (t) = a sin(2π(f + ḟt)t + φ) ( 4 7 )
where f and ḟ are the frequency and frequency derivative, a is the amplitude, and φ is the phase. The noise again is generated the same way as in the previous example, The three densities in each panel result from three different models applied for the noise term. The "unknown, coloured" case corresponds to the model introduced above, the "fixed, coloured" case assumes the true noise spectrum to be a priori known, and in the "unknown white" case the noise is modelled as white with unknown variance. The vertical dashed lines indicate the true parameter values.

by simulating an autoregressive process. The signal's "size" relative to the noise is given by the signal-to-noise ratio (SNR: = 4 j

∆ t N |g(fj )| 2
S1(fj ) ), which here is at 15. We define the signal parameters' prior as uniform for phase, frequency and amplitude (φ ∈ [0, 2π], f ∈ [START_REF] Thorne | Gravitational radiation[END_REF][START_REF] Arnaud | An overview of the Mock LISA Data Challenges[END_REF] and a ∈ [0, 10]) and normal for the frequency derivative ḟ (zero mean and standard deviation 5). The prior distribution for the noise parameters σ 2 0 , . . . , σ 2 50 is set exactly as in the previous section 3.2. For comparison, we analyze the data two more times, once assuming the true noise spectrum to be known, and once assuming the noise to be white, but with an unknown variance. In the latter case, we again assume a (conjugate) Inv-χ 2 prior with an expectation of 2.5 and 3 degrees of freedom for the variance parameter.

The posterior distributions of signal and noise parameters may now be derived via Monte Carlo integration; we implemented a Metropolis sampler to simulate draws from the joint posterior probability distribution of all parameters [START_REF] Gelman | Bayesian data analysis[END_REF]. For the model including the coloured noise spectrum as unknown, the signal parameters (f , ḟ , a, φ) may be sampled based on the marginal likelihood expression [START_REF] Whittle | Curve and periodogram smoothing[END_REF], while samples of the 51 noise parameters (σ 2 0 , . . . , σ 2 50 ) may then be sampled in an additional step via the conditional distribution of noise parameters for given signal parameters and the corresponding implied vector of noise residuals [START_REF] Jenkins | Spectral analysis and its applications[END_REF]. Similarly, for the known spectrum model, sampling of the four signal parameters may be based on the likelihood expression [START_REF] Taqqu | Weak stationarity[END_REF], while sampling for the white noise model may be done using a Gibbs sampler alternately sampling from the conditional distributions of four signal parameters and the noise parameter [START_REF] Gelman | Bayesian data analysis[END_REF]. Figure 4 shows the resulting marginal posterior distributions of the four signal parameters in comparison. One can see that application of the more flexible coloured noise model for the error term yields a more precise posterior distribution than the white noise model, and in fact the posterior is more similar to the case when the true noise spectrum is plugged into the model. Looking at the (marginal) posterior distributions of the noise parameters in figure 5, one can see that although the noise spectrum is only recovered with great uncertainty, this does not seem to harm the estimation of the signal parameters of actual interest. The overall variance is better recovered by the single-parameter white noise model, whereas the adaptability of the coloured noise model to the predominantly low-frequency noise (which is reflected in the posterior) seems to be the greater advantage. While the relative accuracy of the different methods is subject to a multitude of circumstances, one can see that depite its considerably greater complexity, the coloured noise model seems to perform competitively here.

Discussion

This work originated out of the Mock LISA Data Challenges (MLDC) [START_REF] Arnaud | An overview of the Mock LISA Data Challenges[END_REF], a gravitational wave parameter estimation effort in the context of the planned Laser Interferometer Space Antenna (LISA). Here parameters of a signal were to be inferred, where the signal was buried in noise known to be non-white and interspersed with a host of individual emission lines [START_REF] Barack | Confusion noise from LISA capture sources[END_REF]. So the problem was to model a non-white, non-continuous spectrum where the shape of the spectrum was only vaguely known in advance [START_REF] Röver | Bayesian inference on astrophysical binary inspirals based on gravitationalwave measurements[END_REF][START_REF] Röver | Inference on inspiral signals using LISA MLDC data[END_REF]. Also, the spectrum itself was not of primary interest, but rather a nuisance parameter that still needed to be accounted for along with the actual signal. Application of the method described here solved the problem and allowed the implementation of a Markov chain Monte Carlo (MCMC) algorithm based on a straightforward generalization of the commonly utilized Gaussian noise assumption, where modelling the spectrum complicated the analysis only slightly.

Several approaches to the simultaneous estimation of signal and noise parameters have been proposed before. An obvious way to model non-white noise would be to assume a particular time series formulation that allows for some flexibility in the resulting spectrum, for example an autoregressive (AR) model. Representations of this kind have been applied in various signal processing contexts, e.g. for detecting sinusoidal signals and estimating their parameters when these are buried in coloured noise [START_REF] Chatterjee | Estimation of close sinusoids in colored noise and model discrimination[END_REF][START_REF] Cho | Bayesian detection and estimation of cisoids in colored noise[END_REF], or in the context of musical pitch estimation, where the signals to be modelled again are sinusoidal, but include higher harmonics and time-varying amplitudes [START_REF] Godsill | Bayesian harmonic models for musical pitch estimation and analysis[END_REF]. If one is only interested in a narrow frequency band of the data, it may also be appropriate to model the noise spectrum as constant across the concerned range [START_REF] Dupuis | Bayesian estimation of pulsar parameters from gravitational wave data[END_REF]. Regarding data treated in their Fourier domain representation, there is a rich literature concerned with the measurement of spectral densities (see e.g. [START_REF] Brillinger | Time series: Data analysis and theory[END_REF] and references therein). However, there the aim is usually to produce consistent and smooth estimates of the spectrum, and the approaches applied include e.g. averaging [START_REF] Welch | The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms[END_REF][START_REF] Harris | On the use of windows for harmonic analysis with the discrete Fourier transform[END_REF], smoothing via splines [START_REF] Wahba | Automatic smoothing of the log periodogram[END_REF] or Bayesian model fitting [START_REF] Gangopadhyay | Estimation of spectral density of a stationary time series via an asymptotic representation of the periodogram[END_REF][START_REF] Choudhuri | Bayesian estimation of the spectral density of a time series[END_REF]. The approach to modelling the noise spectrum introduced here is different in that the spectrum per se is not of interest, or only of interest as far as it enters into likelihood computations. Smoothness or interpolation therefore are not primarily aimed for. What in fact is of concern is properly accounting for an a priori uncertainty in the (discretized, convolved) spectrum, in the frequency resolution that is given by the numerically Fourier transformed data, which then consequently entails the necessity for updating our knowledge about the spectrum as data is being processed. The resulting approach generalizes the model underlying the commonly used Whittle likelihood [START_REF] Whittle | Curve and periodogram smoothing[END_REF][START_REF] Finn | Detection, measurement, and gravitational radiation[END_REF][START_REF] Choudhuri | Contiguity of the Whittle measure for a Gaussian time series[END_REF], and is hence essentially based on the "plain" periodogram of the noise time series, which for other purposes is commonly dismissed due to its unfavourable large-sample convergence behaviour. The model is very flexible as it is built upon a "binned" spectrum estimate without introducing any extra assumptions on the shape of the underlying noise spectral density. Its generality and simplicity make it useful for modelling residual noise at very little computational cost. In the way it is defined, the model is very general; it is a generalization of the model that constitutes the basis for matched filtering (e.g. [START_REF] Turin | An introduction to matched filters[END_REF][START_REF] Wainstein | Extraction of signals from noise[END_REF]) and that is commonly applied in signal processing problems (e.g. [START_REF] Finn | Detection, measurement, and gravitational radiation[END_REF]), which then in turn constitutes the special case of an a priori known spectrum. The additional feature of the approach introduced here is that it allows to specify corresponding uncertainties in addition to the (prior) scale of the noise spectrum. Marginalization over the uncertainty in the noise spectrum yields a Student-t model for the Fourier-domain data as a natural generalization of the common normal model. Due to the straightforward interpretability of the model and its computationally convenient form, we expect it to be particularly useful for modelling an unknown noise spectrum that constitutes a nuisance parameter rather than being of interest in itself. In that sense, the approach is not primarily aimed at gaining information about an unknown spectrum, but rather at properly accounting for uncertainty, and avoiding bias from supposedly precise a priori knowledge of the spectrum.

Alternatively, the Student-t model may also be viewed as a generalized, robust model accommodating for heavier-tailed, non-Gaussian, or non-stationary noise. Related approaches are commonly used in many other applications in the context of robust statistics [START_REF] Lange | Robust statistical modeling using the t distribution[END_REF][START_REF] Geweke | Bayesian treatment of the independent Student-t linear model[END_REF][START_REF] Divgi | Robust estimation using Student's t distribution[END_REF][START_REF] Mcdonald | Partially adaptive estimation of regression models via the generalized t distribution[END_REF][START_REF] Hampel | Robust statistics: The approach based on influence functions[END_REF][START_REF] Huber | Robust statistics[END_REF]. In fact, the use of similar methods for accommodating noise outliers have already been advocated in the context of gravitational wave signal processing [START_REF] Creighton | Data analysis strategies for the detection of gravitational waves in non-Gaussian noise[END_REF][START_REF] Allen | Robust statistics for deterministic and stochastic gravitational waves in non-Gaussian noise: Frequentist analyses[END_REF]. Along similar lines, Clark et al. [START_REF] Clark | Evidence-based search method for gravitational waves from neutron star ring-downs[END_REF] implemented a (sine-Gaussian) noise-glitch component into the noise model. Principe and Pinto [START_REF] Principe | Locally optimum network detection of unmodelled gravitational wave bursts in an impulsive noise background[END_REF] suggested the use of a dictionary of glitch "atoms" in order to account for burst-like noise events. Veitch and Vecchio [START_REF] Veitch | Bayesian coherent analysis of in-spiral gravitational wave signals with a detector network[END_REF] implemented a coherence test based on a model selection procedure that is able to discriminate noise glitches (that appear independently in the data) from actual astrophysical events (that appear coherently in several data streams). Similarly, Littenberg and Cornish [START_REF] Littenberg | Separating gravitational wave signals from instrument artifacts[END_REF] extend their model to allow for excess noise that is isolated in time and frequency via a wavelet approach. Middleton [START_REF] Middleton | On the theory of random noise. Phenomenological models[END_REF] treated the general case of Gaussian noise with superimposed Poissondistributed impulse noise bursts. We are planning to investigate the performance and sensitivity of a Student-t model for robust detection and parameter estimation in real interferometer noise.

We expect the approach to be especially useful as a model component properly accounting for non-white residual noise, at little computational cost and without introducing overly restrictive assumptions about the noise. The interpretability of the Student-t model in the context of robust modeling and M-estimation also makes it straightforwardly useful for robust inference. In that way, it may particularly be useful in signal processing contexts [START_REF] Mcdonough | Detection of signals in noise[END_REF][START_REF] Wainstein | Extraction of signals from noise[END_REF], but it should be applicable in any case where a model for (residual) noise is needed [START_REF] Jenkins | Spectral analysis and its applications[END_REF][START_REF] Bloomfield | Fourier analysis of time series: An introduction[END_REF]. The fact that the model represents noise properties in terms of its power spectral density and is specified through physically meaningful parameters may make it particularly appealing in physical or engineering applications, where modelling is commonly based on Fourier-domain descriptions, and a time-domain formulation might be hard to incorporate or motivate. The advantage from application of the more flexible Student-t model will very much depend on the particular inference problem at hand. The common normal model will obviously be optimal when its assumptions are met, and the degree to which one will outperform the other will depend on the particular departure(s) from Gaussianity, or on the imperfect prior knowledge of the spectrum, even if the data are perfectly Gaussian. While the range of possible deviations from the standard Gaussian model is infinite, some insight into the behaviour of a Student-t model may be gained from the robust statistics literature; the properties of such an approach for location estimation has been investigated via Monte Carlo studies [START_REF] Andrews | Robust estimates of location: survey and advances[END_REF][START_REF] Divgi | Robust estimation using Student's t distribution[END_REF] as well as theoretical considerations of key figures like relative efficiency or breakdown point [START_REF] Divgi | Robust estimation using Student's t distribution[END_REF][START_REF] Hampel | Robust statistics: The approach based on influence functions[END_REF][START_REF] Huber | Robust statistics[END_REF]. A case study in the regression context may be found in [START_REF] Geweke | Bayesian treatment of the independent Student-t linear model[END_REF].

The above basic model may in future be extended by introducing smoothness constraints on the spectrum. This might for example be approached by considering correlations between neighbouring spectral bins, or by assuming a piecewise constant spectrum, similar to what was done in [START_REF] Dupuis | Bayesian estimation of pulsar parameters from gravitational wave data[END_REF]; the latter approach would in fact be a compromise between two models used in the example discussed above (Sec. 3.2), namely a flat spectrum and individually modelled frequency bins. Another interesting extension would be the incorporation of cross-spectra [START_REF] Jenkins | Spectral analysis and its applications[END_REF][START_REF] Nofrarias | Bayesian parameter estimation in the second LISA Pathfinder Mock Data Challenge[END_REF]. Some of the methods described here have been coded as an R software extension package and are available for free download [START_REF] Röver | bspec: Bayesian spectral inference[END_REF].

Figure 1 .

 1 Figure 1. The example data n i in the time domain, and its (empirical) power ∆t N |ñ j | 2 in the Fourier domain. The solid line in the right panel shows the theoretical (2-sided) power spectral density in comparison.

Figure 2 .

 2 Figure 2. Posterior distributions of the 51 spectrum parameters σ 2 j based on the data shown in figure 1. The left plot shows the posterior corresponding to the uninformative (and improper) Jeffreys prior. Posterior expectations do not exist in this case. The right plot corresponds to assuming a priori white noise with ν j = 3 degrees of freedom for each frequency bin; the dashed line marks the prior expectation value.

Figure 4 .

 4 Figure 4. Marginal posterior densities of the four individual signal parameters.The three densities in each panel result from three different models applied for the noise term. The "unknown, coloured" case corresponds to the model introduced above, the "fixed, coloured" case assumes the true noise spectrum to be a priori known, and in the "unknown white" case the noise is modelled as white with unknown variance. The vertical dashed lines indicate the true parameter values.

Figure 5 .

 5 Figure 5. Marginal posterior distributions of the noise parameters. The left panel shows the posterior distribution for the 51 individual parameters of the coloured noise model. The right panel shows the prior and posterior distributions of the integrated power in comparison with the corresponding variance parameter in the white noise model.

Acknowledgments

This work was supported by the Max-Planck-Society, The Royal Society of New Zealand Marsden Fund grant UOA-204, and National Science Foundation grants PHY-0553422 and PHY-0854790.

Appendix

Appendix A. [START_REF] Thorne | Gravitational radiation[END_REF]

. Discrete Fourier transform (DFT)

The Fourier transform convention used in this paper is specified below; it is defined for a real-valued function h of time t, sampled at N discrete time points, at a sampling rate of 1 ∆t , and it maps from

to a function of frequency f

where

The inverse DFT then is given by

Appendix A.2. Relationship between DFT and time series model

Let α j = Re h(f j ) and β j = Im h(f j ) , (A.5)

i.e.: h(f j ) = α j + β j i. For simplicity, in the following N is assumed to be even; for uneven N the derivation is similar. The inverse DFT was defined as (A.4):

2α j cos(2πf j t) + 2(-β j ) sin(2πf j t)

where t ∈ {0, ∆ t , 2∆ t , . . . , (N -1)∆ t }, and f j = j∆ f = j N ∆t are the Fourier frequencies. So, comparing to (1), one can see that the realizations of a 0 , . . . , a N/2 and b 0 , . . . , b N/2 are derived from a given time series by Fourier-transforming and then setting

for j = 0, . . . , N/2 , which especially implies that

(A.9)