Bayesian cross validation for gravitational-wave searches in pulsar-timing array data

TL;DR

This paper introduces Bayesian cross validation as a robust, prior-insensitive method for model selection in pulsar-timing-array gravitational-wave data analysis, demonstrating its effectiveness on simulated and real datasets.

Contribution

It presents Bayesian cross validation as an alternative to Bayesian model comparison, with advantages in simplicity, prior insensitivity, and overfitting detection in gravitational-wave searches.

Findings

Successfully identifies the correct spectral index in simulations

Favors models with transient features in real pulsar data

Offers a practical, prior-insensitive model selection approach

Abstract

Gravitational-wave data analysis demands sophisticated statistical noise models in a bid to extract highly obscured signals from data. In Bayesian model comparison, we choose among a landscape of models by comparing their marginal likelihoods. However, this computation is numerically fraught and can be sensitive to arbitrary choices in the specification of parameter priors. In Bayesian cross validation, we characterize the fit and predictive power of a model by computing the Bayesian posterior of its parameters in a training dataset, and then use that posterior to compute the averaged likelihood of a different testing dataset. The resulting cross-validation scores are straightforward to compute; they are insensitive to prior tuning; and they penalize unnecessarily complex models that overfit the training data at the expense of predictive performance. In this article, we discuss cross validation in the context of pulsar-timing-array data analysis, and we exemplify its application to simulated pulsar data (where it successfully selects the correct spectral index of a stochastic gravitational-wave background), and to a pulsar dataset from the NANOGrav 11-year release (where it convincingly favors a model that represents a transient feature in the interstellar medium). We argue that cross validation offers a promising alternative to Bayesian model comparison, and we discuss its use for gravitational-wave detection, by selecting or refuting models that include a gravitational-wave component.

Equations

$$p(y|M)=\int p(y|{\theta }_M)p({\theta }_M)\mathrm{d}{\theta }_M,$$

$$\log p(\tilde{y}|y;M)=\log \int p(\tilde{y}|{\theta }_M)p({\theta }_M|y)\mathrm{d}{\theta }_M.$$

$$\begin{gathered} \log \int p(y^{(k)}|{\theta }_M)p({\theta }_M|y^{(-k)})\mathrm{d}\theta \\ \simeq \log \frac{1}{N}\sum _{i=1}^{N}p(y^{(k)}|{\theta }_{M,i}^{(-k)}); \end{gathered}$$

$$\begin{gathered} p(y|{\eta }_N,{\eta }_{\mathrm{GP}},c_{\mathrm{GP}})=p(y|{\eta }_N,c_{\mathrm{GP}})\times p(c_{\mathrm{GP}}|{\eta }_{\mathrm{GP}}) \\ =\frac{{\mathrm{e}}^{-(y-F_{\mathrm{GP}}{c}_{\mathrm{GP}}{)}^T{N}^{-1}(y-F_{\mathrm{GP}}{c}_{\mathrm{GP}})/2}}{\sqrt{(2\pi {)}^n|N|}}\times \frac{{\mathrm{e}}^{-c_{\mathrm{GP}}^T{\Phi }_{\mathrm{GP}}^{-1}{c}_{\mathrm{GP}}/2}}{\sqrt{(2\pi {)}^m|\Phi |}}, \end{gathered}$$

$$\begin{gathered} p(y|{\eta }_N,{\eta }_{\mathrm{GP}})=\int p(y|{\eta }_N,c_{\mathrm{GP}},{\eta }_{\mathrm{GP}})\mathrm{d}{c}_{\mathrm{GP}} \\ =\frac{{\mathrm{e}}^{-y^T(N+F_{\mathrm{GP}}{\Phi }_{\mathrm{GP}}{F}_{\mathrm{GP}}^T{)}^{-1}y/2}}{\sqrt{(2\pi {)}^n|N+F_{\mathrm{GP}}{\Phi }_{\mathrm{GP}}{F}_{\mathrm{GP}}^T|}}. \end{gathered}$$

$$\bar{c}(y;{\eta }_N,{\eta }_{\mathrm{GP}})=\Sigma F_{\mathrm{GP}}^T{N}^{-1}y$$

$$\Sigma ({\eta }_N,{\eta }_{\mathrm{GP}})=({\Phi }_{\mathrm{GP}}^{-1}+F_{\mathrm{GP}}^T{N}^{-1}{F}_{\mathrm{GP}}{)}^{-1}.$$

$$P(f)=A_{\mathrm{red}}^2{\left(\frac{f}{f_{\mathrm{yr}}}\right)}^{-{\gamma }_{\mathrm{red}}},$$

$$\begin{gathered} \log \int p(Y^{(k)}|{\eta }^{(k)},{\theta }_{\mathrm{GW}})p({\theta }_{\mathrm{GW}}|Y^{(k)})p({\eta }^{(k)})\mathrm{d}{\eta }^{(k)}\mathrm{d}{\theta }_{\mathrm{GW}} \\ \simeq \log \frac{1}{N}\sum _{i=1}^N\int p(Y^{(k)}|{\eta }^{(k)},{\theta }_{\mathrm{GW},\mathrm{i}}^{(-k)})p({\eta }^{(k)})\mathrm{d}{\eta }^{(k)}, \end{gathered}$$

$$\begin{gathered} \log \sum _{i=1}^N\int p(Y^{(2)}|{\eta }^{(2)},c_{\mathrm{GW}}^{(2)})p(c_{\mathrm{GW}}^{(2)}|c_{\mathrm{GW}}^{(1)},{\theta }_{\mathrm{GW},i})\times \\ p(c_{\mathrm{GW}}^{(1)}|Y^{(1)};{\eta }_i^{(1)},{\theta }_{\mathrm{GW},i})\mathrm{d}{\eta }^{(2)}\mathrm{d}{c}_{\mathrm{GW}}^{(2)}\mathrm{d}{c}_{\mathrm{GW}}^{(1)}. \end{gathered}$$

$$p(c_{\mathrm{GW}}^{(2)}|Y^{(1)};{\eta }_i^{(1)},{\theta }_{\mathrm{GW},i})=\mathcal{N}({\bar{c}}^{(2)|(1)},{\Sigma }^{(2)|(1)})$$

$${\bar{c}}^{(2)|(1)}(Y^{(1)};{\eta }_i^{(1)},{\theta }_{\mathrm{GW},i})={\Phi }_{21}{\Phi }_{11}^{-1}{\bar{c}}^{(1)},$$

$$\begin{gathered} {\Sigma }^{(2)|(1)}({\eta }_i^{(1)},{\theta }_{\mathrm{GW},i}) \\ ={\Phi }_{22}-{\Phi }_{21}({\Phi }_{11}^{-1}-{\Phi }_{11}^{-1}{\Sigma }^{(1)}{\Phi }_{11}^{-1}){\Phi }_{21}, \end{gathered}$$

Full text

1 Introduction

Searches for gravitational waves (GWs) in pulsar-timing-array (PTA) data (Lommen, 2015; Burke-Spolaor, 2015) seek to identify weak GW signals among a plethora of other effects, including deterministic delays due to the relative motion of pulsar and observatory and to pulsar binary dynamics, stochastic delays due to the interplanetary and interstellar media, as well as intrinsic irregularities in the pulsar’s period emission (Cordes, 2013; Stinebring, 2013). These searches are commonly formulated as Bayesian-inference problems (Gregory, 2010), whereby we derive the joint posterior probability density of the GW parameters and of the noise parameters of all analyzed pulsars. Choosing appropriate probabilistic models for pulsar noise is therefore crucial to reliable PTA searches (Taylor et al., 2013; Lentati et al., 2016; Cordes & Shannon, 2010): unmodeled noise components may be interpreted as GWs, while overgenerous noise assumptions may reduce GW sensitivity. In current practice, pulsar noise models are informed by the physics of millisecond pulsars and of the interplanetary/interstellar medium, but they are largely driven by inference from PTA datasets, since these often represent the best observations to date for PTA pulsars.

2 Model comparison

Within the data-analysis practice of the NANOGrav collaboration (NANOGrav, 2019), Bayesian model comparison (Gregory, 2010) is used to select the noise model most appropriate to each pulsar (Arzoumanian et al., 2019). The goal is not only to improve the physical characterization of the processes affecting pulse times of arrival (TOAs), but also to isolate these processes from a putative GW signal with greater confidence. In this framework, we evaluate the fully marginalized likelihood (a.k.a. evidence) for each model $M$:

[TABLE]

where $y$ denotes the observed data, $\theta_{M}$ the parameters of model $M$, $p(y|\theta_{M})$ the likelihood (the probability of $y$ given $\theta_{M}$), and $p(\theta_{M})$ the prior probability density assigned to the parameters. We then compare models by evaluating Bayes ratios***Bayesian model comparison calls for the computation of odds ratios, which account for the prior relative probability of entire models. Since it is very difficult to attribute such priors on physical grounds, we generally work directly with Bayes ratios. $B_{21}=p(y|M_{2})/p(y|M_{1})$, either directly through Eq. (1) (often requiring significant numerical sophistication, see Trotta 2008) or by the Monte Carlo exploration of uber-likelihoods that specialize to individual models depending on the value of an index parameter (in which case the Bayes ratio is given by the ratio of the “time” spent in each model, see Godsill 2001; Sisson 2005). A large Bayes ratio $B_{21}$ implies that the data favors model $M_{2}$ over $M_{1}$. However, it is difficult to give Bayes ratios a principled quantitative interpretation. The exception are cases where alternative models represent exclusive physical outcomes; Bayes ratios can then be calibrated in terms of statistical decision theory, by relating their sampling distribution to false-alarm and false-dismissal probabilities (see, e.g., Vallisneri 2012).

In choosing between alternative models for a dataset, we need to be wary of overfitting: that is, while it is always possible to improve model fit by adding parameters, the enhanced model may end up conforming to contingent noise features instead of highlighting the physical properties of interest. Correspondingly, the model loses predictive power for yet-to-observed data. Bayesian model comparison incorporates a defense against overfitting, in that the evidence integral penalizes fine tunings that restrict parameters to small regions within their prior ranges. Unfortunately, this defense creates a different, significant weakness—Bayes ratios are then sensitive to parameter-prior assignments that may be largely arbitrary, and are not testable from data (Gelman et al., 2013). Consider for instance the case of two nested models $M_{1}\subset M_{2}$, where $M_{2}$ is obtained by adding parameter $\theta^{*}$ to $M_{1}$; let $\theta^{*}$ have uniform prior in $[-a,a]$, with $a$ arbitrarily assigned without cogent physical grounds; finally, let the data constrain $\theta^{*}$ to a small range close to [math]. It is then easy to see that $B_{21}\propto 1/a$.

3 Cross validation

An alternative framework for model comparison is offered by measures of predictive performance, which quantify how well a model that has been fit to dataset $y$ can predict yet-to-be-observed data $\tilde{y}$ (Gelman et al., 2013). These measures penalize overfitting by construction, because a model that conforms to contingent features in $y$ will usually do worse in fitting $\tilde{y}$. In a Bayesian setting, a commonly adopted measure of predictive performance is the log predictive density for the new data $\tilde{y}$, as induced by the posterior $p(\theta_{M}|y)$:

[TABLE]

Ideally, we would average $\log p(\tilde{y}|y;M)$ over the true distribution of future data $\tilde{y}$; doing so is however seldom possible. In practice, we can: a) estimate within-sample predictive accuracy using the data $y$ that we already have, by applying corrections for the overfitting bias, as in the various “information criteria” (Gelman et al., 2013); b) evaluate out-of-sample predictive accuracy on one or more holdout datasets that were not used to infer parameter posteriors. The latter approach is known as cross-validation, and we will pursue it for PTA data in the rest of this article.

Specifically, we adopt $k$-fold cross validation as follows:

We divide the dataset $y$ randomly in $k$ exclusive subsets $y^{(k)}$ (the testing datasets); 2. 2.

For each $k$, we derive the posterior $p(\theta_{M}|y^{(-k)})$, where $y^{(-k)}=\cup_{j\neq k}y^{(j)}$ is the training dataset corresponding by omitting $y^{(k)}$. We represent posteriors as sequences $\{\theta^{(-k)}_{M,i}\}$ of $N$ quasi-independent samples, obtained by Monte Carlo methods; 3. 3.

For each $k$, we evaluate the log predictive density $\log p(y^{(k)}|y^{(-k)})$, given by

[TABLE] 4. 4.

We repeat this procedure for every model under consideration, and then compare the respective log predictive densities, averaged over the $k$ repetitions. The variance of the densities is a measure of their statistical uncertainty.

In Figure 1, we exemplify this process with a simple toy problem: a sinusoidal signal parametrized by amplitude, frequency, and phase, as modeled by that very model and by an expanded model that includes a linear frequency drift. The more complicated model results in significantly lower log predictive density, demonstrating that cross validation can recognize and reject overfitting.

4 Cross validation for single-pulsar noise modeling

In current practice (see van Haasteren & Vallisneri 2014 for a recent review), probabilistic noise models for pulsars are built as the sum of a number of Gaussian processes (GPs; Rasmussen & Williams 2006) representing all sources of correlated noise: errors in the parameters of the deterministic timing model, pulsar-rotation irregularities, dispersion-measure variations along the pulse propagation path, jitter-like noise in multifrequency observations, and more. The TOAs are also subject to “white” radiometer measurement noise, conceptualized as independent and heteroskedastic normal variates. Both the GPs and measurement noise are governed by a set of hyperparameters (e.g., the amplitude and spectral slope of delays due to rotation irregularities) that are estimated from PTA datasets.

In keeping with basis–kernel duality for GPs, the model likelihood can be written in two complementary ways. In hierarchical form, the likelihood is given by

[TABLE]

where $y$ is the vector of $n$ timing residuals obtained by subtracting the best-fit timing model from the observed pulse times of arrival; the $n\times m$ matrix $F_{\mathrm{GP}}$ collects the $m$ GP basis vectors, and the $c_{\mathrm{GP}}$ are the corresponding weights (or coefficients); $N$ (a function of the hyperparameters $\eta_{N}$) is a diagonal matrix expressing measurement-noise variance; and $\Phi_{\mathrm{GP}}$ (a function of the hyperparameters $\eta_{\mathrm{GP}}$) represents the normal priors for the GP weights. In marginalized form, we eliminate the dependence on the GP weights by integrating over them (van Haasteren & Vallisneri, 2014):

[TABLE]

The marginalized form is usually employed for the Monte Carlo exploration of hyperparameter posteriors. The GP weights can still be characterized by way of their conditional posterior given the data and the hyperparameters, which follows the jointly normal distribution $p(c_{\mathrm{GP}}|y;\eta_{N},\eta_{\mathrm{GP}})=\mathcal{N}(\bar{c},\Sigma)$ with mean

[TABLE]

and covariance

[TABLE]

Armed with Eqs. (4)–(7), we perform steps 1–3 of $k$-fold cross-validation for a single-pulsar dataset as follows:

We partition the timing residuals $y$ into exclusive testing datasets $y^{(k)}$, making sure that each training dataset $y^{(-k)}$ depends on all weights and hyperparameters. For instance, for dispersion-measure variations described as piecewise-constant “DMX” functions, each DMX epoch must be populated by at least one residual in every $y^{(k)}$; likewise, for “EFAC” measurement noise that is rescaled differently in each radio backend, each backend must be represented in every $y^{(k)}$; 2. 2.

We sample the marginalized hyperparameter posterior $p(\eta_{N},\eta_{\mathrm{GP}}|y^{(-k)})$ [proportional to Eq. (5) times prior $p(\eta_{N},\eta_{\mathrm{GP}})$], using the PTA data-analysis package Enterprise†††https://github.com/nanograv/enterprise and the Markov Chain Monte Carlo sampler PTMCMCSampler;‡‡‡https://github.com/jellis18/PTMCMCSampler 3. 3.

For each of the $N$ quasi-independent $(\eta^{(-k)}_{N,i},\eta^{(-k)}_{\mathrm{GP},i})$ obtained at step 2, we draw $P$ weight vectors $\{c^{(-k)}_{\mathrm{GP},ij}\}$ from their conditional distribution [Eqs. (6)–(7)], then we evaluate the log predictive density by averaging the hierarchical likelihood (4) over the $N\times P$ triples $(\eta^{(-k)}_{N,i},\eta^{(-k)}_{\mathrm{GP},i},c^{(-k)}_{\mathrm{GP},ij})$.

The “representation” condition imposed in step 1 is necessary so that the parameters for which we derive posteriors at step 2 fully specify the model’s prediction for each testing dataset; this prediction is used in step 3 to evaluate the predictive density. Another important technical subtlety is that the GP weights must conserve the same identity across all $y^{(-k)}$ (for example, Fourier coefficients for the same set of frequencies in the case of correlated spin noise); by contrast, the GP basis vectors would change in value, because they refer to different time-of-arrival measurements (continuing our example, the basis elements would be sines and cosines of the same frequencies for all subdatasets, but would be evaluated at different times).

We note also that in step 2 we could have used the hierarchical likelihood to sample the full parameter set $(\eta_{N},\eta_{\mathrm{GP}},c_{\mathrm{GP}})$, avoiding the conditional $c_{\mathrm{GP}}$ draws at step 3. However, the hierarchical likelihood is considerably harder to explore stochastically. Also, in step 3 the sum over the $c^{(-k)}_{\mathrm{GP},ij}$ for each $i$ could be replaced by analytical integration in terms of $\bar{c}$ and $\Sigma$, at the cost of some algebraic complication (see below for a related development).

5 Results

To demonstrate how cross validation can be applied to PTA noise-model selection, we first consider the simulated TOA residuals for pulsar J1713+0747 from the $1^{\mathrm{st}}$ International Pulsar Timing Array (IPTA) Data Challenge (Hazboun et al., 2018).§§§http://ipta4gw.org/data-challenge These residuals (130 values at 14-day cadence) are generated from a simplified deterministic timing model and a simple noise model: the timing model parameters describe the intrinsic-spin, astrometry, and binary-orbit properties of the pulsar; the noise model includes white measurement noise (described by “EFAC” and “EQUAD” parameters) and red correlated spin noise, specified by a power-law spectrum (Arzoumanian et al., 2018b)

[TABLE]

where $A_{\mathrm{red}}$ is the amplitude of the red-noise process in units of $\mu s\times\mathrm{yr}^{1/2}$, $\gamma_{\mathrm{red}}$ is its spectral index (set to 13/3 to generate the data), and $f_{\mathrm{yr}}=1\,\mathrm{yr}^{-1}$. We compare power-law red-noise models with different $\gamma_{\mathrm{red}}$ by evaluating their respective cross-validation predictive densities, shown in Figure 2 for $\gamma_{\mathrm{red}}$ ranging from 1 to 6. It is clear that a value larger than 4 (and consistent with 13/3) is preferred.

Moving on to real PTA datasets, we perform two-fold cross validation on the TOA residuals of pulsar J1713+0747 from the NANOGrav 11-yr data release (Arzoumanian et al., 2018a). As apparent in Figure 3 (from Arzoumanian et al. 2018a), around 2009 the residuals underwent a dispersion-measure (DM) dip (i.e., an apparent decrease in the electron density experienced by radio pulses traveling to Earth, resulting in reduced “fanning” across frequencies). We compare two noise models: the first (GP+dip) represents DM with a Gaussian process, but it includes also a transient feature with exponential decay to represent the dip; the second (GP-only) represents DM with the Gaussian process alone. In keeping with the representation condition introduced above, when we select the training dataset we need to make sure that it includes a sufficient number of residuals around the dip. To achieve this, we build the training dataset by randomly selecting half of the residuals within a one-year window centered around 2009 (see Figure 3), as well as half of the residuals outside the window.

We perform two-fold cross validation 32 times with different random data partitions. The GP+dip model yields consistently higher predictive densities than the GP-only model, with delta log density $54.8^{+584.2}_{-19.9}$ (quoted as median augmented by the 90% interquantile range). The stronger predictive performance of the GP+dip model offers statistical evidence that the dip is a real physical feature. This is confirmed by inspecting the reconstructed DM. The gray bands in Figure 4 show the 90% interquantile range of reconstructed DM over 100 posterior draws of the model parameters in a single cross-validation run, while the red line traces the “true” DM at each epoch, as obtained by fitting independent dispersion values to multifrequency data.¶¶¶In PTA jargon, the red line shows the best-fit “DMX” parameters. It is apparent that the GP+dip model captures the DM transient more accurately, and that it follows the overall evolution of DM with smaller variance. Both conditions result in higher predictive density.

6 Cross validation for GW detection in multi-pulsar datasets.

To perform cross validation on multi-pulsar datasets and models (the latter possibly including a common GP describing the correlated delays induced by GWs for each pulsar), we may proceed without change if we satisfy the representation condition: that is, if we partition the multi-pulsar residual vector, $Y$, into testing datasets $Y^{(k)}$ in such a way that all parameters of every pulsar are represented in every training dataset $Y^{(-k)}$. However, it seems natural to partition $Y$ instead into subsets that correspond to individual pulsars, or groups of individual pulsars. Such an arrangement may allow us to identify pulsars contaminated by pathological observations, pulsars that are poorly described by the noise model chosen for them, or pulsars that are too noisy to contribute to GW inference. We next discuss how to proceed for this more general partitioning. We describe separately the case of deterministic and stochastic GWs.

For GWs described by deterministic models (e.g., isolated supermassive black-hole binaries), individual pulsars are described by the likelihoods of Eqs. (4) and (5), with the replacement $y\rightarrow Y^{(a)}-d^{(a)}(\theta_{\mathrm{GW}})$, where $Y^{(a)}$ is the vector of timing residuals for pulsar $a$, where the $\theta_{\mathrm{GW}}$ describe the GW parameters (a common set for all pulsars), and where the $d^{(a)}(\theta_{\mathrm{GW}})$ are the delays induced by the GWs on pulsar $a$. In this case, cross validation would proceed as follows: a) for each $Y^{(-k)}$ we would sample the joint posterior distribution of the $\theta_{\mathrm{GW}}$ and of the noise hyperparameters $\eta_{N,\mathrm{GP}}^{(-k)}$ describing the pulsars represented in $Y^{(-k)}$; b) for each corresponding $Y^{(k)}$, we would evaluate the log marginalized predictive likelihood

[TABLE]

where the $\theta^{(-k)}_{\mathrm{GW,i}}$ are Markov Chain (sub-)samples from the training posterior $p(\theta_{\mathrm{GW}},\eta_{N,\mathrm{GP}}^{(-k)}|Y^{(-k)})$, and where we have dropped the $\eta_{N,\mathrm{GP}}$ suffix for compactness. From an implementation standpoint, the nested sum/integral in Eq. (9) may require a dedicated stochastic algorithm similar to those employed to evaluate the Bayesian evidence (Trotta, 2008).

It is important to notice that Eq. (9) depends directly on the noise-hyperparameter priors $p(\eta^{(k)})$, which invalidates some of our motivation for computing predictive likelihoods in the first place. In practice, we may worry that we cannot compare predictive likelihoods for different $Y^{(k)}$ because they have different prior “calibrations.” With respect to this objection, it seems then natural to consider the ratio of Eq. (9) to the noise-only evidence $\int p(Y^{(k)}|\eta^{(k)})p(\eta^{(k)})\,\mathrm{d}\eta^{(k)}$ (which in fact factorizes over the pulsars in $Y^{(k)}$). We leave to future work the exploration of marginalized predictive likelihoods as GW detection statistics, as well as the development of an efficient sampling method for Eq. (9).

Last, for stochastic GWs described by their spectrum and by their correlations across pulsars, we need to account not only for the common GW parameters $\theta_{\mathrm{GW}}$, but also for the correlations between the GW GP weights in each pulsar. To sketch the mathematical structure of the problem with more readable notation, let us consider the case of training on pulsar 1 and testing on pulsar 2; formulas generalize readily to $k$-fold–validation testing pairs $(Y^{(-k)},Y^{(k)})$. We do as follows: We first obtain samples $(\eta^{(1)}_{i},\theta_{\mathrm{GW},i})$ from the posterior $p(\eta^{(1)},\theta_{\mathrm{GW}}|Y^{(1)})$; we then evaluate the log marginalized predictive likelihood in the form

[TABLE]

The integral over the $c^{(1)}_{\mathrm{GW}}$ can be performed analytically, and the resulting conditional prior for the $c^{(1)}_{\mathrm{GW}}$ expressed as

[TABLE]

with

[TABLE]

and

[TABLE]

where $\bar{c}^{(1)}(Y^{(1)};\eta^{(1)}_{i},\theta_{\mathrm{GW},i})$ and $\Sigma^{(1)}(\eta^{(1)}_{i},\theta_{\mathrm{GW},i})$ are given in Eqs. (6) and (7), and the $\Phi_{ij}(\theta_{\mathrm{GW},i})$ denote the blocks of the joint normal prior for the GW GP weights. The integral over the $c^{(2)}_{\mathrm{GW}}$ in Eq. (10) may also be performed analytically by way of Eq. (5), with the replacement $y\rightarrow Y^{(2)}-F^{(2)}_{\mathrm{GW}}\bar{c}^{(2)|(1)}$ and $\Phi_{\mathrm{GW}}\rightarrow\Sigma^{(2)|(1)}$ (in this notation, $F^{(2)}_{\mathrm{GW}}$ and $\Phi_{\mathrm{GW}}$ represent the GW blocks of $F_{\mathrm{GP}}$ and $\Phi_{\mathrm{GP}}$). Again, we leave to future work the investigation of marginalized predictive likelihoods as detection statistics for stochastic GWs in PTA data.

7 Discussion

The development of sophisticated noise models is of paramount importance to ongoing PTA searches for nanohertz GWs, especially so because much of the PTA “instrument” was not engineered by humans. Contrast LIGO’s precisely fabricated arms with a PTA Earth–pulsar “arm”, which comprises a radiotelescope, a distant ($\sim$ kpc) pulsar, and the expanse of spacetime in between. The telescope admits some experimental control, since we can test its signal response and mitigate noise in the receiver, but we still have to contend with poorly constrained physical processes in pulsar interiors and emission regions (Lasky et al., 2015), not to mention dispersive effects as radio pulses propagate through the ionized interstellar medium to the Earth (Cordes & Shannon, 2010; Lam et al., 2019, 2018; Lam et al., 2016). Thus PTA noise models must be well motivated physically, and yet flexible enough to accommodate unknown unknowns, if they are to allow the identification of subtle GW-induced delays.

Current techniques to test the relative aptness and robustness of noise models include basic checks (i.e., evaluating $\chi^{2}$ fit residuals under different models), frequentist approaches (computing receiver operating characteristic curves to maximize the signal detection probability at a given false-alarm probability), as well as Bayesian model comparison (Taylor et al., 2013; Taylor et al., 2017; Sampson et al., 2015; Taylor et al., 2014; Cornish & Sampson, 2016). In this last technique, we marginalize likelihoods for different models over their respective prior volumes, and use the ratios of marginal likelihoods (usually under the assumption of equal prior probabilities for each model) to make statements about the posterior odds with which one model is favored over the other. This approach is perfectly valid, but it can be troubled by practical issues such as accurately integrating into the tails of the likelihood distribution, as well as assigning the ranges of parameter priors in the first place.

Bayesian cross validation addresses some of these issues by partitioning datasets into training and testing samples. Models are conditioned on the training set, producing posterior probability distributions for the parameters, over which we average the likelihood of the testing set. The cross-validation score is then the probability of the test data under the trained model. For reasonably informative training data, posteriors will be more compact than the priors, shielding the integration from ad hoc prior choices. To evaluate the integrals, it is convenient to average the likelihood over posterior samples from a conventional MCMC analysis of the training set.

In this article we argued for the power of Bayesian cross validation as applied to PTA data. Our case studies included the spectral characterization of a GW background, and the modeling of dispersive noise from the interstellar medium. For the former, we performed two-fold cross-validation on a dataset from the $1^{\mathrm{st}}$ IPTA mock data challenge (Hazboun et al., 2018), showing that the data favored a correlated process consistent with a stochastic GW background from supermassive binary black holes. For the latter, two-fold cross-validation on $11$ years of NANOGrav data for pulsar J$1713$+$0747$ illustrated the necessity to include a transient dispersive noise feature around the year $2009$, consistent with a void in the electron density along the line of sight (Lentati et al., 2016). We also introduced a formalism for multi-pulsar cross validation, where GW models conditioned on training data from a subarray are assessed for predictive performance on left-out pulsars. In future work we will investigate this approach as a tool to validate claims of GW detections in real PTA datasets.

We expect cross validation to be similarly useful in analyzing data from other GW detectors. For instance, within the global network of ground-based GW interferometers, support for a GW signal in one detector could be validated using data from a different widely separated detector. Furthermore, signals such as GW150914, which could have been observed by a LISA-like detector years before it was seen by LIGO (Sesana, 2016), raise the prospect of multi-band, multi-detector cross-validation of GW signals.

Acknowledgments.

We thank Joseph Simon and Michael Lam for discussions regarding the dispersion-measure variation of PSR J$1713$+$0747$. This research was performed in part using the Zwicky computer cluster at Caltech supported by NSF under MRI-R2 award No. PHY0960291 and by the Sherman Fairchild Foundation. Portions of this research were carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. This work was supported in part by National Science Foundation Grant No. PHYS-1066293 and by the hospitality of the Aspen Center for Physics. MV was supported by the Jet Propulsion Laboratory RTD program. SRT was supported by the NANOGrav National Science Foundation Physics Frontier Center, award number 1430284. Parts of this work were carried out at the Jet Propulsion Laboratory, California Institute of Technology, under contract to the National Aeronautics and Space Administration. Copyright 2019 California Institute of Technology. Government sponsorship acknowledged.