NANOGrav Collaboration

# The NANOGrav 15-year Gravitational-Wave Background Analysis Pipeline

Aaron D. Johnson  0000-0002-7445-8423 Center for Gravitation, Cosmology and
Astrophysics, Department of Physics, University of Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  Division of Physics, Mathematics, and
Astronomy, California Institute of Technology, Pasadena, CA 91125, USA
Patrick M. Meyers  0000-0002-2689-0190 Division of Physics, Mathematics, and
Astronomy, California Institute of Technology, Pasadena, CA 91125, USA  Paul
T. Baker  0000-0003-2745-753X Department of Physics and Astronomy, Widener
University, One University Place, Chester, PA 19013, USA  Neil J. Cornish
0000-0002-7435-0869 Department of Physics, Montana State University, Bozeman,
MT 59717, USA  Jeffrey S. Hazboun  0000-0003-2742-3321 Department of Physics,
Oregon State University, Corvallis, OR 97331, USA  Tyson B. Littenberg
0000-0002-9574-578X NASA Marshall Space Flight Center, Huntsville, AL 35812,
USA  Joseph D. Romano  0000-0003-4915-3246 Department of Physics, Texas Tech
University, Box 41051, Lubbock, TX 79409, USA  Stephen R. Taylor
0000-0003-0264-1453 Department of Physics and Astronomy, Vanderbilt
University, 2301 Vanderbilt Place, Nashville, TN 37235, USA  Michele
Vallisneri  0000-0002-4162-0033 Jet Propulsion Laboratory, California
Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA
Division of Physics, Mathematics, and Astronomy, California Institute of
Technology, Pasadena, CA 91125, USA  Sarah J. Vigeland  0000-0003-4700-9072
Center for Gravitation, Cosmology and Astrophysics, Department of Physics,
University of Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  Ken D. Olum  0000-0002-2027-3714
Institute of Cosmology, Department of Physics and Astronomy, Tufts University,
Medford, MA 02155, USA  Xavier Siemens  0000-0002-7778-2990 Department of
Physics, Oregon State University, Corvallis, OR 97331, USA  Center for
Gravitation, Cosmology and Astrophysics, Department of Physics, University of
Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  Justin A. Ellis  Department of
Astronomy, University of Maryland, College Park, MD 20742  Center for Research
and Exploration in Space Science and Technology, NASA/GSFC, Greenbelt, MD
20771  NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA  Infinia ML,
202 Rigsbee Avenue, Durham NC, 27701  Rutger van Haasteren
0000-0002-6428-2620 Max-Planck-Institut für Gravitationsphysik (Albert-
Einstein-Institut), Callinstrasse 38, D-30167, Hannover, Germany  Sophie
Hourihane  missing Division of Physics, Mathematics, and Astronomy, California
Institute of Technology, Pasadena, CA 91125, USA  Gabriella Agazie
0000-0001-5134-3925 Center for Gravitation, Cosmology and Astrophysics,
Department of Physics, University of Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  Akash Anumarlapudi
0000-0002-8935-9882 Center for Gravitation, Cosmology and Astrophysics,
Department of Physics, University of Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  Anne M. Archibald  0000-0003-0638-3340
Newcastle University, NE1 7RU, UK  Zaven Arzoumanian  missing X-Ray
Astrophysics Laboratory, NASA Goddard Space Flight Center, Code 662,
Greenbelt, MD 20771, USA  Laura Blecha  0000-0002-2183-1087 Physics
Department, University of Florida, Gainesville, FL 32611, USA  Adam Brazier
0000-0001-6341-7178 Cornell Center for Astrophysics and Planetary Science and
Department of Astronomy, Cornell University, Ithaca, NY 14853, USA  Cornell
Center for Advanced Computing, Cornell University, Ithaca, NY 14853, USA  Paul
R. Brook  0000-0003-3053-6538 Institute for Gravitational Wave Astronomy and
School of Physics and Astronomy, University of Birmingham, Edgbaston,
Birmingham B15 2TT, UK  Sarah Burke-Spolaor  0000-0003-4052-7838 Department of
Physics and Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV
26506, USA  Center for Gravitational Waves and Cosmology, West Virginia
University, Chestnut Ridge Research Building, Morgantown, WV 26505, USA  Bence
Bécsy  0000-0003-0909-5563 Department of Physics, Oregon State University,
Corvallis, OR 97331, USA  J. Andrew Casey-Clyde  0000-0002-5557-4007
Department of Physics, University of Connecticut, 196 Auditorium Road, U-3046,
Storrs, CT 06269-3046, USA  Maria Charisi  0000-0003-3579-2522 Department of
Physics and Astronomy, Vanderbilt University, 2301 Vanderbilt Place,
Nashville, TN 37235, USA  Shami Chatterjee  0000-0002-2878-1502 Cornell Center
for Astrophysics and Planetary Science and Department of Astronomy, Cornell
University, Ithaca, NY 14853, USA  Katerina Chatziioannou  missing Division of
Physics, Mathematics, and Astronomy, California Institute of Technology,
Pasadena, CA 91125, USA  Tyler Cohen  0000-0001-7587-5483 Department of
Physics, New Mexico Institute of Mining and Technology, 801 Leroy Place,
Socorro, NM 87801, USA  James M. Cordes  0000-0002-4049-1882 Cornell Center
for Astrophysics and Planetary Science and Department of Astronomy, Cornell
University, Ithaca, NY 14853, USA  Fronefield Crawford  0000-0002-2578-0360
Department of Physics and Astronomy, Franklin & Marshall College, P.O. Box
3003, Lancaster, PA 17604, USA  H. Thankful Cromartie  0000-0002-6039-692X
NASA Hubble Fellowship: Einstein Postdoctoral Fellow  Cornell Center for
Astrophysics and Planetary Science and Department of Astronomy, Cornell
University, Ithaca, NY 14853, USA  Kathryn Crowter  0000-0002-1529-5169
Department of Physics and Astronomy, University of British Columbia, 6224
Agricultural Road, Vancouver, BC V6T 1Z1, Canada  Megan E. DeCesar
0000-0002-2185-1790 George Mason University, resident at the Naval Research
Laboratory, Washington, DC 20375, USA  Paul B. Demorest  0000-0002-6664-965X
National Radio Astronomy Observatory, 1003 Lopezville Rd., Socorro, NM 87801,
USA  Timothy Dolch  0000-0001-8885-6388 Department of Physics, Hillsdale
College, 33 E. College Street, Hillsdale, MI 49242, USA  Eureka Scientific,
2452 Delmer Street, Suite 100, Oakland, CA 94602-3017, USA  Brendan Drachler
missing School of Physics and Astronomy, Rochester Institute of Technology,
Rochester, NY 14623, USA  Laboratory for Multiwavelength Astrophysics,
Rochester Institute of Technology, Rochester, NY 14623, USA  Elizabeth C.
Ferrara  0000-0001-7828-7708 Department of Astronomy, University of Maryland,
College Park, MD 20742  Center for Research and Exploration in Space Science
and Technology, NASA/GSFC, Greenbelt, MD 20771  NASA Goddard Space Flight
Center, Greenbelt, MD 20771, USA  William Fiore  0000-0001-5645-5336
Department of Physics and Astronomy, West Virginia University, P.O. Box 6315,
Morgantown, WV 26506, USA  Center for Gravitational Waves and Cosmology, West
Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505,
USA  Emmanuel Fonseca  0000-0001-8384-5049 Department of Physics and
Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA
Center for Gravitational Waves and Cosmology, West Virginia University,
Chestnut Ridge Research Building, Morgantown, WV 26505, USA  Gabriel E.
Freedman  0000-0001-7624-4616 Center for Gravitation, Cosmology and
Astrophysics, Department of Physics, University of Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  Nate Garver-Daniels
0000-0001-6166-9646 Department of Physics and Astronomy, West Virginia
University, P.O. Box 6315, Morgantown, WV 26506, USA  Center for Gravitational
Waves and Cosmology, West Virginia University, Chestnut Ridge Research
Building, Morgantown, WV 26505, USA  Peter A. Gentile  0000-0001-8158-683X
Department of Physics and Astronomy, West Virginia University, P.O. Box 6315,
Morgantown, WV 26506, USA  Center for Gravitational Waves and Cosmology, West
Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505,
USA  Joseph Glaser  0000-0003-4090-9780 Department of Physics and Astronomy,
West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA  Center for
Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge
Research Building, Morgantown, WV 26505, USA  Deborah C. Good
0000-0003-1884-348X Department of Physics, University of Connecticut, 196
Auditorium Road, U-3046, Storrs, CT 06269-3046, USA  Center for Computational
Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA
Kayhan Gültekin  0000-0002-1146-0198 Department of Astronomy and Astrophysics,
University of Michigan, Ann Arbor, MI 48109, USA  Ross J. Jennings
0000-0003-1082-2342 NANOGrav Physics Frontiers Center Postdoctoral Fellow
Department of Physics and Astronomy, West Virginia University, P.O. Box 6315,
Morgantown, WV 26506, USA  Center for Gravitational Waves and Cosmology, West
Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505,
USA  Megan L. Jones  0000-0001-6607-3710 Center for Gravitation, Cosmology and
Astrophysics, Department of Physics, University of Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  Andrew R. Kaiser  0000-0002-3654-980X
Department of Physics and Astronomy, West Virginia University, P.O. Box 6315,
Morgantown, WV 26506, USA  Center for Gravitational Waves and Cosmology, West
Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505,
USA  David L. Kaplan  0000-0001-6295-2881 Center for Gravitation, Cosmology
and Astrophysics, Department of Physics, University of Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  Luke Zoltan Kelley
0000-0002-6625-6450 Department of Astronomy, University of California,
Berkeley, 501 Campbell Hall #3411, Berkeley, CA 94720, USA  Matthew Kerr
0000-0002-0893-4073 Space Science Division, Naval Research Laboratory,
Washington, DC 20375-5352, USA  Joey S. Key  0000-0003-0123-7600 University of
Washington Bothell, 18115 Campus Way NE, Bothell, WA 98011, USA  Nima Laal
0000-0002-9197-7604 Department of Physics, Oregon State University, Corvallis,
OR 97331, USA  Michael T. Lam  0000-0003-0721-651X School of Physics and
Astronomy, Rochester Institute of Technology, Rochester, NY 14623, USA
Laboratory for Multiwavelength Astrophysics, Rochester Institute of
Technology, Rochester, NY 14623, USA  William G. Lamb  0000-0003-1096-4156
Department of Physics and Astronomy, Vanderbilt University, 2301 Vanderbilt
Place, Nashville, TN 37235, USA  T. Joseph W. Lazio  missing Jet Propulsion
Laboratory, California Institute of Technology, 4800 Oak Grove Drive,
Pasadena, CA 91109, USA  Natalia Lewandowska  0000-0003-0771-6581 Department
of Physics, State University of New York at Oswego, Oswego, NY, 13126, USA
Tingting Liu  0000-0001-5766-4287 Department of Physics and Astronomy, West
Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA  Center for
Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge
Research Building, Morgantown, WV 26505, USA  Duncan R. Lorimer
0000-0003-1301-966X Department of Physics and Astronomy, West Virginia
University, P.O. Box 6315, Morgantown, WV 26506, USA  Center for Gravitational
Waves and Cosmology, West Virginia University, Chestnut Ridge Research
Building, Morgantown, WV 26505, USA  Jing Luo  0000-0001-5373-5914 Deceased
Department of Astronomy & Astrophysics, University of Toronto, 50 Saint George
Street, Toronto, ON M5S 3H4, Canada  Ryan S. Lynch  0000-0001-5229-7430 Green
Bank Observatory, P.O. Box 2, Green Bank, WV 24944, USA  Chung-Pei Ma
0000-0002-4430-102X Department of Astronomy, University of California,
Berkeley, 501 Campbell Hall #3411, Berkeley, CA 94720, USA  Department of
Physics, University of California, Berkeley, CA 94720, USA  Dustin R. Madison
0000-0003-2285-0404 Department of Physics, University of the Pacific, 3601
Pacific Avenue, Stockton, CA 95211, USA  Alexander McEwen  0000-0001-5481-7559
Center for Gravitation, Cosmology and Astrophysics, Department of Physics,
University of Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  James W. McKee  0000-0002-2885-8485
E.A. Milne Centre for Astrophysics, University of Hull, Cottingham Road,
Kingston-upon-Hull, HU6 7RX, UK  Centre of Excellence for Data Science,
Artificial Intelligence and Modelling (DAIM), University of Hull, Cottingham
Road, Kingston-upon-Hull, HU6 7RX, UK  Maura A. McLaughlin
0000-0001-7697-7422 Department of Physics and Astronomy, West Virginia
University, P.O. Box 6315, Morgantown, WV 26506, USA  Center for Gravitational
Waves and Cosmology, West Virginia University, Chestnut Ridge Research
Building, Morgantown, WV 26505, USA  Natasha McMann  0000-0002-4642-1260
Department of Physics and Astronomy, Vanderbilt University, 2301 Vanderbilt
Place, Nashville, TN 37235, USA  Bradley W. Meyers  0000-0001-8845-1225
Department of Physics and Astronomy, University of British Columbia, 6224
Agricultural Road, Vancouver, BC V6T 1Z1, Canada  International Centre for
Radio Astronomy Research, Curtin University, Bentley, WA 6102, Australia
Chiara M. F. Mingarelli  0000-0002-4307-1322 Center for Computational
Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA
Department of Physics, University of Connecticut, 196 Auditorium Road, U-3046,
Storrs, CT 06269-3046, USA  Department of Physics, Yale University, New Haven,
CT 06520, USA  Andrea Mitridate  0000-0003-2898-5844 Deutsches Elektronen-
Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany  Cherry Ng
0000-0002-3616-5160 Dunlap Institute for Astronomy and Astrophysics,
University of Toronto, 50 St. George St., Toronto, ON M5S 3H4, Canada  David
J. Nice  0000-0002-6709-2566 Department of Physics, Lafayette College, Easton,
PA 18042, USA  Stella Koch Ocker  0000-0002-4941-5333 Cornell Center for
Astrophysics and Planetary Science and Department of Astronomy, Cornell
University, Ithaca, NY 14853, USA  Timothy T. Pennucci  0000-0001-5465-2889
Institute of Physics and Astronomy, Eötvös Loránd University, Pázmány P. s.
1/A, 1117 Budapest, Hungary  Benetge B. P. Perera  0000-0002-8509-5947 Arecibo
Observatory, HC3 Box 53995, Arecibo, PR 00612, USA  Nihan S. Pol
0000-0002-8826-1285 Department of Physics and Astronomy, Vanderbilt
University, 2301 Vanderbilt Place, Nashville, TN 37235, USA  Henri A. Radovan
0000-0002-2074-4360 Department of Physics, University of Puerto Rico,
Mayagüez, PR 00681, USA  Scott M. Ransom  0000-0001-5799-9714 National Radio
Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA  Paul
S. Ray  0000-0002-5297-5278 Space Science Division, Naval Research Laboratory,
Washington, DC 20375-5352, USA  Shashwat C. Sardesai  0009-0006-5476-3603
Center for Gravitation, Cosmology and Astrophysics, Department of Physics,
University of Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  Carl Schmiedekamp  0000-0002-1283-2184
Department of Physics, Penn State Abington, Abington, PA 19001, USA  Ann
Schmiedekamp  0000-0003-4391-936X Department of Physics, Penn State Abington,
Abington, PA 19001, USA  Kai Schmitz  0000-0003-2807-6472 Institute for
Theoretical Physics, University of Münster, 48149 Münster, Germany  Brent J.
Shapiro-Albert  0000-0002-7283-1124 Department of Physics and Astronomy, West
Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA  Center for
Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge
Research Building, Morgantown, WV 26505, USA  Giant Army, 915A 17th Ave,
Seattle WA 98122  Joseph Simon  0000-0003-1407-6607 NSF Astronomy and
Astrophysics Postdoctoral Fellow  Department of Astrophysical and Planetary
Sciences, University of Colorado, Boulder, CO 80309, USA  Magdalena S. Siwek
0000-0002-1530-9778 Center for Astrophysics, Harvard University, 60 Garden St,
Cambridge, MA 02138  Ingrid H. Stairs  0000-0001-9784-8670 Department of
Physics and Astronomy, University of British Columbia, 6224 Agricultural Road,
Vancouver, BC V6T 1Z1, Canada  Daniel R. Stinebring  0000-0002-1797-3277
Department of Physics and Astronomy, Oberlin College, Oberlin, OH 44074, USA
Kevin Stovall  0000-0002-7261-594X National Radio Astronomy Observatory, 1003
Lopezville Rd., Socorro, NM 87801, USA  Abhimanyu Susobhanan
0000-0002-2820-0931 Center for Gravitation, Cosmology and Astrophysics,
Department of Physics, University of Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  Joseph K. Swiggum  0000-0002-1075-3837
NANOGrav Physics Frontiers Center Postdoctoral Fellow  Department of Physics,
Lafayette College, Easton, PA 18042, USA  Jacob E. Turner  0000-0002-2451-7288
Department of Physics and Astronomy, West Virginia University, P.O. Box 6315,
Morgantown, WV 26506, USA  Center for Gravitational Waves and Cosmology, West
Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505,
USA  Caner Unal  0000-0001-8800-0192 Department of Physics, Ben-Gurion
University of the Negev, Be’er Sheva 84105, Israel  Feza Gursey Institute,
Bogazici University, Kandilli, 34684, Istanbul, Turkey  Haley M. Wahl
0000-0001-9678-0299 Department of Physics and Astronomy, West Virginia
University, P.O. Box 6315, Morgantown, WV 26506, USA  Center for Gravitational
Waves and Cosmology, West Virginia University, Chestnut Ridge Research
Building, Morgantown, WV 26505, USA  Caitlin A. Witt  0000-0002-6020-9274
Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA),
Northwestern University, Evanston, IL 60208  Adler Planetarium, 1300 S.
DuSable Lake Shore Dr., Chicago, IL 60605, USA  Olivia Young
0000-0002-0883-0688 School of Physics and Astronomy, Rochester Institute of
Technology, Rochester, NY 14623, USA  Laboratory for Multiwavelength
Astrophysics, Rochester Institute of Technology, Rochester, NY 14623, USA

###### Abstract

This paper presents rigorous tests of pulsar timing array methods and
software, examining their consistency across a wide range of injected
parameters and signal strength. We discuss updates to the 15-year isotropic
gravitational-wave background analyses and their corresponding code
representations. Descriptions of the internal structure of the flagship
algorithms Enterprise and PTMCMCSampler are given to facilitate understanding
of the PTA likelihood structure, how models are built, and what methods are
currently used in sampling the high-dimensional PTA parameter space. We
introduce a novel version of the PTA likelihood that uses a two-step
marginalization procedure that performs much faster when the white noise
parameters remain fixed. We perform stringent tests of consistency and
correctness of the Bayesian and frequentist analysis software. For the
Bayesian analysis, we test prior recovery, injection recovery, and Bayes
factors. For the frequentist analysis, we test that the cross-correlation-
based optimal statistic, when modified to account for a non-negligible
gravitational-wave background, accurately recovers the amplitude of the
background. We also summarize recent advances and tests performed on the
optimal statistic in the literature from both GWB detection and parameter
estimation perspectives. The tests presented here validate current and future
analyses of PTA data.

††preprint: APS/123-QED

##  I Introduction

Since the first detection of gravitational waves (GWs) from a stellar mass
black hole binary in 2015 [1], the field of GW astronomy has flourished with
dozens more detections of transient signals [2]. Besides these signals at
${\cal{O}}(100)$ Hz frequencies, detected via ground-based laser
interferometry, other methods can detect GWs across a wide range of
frequencies. Pulsar timing arrays (PTAs) create a galactic-scale GW detector
using radio telescopes to collect pulse times of arrival (TOAs) from an array
of millisecond pulsars. These pulsars exhibit exceptional long-timescale
arrival-time stability, allowing PTAs to use them as galactic scale clocks
that are sensitive to perturbations at frequencies $1$–$100$ nHz. The expected
target signal is the stochastic GW background (GWB), possibly from an ensemble
of merging supermassive black-hole binaries that could result from galactic
mergers [3].

Previously, the European Pulsar Timing Array (EPTA) [4], the Parkes Pulsar
Timing Array in Australia (PPTA) [5], and the North American Nanohertz
Observatory for Gravitational waves (NANOGrav) [6, 7] all reported detection
of a red-noise process with common spectrum among pulsars, but no evidence
either way for inter-pulsar correlations [8, 9, 10]. Such a process is known
as CURN, for Common-spectrum Uncorrelated Red Noise. The International Pulsar
Timing Array (IPTA) [11] consists of these collaborations along with the
Indian Pulsar Timing Array (InPTA) [12] and the MeerKAT Pulsar Timing Array
(MPTA) [13]. Combining data from older data sets from NANOGrav, the EPTA, and
the PPTA, the IPTA also found a consistent CURN in their second data release
[14]. Such a signal might arise from some currently unknown noise processes
shared by the measured pulsars [10, 15, 16]. Thus, to claim a detection of GWs
we require evidence of the telltale correlations between pulsar pairs, known
as Hellings and Downs (HD) correlations [17].

Among the 67 pulsars used in the GWB analysis in the 15-year data [18],
NANOGrav reports a stronger detection of the CURN process. Additionally,
NANOGrav reports evidence for HD correlations suggesting that the common
signal seen by the three PTAs is indeed a gravitational-wave background [19].
While the parameter estimation and subsequent astrophysical interpretation
remains somewhat dependent on noise models (see, e.g., [20] for more
information about the noise models used), the GWB remains consistent with an
origin of an ensemble of supermassive black hole binaries [21]. In addition to
searching for a GWB, other analyses have been performed looking for single
continuous-wave sources [22], anisotropy in the GWB [23], and possible hints
of new physics [24]. Future joint analyses in the IPTA will work toward
combining data to improve sensitivity to GW sources and their corresponding
parameters for astrophysical interpretation.

In performing analyses for HD correlations, PTA collaborations employ a
variety of techniques to process and analyze their data. NANOGrav uses a
modular pipeline that starts with radio telescope data and ends with
astrophysical GWB inference. Radio telescope data are processed into TOAs,
analyzed for outliers, and fitted with a timing model [18]. Subtracting the
predicted TOAs from the observed TOAs results in the TOA _residuals_. These
residuals are broken down into terms associated with deterministic signals, as
one might expect from an individual binary, as well as stochastic signals such
as those from a GWB. Stochastic signals may be further broken into different
types of noise processes: white and red. White noise has a flat frequency
spectrum and is associated with noise in telescope observations. Red noise has
a larger amplitude at lower frequencies and originates from pulsars themselves
in the form of spin noise [25], fluctuations in the dispersion measure (DM)
caused by the interstellar medium in between us and the pulsar, and a
potential GWB [26]. A single-pulsar white and red noise analysis is performed
at the end of this part in the pipeline. The TOAs, timing model, residuals,
and noise analyses comprise the initial input for any GW analysis.

Unlike other regimes of GW data analysis, the uneven sampling and large gaps
in NANOGrav data make directly transforming to the frequency domain difficult.
Most analyses are therefore performed in the time domain [27, 28]. In this
paper, we specialize to the case of a GWB and describe the analysis
implementations as they currently exist. The Bayesian GWB analysis started
with the pioneering work of [28]. While the brute force inversion of a full
$N_{\text{TOA}}\times N_{\text{TOA}}$ matrix was possible with $\sim 10^{3}$
TOAs at the time of publication, it would not work today with $\sim 5\times
10^{5}$ TOAs. The modern PTA likelihood was introduced in [29] where the
authors expanded the red noise in a set of Fourier coefficients. Lentati  _et
al._ [29] marginalized over the timing model as in [28] and added an
additional marginalization over the Fourier coefficients. van Haasteren and
Vallisneri [30, 31] used the connection between these methods and the theory
of Gaussian processes leading to further optimizations in the likelihood
computation and sampling.

Complementary to Bayesian approaches, the NANOGrav frequentist GWB analysis
uses the so-called optimal statistic (OS), an estimator of the amplitude and
significance of a GWB. Initially, this statistic was formulated in the
frequency domain [32], and later it was re-implemented in the time domain
[27]. Traditionally, it has been formulated in the regime where the amplitude
of the GWB is much lower than the amplitude of the intrinsic red noise in the
pulsars, an assumption that has been relaxed recently [33]. Additionally, the
optimal statistic relies upon an estimate of the intrinsic red noise in each
pulsar. Our estimates of intrinsic red noise are uncertain, and using a
specific choice of the red noise can result in a biased statistic [34]. The
solution is to create a hybrid Bayesian-frequentist “noise-marginalized”
optimal statistic, in which the optimal statistic is computed over many
posterior draws for our red noise model, creating a distribution for the
amplitude estimate of the GWB and its associated significance. It is common
practice to average the optimal statistic over this distribution [34],
although recently alternative approaches have also been proposed [35].

The traditional optimal statistic assumes that there is only one type of
spatial correlation in the data. Monopolar and dipolar correlations are also
possible, resulting from systematic issues such as clock errors [36] or solar
system ephemeris errors [37] respectively, and so potentially there will be
multiple spatially correlated signals in the data simultaneously. When
searching for only one type of spatially correlated signal at a time, one
might make spurious detections due to the contribution of other spatially
correlated signals. The solution to this is to simultaneously fit for multiple
spatial correlation patterns using the multiple-component optimal statistic
(MCOS) [38].

These methods have evolved conceptually over the past decade, and so have
their software implementations. NANOGrav maintains a set of publicly available
software packages that are written in the Python programming language and make
extensive use of numpy and scipy [39, 40]. The packages that we focus on in
this work include Enterprise , enterprise_extensions , PTMCMCSampler , and
la_forge [41, 42, 43, 44]. In this paper, we review these packages and perform
a campaign of simulations to validate them for the specific case of searching
for a GWB with NANOGrav data, as performed in NANOGrav’s 15-year GWB analysis
[19]. We base all simulations in this paper on the NANOGrav 15-year data set,
using all 67 pulsars that have been timed for more than three years [18].

Inspired by probabilistic programming packages such as Stan [45] and PyMC
[46], Enterprise allows the specification of probabilistic data models for
PTAs, and the evaluation of the resulting priors and likelihoods. By contrast,
enterprise_extensions contains pre-built models that are commonly used in PTA
analyses. It also includes _hypermodels_ used in Bayesian model selection by
way of product-space sampling, as well as implementations of the OS, MCOS, and
a noise-averaged version of these known as the “noise-marginalized” optimal
statistic (NMOS). We use PTMCMCSampler, a parallel-tempering enabled Markov-
chain Monte Carlo (MCMC) sampler, to approximate the posterior of the models
created with Enterprise. Finally, once we finish sampling, we use la_forge to
compress and post-process the resulting chains.

In Section II, we describe Enterprise, including the structure of the
software, the traditional likelihood computation, and a new faster
implementation for situations where the white-noise parameters remain
constant. We also describe the computational scaling associated with each
calculation. In Section III we discuss using Enterprise along with
PTMCMCSampler to explore the high dimensional parameter spaces in GWB
analyses. Frequentist methods implemented in enterprise_extensions are
discussed in Section IV. Next, we present tests on the Bayesian methods in
Section V and frequentist methods in Section VI. Finally, we discuss possible
future directions and conclude in Section VII.

##  II Enterprise

The Enhanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE
(Enterprise) [41] is a Python package built to analyze pulsar noise and timing
models, and to search for GWs in PTA data. It grew out of previous PTA data
analysis packages such as NX01 [47], PAL2 [48], and piccard [49]. In the
context of this paper, we are using Enterprise to search for an isotropic GWB
in simulations of the NANOGrav 15-year data set. Development of Enterprise
began in the interest of making a suite of tools that could analyze data from
any PTA consortium, and support international collaboration without requiring
significant knowledge of programming. As such, Enterprise includes the
flexibility to work with many different models, and even create new models of
noise and GW signals. We will first discuss the structure of the likelihood
implemented in Enterprise, and then describe its representation in the code.

###  II.1 The PTA Likelihood

Our TOAs $\bm{t}$ can be written as

$\bm{t}=\bm{t}_{\text{det}}+\bm{t}_{\text{stoc}}\,,$ (1)

where $\bm{t}_{\text{det}}$ is the deterministic part of the TOAs,
$\bm{t}_{\text{stoc}}$ is the stochastic part of the TOAs. After fitting the
deterministic part of the signal with a least squares fit,
$\bm{t}_{\text{M}}$,

$\bm{t}_{\text{det}}\approx\bm{t}_{\text{M}}+\bm{M}\bm{\epsilon}\,,$ (2)

where $\bm{M}\bm{\epsilon}$ is a Taylor expansion of the residual
deterministic part where $\bm{M}$ consists of derivatives in the Taylor
expansion. The details of the process used to produce the timing model in the
NANOGrav 15-year data set can be found in the data set paper [18]. The
stochastic part of the TOAs

$\bm{t}_{\text{stoc}}=\bm{Fc}+\bm{n}\,,$ (3)

where $\bm{c}$ represent the Fourier coefficients, $\bm{F}$ represents a
discrete Fourier transform of the red noise processes in the data, and
$\bm{n}$ consists of white noise. Combining all of the above and subtracting
off the timing model fit, we find

$\bm{\delta
t}=\bm{t}-\bm{t}_{\text{M}}\approx\bm{M\epsilon}+\bm{Fc}+\bm{n}\,.$ (4)

As in [29], we expand red noise in a set of Fourier coefficients,

$\bm{Fc}=\sum_{j=1}^{N}\left[X_{j}\sin\left(2\pi
f_{j}t\right)+Y_{j}\cos\left(2\pi f_{j}t\right)\right]\,,$ (5)

where alternating $\bm{X},\bm{Y}$ make up $\bm{c}$, $\bm{F}$ contains
alternating columns of sine and cosine components, and $f_{i}=i/T$ with $T$
the observing time span of the entire data set (16.03 years in the 15-year
data set 111The data set has been so-named, because the pulsar with the
longest observation time span contains 15.8 years of data.). Red noise may
consist of components intrinsic to the pulsar, such as spin noise, or even a
gravitational-wave background. We limit the number of frequency bins $N$ used
based on the model of each red noise signal. For pulsar intrinsic red noise,
we limit ourselves to 30 frequency bins which is sufficient to capture the
high-frequency content of the data. This corresponds to frequencies from about
2 nHz to about 60 nHz. The number of frequencies used in the GWB analyses
depends on how we model it. For example, a two-parameter power law with
amplitude and spectral index in the 15-year data set shows that the red noise
dips below the white noise at around 14 frequency bins or around 28 nHz [19].
While we could include more frequency bins (50 were included in [50]), adding
more frequencies on this model would bias the common power law spectral index
and amplitude, unless a more advanced red noise model is used [51]. Using a
model that allows each frequency bin to vary independently is known as a “free
spectrum model”. For this model, bins are not affected significantly by
adjacent bins, and we use the same 30 frequency bins on the GWB that we use on
the intrinsic red noise. Therefore, including more frequency bins on common
signals has negligible influence on a free spectrum analysis, but this model
also includes an additional parameter per frequency bin.

Several types of white noise parameters are searched over that adjust the TOA
uncertainties found during template fitting, e.g., [52]. EQUAD $\bm{Q}$ is an
extra factor added in quadrature to the TOA uncertainties. EFAC $\bm{G}$
rescales the TOA uncertainties and EQUAD together. Finally, ECORR $\bm{J}$ is
an extra term that correlates different frequency bands within the same epoch,
a term which here means a single observation. A single observation (epoch)
could be a single day, or it could be a combination of multi-band data from
separate days combined. Separate epochs remain completely uncorrelated. ECORR
can account for pulse jitter [20]. Each of these white noise signals add one
dimension per pulsar per observing backend per frequency band used to acquire
data. The white noise covariance matrix can be written as

$\bm{N}=\left<\bm{n}\bm{n}^{T}\right>=\sum_{\mu}\left[G_{\mu}^{2}\left(\sigma_{i}^{2}+Q_{\mu}^{2}\right)\delta_{ij}+J_{\mu}^{2}\delta_{e(i)e(j)}\right]\,,$
(6)

where the $i$th TOA belongs to the backend $\mu$, $\delta_{ij}$ denotes the
Kronecker delta, and $\delta_{e(i)e(j)}$ denotes another Kronecker delta that
equals 1 only when the epochs are the same for both TOAs considered and 0
otherwise.

Subtracting the deterministic and stochastic models results in the residual,

$\bm{r}=\bm{\delta t}-\bm{M\epsilon}-\bm{Fc}\,.$ (7)

Under a multi-variate Gaussian assumption for the noise the full likelihood
can then be written as

$p(\bm{\delta
t}\mid\bm{c},\bm{\epsilon})=\frac{1}{\sqrt{\det(2\pi\bm{N})}}\exp\left({-\frac{1}{2}\bm{r}^{T}\bm{N}^{-1}\bm{r}}\right)\,.$
(8)

We then group the matrices and vectors into a more compact notation,

$\displaystyle\bm{T}=\begin{bmatrix}\bm{M}&\bm{F}\end{bmatrix}\,,$
$\displaystyle\bm{b}=\begin{bmatrix}\bm{\epsilon}\\\ \bm{c}\end{bmatrix}\,,$
(9)

and place a Gaussian prior on the Fourier coefficients with covariance

$\bm{B}(\bm{\eta})=\begin{bmatrix}\bm{E}&0\\\
0&\bm{\phi}(\bm{\eta})\end{bmatrix}=\begin{bmatrix}\bm{\infty}&0\\\
0&\bm{\phi}(\bm{\eta})\end{bmatrix}\,,$ (10)

where $\bm{\eta}$ contains the hyper-parameters, i.e., the parameters of the
prior such as amplitudes or spectral indices of the GWB or red noise power
laws. We let $\bm{E}$ be a diagonal matrix of very large values ($10^{40}$ by
default) in $\bm{B}$. This places an improper, almost-infinite variance
Gaussian prior on the timing model parameters. However, these parameters are
well determined by pulsar timing observations and thus likelihood-dominated.
Upon inversion, this choice marginalizes over the timing model uncertainties
[29, 30, 31].

The covariance matrix for the red noise coefficients,

$\bm{\phi}=\left<\bm{c}\bm{c}^{T}\right>\,,$ (11)

can be constructed via blocks. Each block contains

$N_{\text{freq}}=\max\left(N_{\text{IRN}},N_{\text{GWB}}\right)\,,$ (12)

frequencies where $N_{\text{IRN}}$ is the number of frequencies used on the
intrinsic red noise, and $N_{\text{GWB}}$ is the number of frequencies used on
a red noise process common among pulsars (correlated or not). If we denote the
block containing the frequencies of pulsars as $(a,b)$, and we further specify
each frequency as $(i,j)$, then

$\left[\phi\right]_{(ai)(bj)}=\left\langle
c_{ai}c_{bj}\right\rangle=\delta_{ij}\left(\delta_{ab}\varphi_{ai}+\Phi_{ab,i}\right),$
(13)

with

$\Phi_{ab,i}=\Gamma_{ab}\Phi_{i},$ (14)

where the overlap reduction function (ORF), $\Gamma_{ab}$, describes the
correlations between pulsar pairs. In the isotropic GWB analysis, the ORFs we
can search for include a no-correlation ORF for CURN, monopolar correlations
which could be caused by clock corrections [36], dipolar correlations which
could be caused by solar system ephemeris errors [37], and the HD correlations
which are characteristic of a GWB:

$\displaystyle\Gamma_{ab}^{\text{CURN}}$ $\displaystyle=\delta_{ab}\,,$ (15)
$\displaystyle\Gamma_{ab}^{\text{MON}}$ $\displaystyle=1\,,$ (16)
$\displaystyle\Gamma_{ab}^{\text{DIP}}$ $\displaystyle=\cos{\theta_{ab}}\,,$
(17) $\displaystyle\Gamma_{ab}^{\text{HD}}$
$\displaystyle=\frac{1}{2}\delta_{ab}+\frac{3}{2}x_{ab}\ln
x_{ab}-\frac{1}{4}x_{ab}+\frac{1}{2}\,,$ (18)

where $x_{ab}=\left(1-\cos{\xi_{ab}}\right)/2$ and $\xi_{ab}$ is the angle
between pulsars $a$ and $b$ on the sky. The terms $\rho$ and $\kappa_{a}$ are
related to the power spectral density, $S(f)$ of the time delay caused by
intrinsic red noise or a GWB respectively as

$\varphi_{ai},\Phi_{i}=S(f_{i})\Delta f=S(f)/T\,.$ (19)

A typical model choice for these spectra is a power law,

$\varphi_{ai},\Phi_{i}=\frac{A^{2}}{12\pi^{2}}\frac{1}{T}\left(\frac{f_{i}}{f_{\text{ref}}}\right)^{-\gamma}\text{yr}^{2}\,,$
(20)

where the GWB characteristic strain is,

$h_{c}=A\left(\frac{f}{f_{\text{ref}}}\right)^{\alpha}\,,$ (21)

with $\gamma=3-2\alpha$. Provided that the GWB is made up of signals from an
ensemble of supermassive black hole binaries, we expect $\gamma=13/3$
($\alpha=-2/3$) [53]. The reference frequency, $f_{\text{ref}}$, traditionally
has been set to $1/\text{yr}^{-1}$, and this is the value that we use in our
simulations. For the intrinsic red noise and common uncorrelated red process,
$A$ and $\gamma$, $A$ and $\gamma$ are the hyperparameters, $\bm{\eta}$, for a
power law spectrum model.

In the traditional form of the likelihood, we analytically marginalize over
$\bm{b}$. This is a simultaneous marginalization over the timing model
uncertainty and red noise coefficients,

$p(\bm{\delta t}\mid\bm{\eta})=\int p(\bm{\delta
t}\mid\bm{b})p(\bm{b}\mid\bm{\eta})d\bm{b}\,,$ (22)

where

$p(\bm{b}\mid\bm{\eta})=\frac{\exp\left(-\frac{1}{2}\bm{b}^{T}\bm{B}^{-1}\bm{b}\right)}{\sqrt{\det(2\pi\bm{B})}}\,,$
(23)

is a Gaussian prior with hyper-parameters $\bm{\eta}$. Evaluating the integral
gives

$p(\bm{\delta
t}\mid\bm{\eta})=\frac{1}{\sqrt{\det(2\pi\bm{C})}}\exp\left(-\frac{1}{2}\bm{\delta
t}^{\mathrm{T}}\bm{C}^{-1}\bm{\delta t}\right)\,,$ (24)

where the covariance matrix is now

$\bm{C}=\bm{N}+\bm{TBT}^{T}\,.$ (25)

Such an inversion can be efficiently evaluated with the Woodbury matrix
identity,

$\bm{C}^{-1}=\bm{N}^{-1}-\bm{N}^{-1}\bm{T}\bm{\Sigma}\bm{T}^{T}\bm{N}^{-1}\,,$
(26)

with

$\bm{\Sigma}=\left(\bm{B}^{-1}+\bm{T}^{T}\bm{N}^{-1}\bm{T}\right)^{-1}\,.$
(27)

This formulation reduces the number of operations required for the
computational bottleneck from an inversion of the full covariance matrix
$\bm{C}$ which takes
$\mathcal{O}\left(N_{\text{p}}^{3}N_{\text{TOA}}^{3}\right)$ operations to the
inversion of $\bm{\Sigma}$ which takes
$\mathcal{O}\left(N_{\text{p}}^{3}(2N_{\text{freq}}+N_{\text{M}})^{3}\right)$.
Given that the number of TOAs is $\sim 10^{5}$, $N_{\text{freq}}=30$,
$N_{\text{M}}\sim 10^{2}$, the savings are significant.

The white noise covariance matrix is block-diagonal with a block for each
observation epoch. We invert each block efficiently using the Sherman-Morrison
formula

$\bm{N}_{b}^{-1}=\bm{H}_{b}^{-1}-\frac{\bm{H}_{b}^{-1}\bm{u}\bm{u}^{\mathrm{T}}\bm{H}_{b}^{-1}}{J^{-2}+\bm{u}^{\mathrm{T}}\bm{N}_{b}^{-1}\bm{u}}\,,$
(28)

where $\bm{H}_{b}$ contains the diagonal elements of the white noise
covariance matrix, and

$\bm{u}^{T}=(1,1,\cdots,1),$ (29)

with length of the number of TOAs in the epoch. While this could be a non-
negligible part of the computation, we fix the white noise in the GWB
analysis, and therefore we do not consider the speed of evaluation here.

###  II.2 Faster Likelihood for GWB Analyses

White-noise parameters are varied in noise runs prior to performing any GWB
analyses. After these initial runs, all white noise parameters are set to
their median marginalized posterior values to reduce the total number of
parameters from $\sim 10^{3}$ to $\sim 10^{2}$. When modeling the stochastic
processes as a power law, this results in 136 parameters: two parameters for
each of the 67 pulsar’s intrinsic red noise parameters, and two for the GWB.
White noise parameters in Enterprise are cached to speed up evaluation time
significantly after the initial likelihood evaluation.

Starting with the same multi-variate Gaussian in Equation 8, we now
marginalize in two steps instead of simultaneously. Marginalizing over the
timing model,

$p(\bm{\delta t}\mid\bm{c},\bm{E})=\int p(\bm{\delta
t}\mid\bm{\epsilon},\bm{c})p(\bm{\epsilon}\mid\bm{\eta})d\bm{\epsilon}\,,$
(30)

with the Gaussian prior

$p(\bm{\epsilon}\mid\bm{\eta})=\frac{\exp\left(-\frac{1}{2}\bm{\epsilon}^{T}\bm{E}^{-1}\bm{\epsilon}\right)}{\sqrt{\operatorname{det}(2\pi\bm{E})}}\,,$
(31)

and $\bm{E}=\left<\bm{\epsilon}\bm{\epsilon}^{T}\right>$. This integral
results in

$p(\bm{\delta t}\mid\bm{c})=\frac{\exp{\left(-\frac{1}{2}(\bm{\delta
t}-\bm{F}\bm{c})^{T}\bm{D}^{-1}(\bm{\delta
t}-\bm{F}\bm{c})\right)}}{\sqrt{\det(2\pi\bm{D})}}\,,$ (32)

with

$\bm{D}=\bm{N}+\bm{M}\bm{E}\bm{M}^{T}\,.$ (33)

We use the Woodbury matrix identity for this inversion,

$\bm{D}^{-1}=\bm{N}^{-1}-\bm{N}^{-1}\bm{M}\bm{\Lambda}^{-1}\bm{M}^{T}\bm{N}^{-1}\,,$
(34)

with

$\bm{\Lambda}=\bm{E}^{-1}+\bm{M}^{T}\bm{N}^{-1}\bm{M}=\bm{M}^{T}\bm{N}^{-1}\bm{M}\,,$
(35)

because $\operatorname{diag}(\bm{E})\rightarrow\infty$222Formally, we use the
matrix determinant lemma to evaluate the determinant in Equation 32 as
$\det(\bm{D})=\det\left(\bm{E}^{-1}+\bm{M}^{T}\bm{N}^{-1}\bm{M}\right)\det(\bm{E})\det(\bm{N}),$
(36) where $\operatorname{diag}(\bm{E})=\infty$. However, in practice we can
use $10^{40}$ as an effectively infinite value and avoid dealing with the
infinite determinant term.. Now, we complete the two-step marginalization
procedure by marginalizing over the Fourier coefficients,

$p(\bm{\delta t}\mid\bm{\eta})=\int p(\bm{\delta
t}\mid\bm{c})p(\bm{c}\mid\bm{\eta})d\bm{c}\,,$ (37)

with the Gaussian prior

$p(\bm{c}\mid\bm{\eta})=\frac{\exp\left(-\frac{1}{2}\bm{c}^{T}\bm{\phi}^{-1}\bm{c}\right)}{\sqrt{\operatorname{det}(2\pi\bm{\phi})}}\,,$
(38)

resulting in

$p(\bm{\delta
t}\mid\bm{\eta})=\frac{1}{\sqrt{\det(2\pi\bm{K})}}\exp\left(-\frac{1}{2}\bm{\delta
t}^{T}\bm{K}^{-1}\bm{\delta t}\right)\,,$ (39)

where

$\bm{K}=\bm{D}+\bm{F}\bm{\phi}\bm{F}^{T},$ (40)

and once again we use the Woodbury matrix identity to invert this covariance
matrix. We find

$\bm{K}^{-1}=\bm{D}^{-1}-\bm{D}^{-1}\bm{F}\bm{\Theta}\bm{F}^{T}\bm{D}^{-1}\,,$
(41)

with

$\bm{\Theta}=\left(\bm{\phi}^{-1}+\bm{F}^{T}\bm{D}^{-1}\bm{F}\right)^{-1}\,.$
(42)

As long as $\bm{D}^{-1}$ can be cached in its entirety, the computation will
be faster with the two-step marginalization, Equation 39, over the
simultaneous marginalization procedure, Equation 24. Now, the $\bm{\Theta}$
inversion in Equation 42 dominates the likelihood computation with
$\mathcal{O}\left(N_{p}^{3}(2N_{\text{freq}})^{3}\right)$ operations.

###  II.3 Empirical Computational Scaling for Likelihood Evaluations

(a) HD Model Comparison

(b) CURN Model Comparison

Figure 1: (a) Comparison between the empirical scaling of HD-correlated models
for both the simultaneous marginalization procedure (labeled “traditional
likelihood”) and the new two-step marginalization procedure (labeled “fast
likelihood”) which runs faster for models where the white noise is not being
varied. In the latter case, we can cache the $\bm{D}$ matrix resulting in a
drastic reduction in computation time at the cost of some memory. Inversions
in both methods use a sparse Cholesky decomposition method on a single CPU
core. This new marginalization procedure results in an average speed increase
per evaluation of $5.65\times$ for 67 pulsars. (b) Comparison between the CURN
models using the different marginalization procedures. Inverting the now
diagonal matrix of these models is trivial. The reduction in number of
elements required to be inverted results in a speed increase of $3.80\times$
for 67 pulsars. Individual points and uncertainties are found by averaging
over 100 evaluations of a model and taking the standard deviation. Fits to the
points have been made with a nonlinear least squares algorithm fitting $A,B,C$
to a function $f(N_{p})=AN_{p}^{B}+C$.

As shown in Figure 1, we find a significant improvement in the fast version of
the likelihood, Equation 39, over the traditional likelihood Equation 24. With
67 pulsars, in the HD-correlated model, the fast likelihood leads to 5.65
times faster evaluation than the traditional computation. The speed up factor
is reduced to 3.8 in the CURN model. Empirically, as the number of pulsars
increases, the difference in evaluation times also increases for the HD-
correlated case.

Empirically, we do not find that all computations currently scale as
$\mathcal{O}(N_{p}^{3})$. By fixing the number of frequency bins used in each
model and including pulsars in alphanumeric order, we can compare the matrix
inversion cost among models and check how the computations scale with number
of pulsars. There are two separate likelihood implementations: a sparse
version that uses scipy.sparse.csc compressed sparse column matrices along
with a sparse Cholesky decomposition scikit-sparse [54, 55, 56] and a version
that uses a dense Cholesky decomposition. Through profiling, we find that the
dense computation indeed scales as $\mathcal{O}(N_{p}^{3})$ and is generally
slower than the sparse method on a single core. The sparse method instead
scales as $\mathcal{O}(N_{p}^{2})$ using the fast likelihood, Equation 39, or
as $\mathcal{O}(N_{p}^{2.3})$ when using the slower version of the likelihood,
Equation 24.

This $\mathcal{O}(N_{p}^{2.3})$ empirical scaling of the traditional
likelihood displayed in 1(a) can be traced to the sparse Cholesky inversion
$\bm{\Sigma}^{-1}$. In general, the sparse Cholesky inversion does not have a
strict scaling, but depends on how the matrix is structured. Because of this,
the scaling of the sparse Cholesky inversion may have some dependence on the
number of frequencies and pulsars used. By profiling the code to see which
parts use the largest fraction of the total evaluation time, we find that in
the sparse cases, changing from a dense to a sparse matrix takes up the
majority of the computation time, and the inversion is performed very quickly
by comparison. This means that our inversion has been sped up to a point where
it is no longer the bottleneck of the computation, and the
$\mathcal{O}(N_{p}^{2})$ pieces of the computation dominate the overall
scaling for Equation 39 evaluations. Future work will aim to reduce the
computation time required in these $\mathcal{O}(N_{p}^{2})$ operations,
perhaps by working with sparse matrices from the start.

###  II.4 Structure of Enterprise

Enterprise is a thoroughly object-oriented package that defines a PTA data
model as a hierarchical structure (see Figure 2). The top-level object is PTA,
which interfaces with the sampler by way of the get_lnlikelihood and
get_lnprior methods: these evaluate the (log) Gaussian-process-marginalized
likelihood of Equation 24 and the total prior respectively, taking as input a
Python dictionary of parameter values, or alternatively a numpy vector with
the parameters in the same order as in the property PTA.params.

Figure 2: The hierarchical structure of Enterprise.

The PTA object is created from a sequence of SignalCollection objects, each of
which corresponds to a pulsar in the array, and represents the pulsar’s
complete data model. SignalCollection implements a set of methods that return
the vector and matrix constituents used by PTA.get_lnlikelihood: for instance,
get_residuals for the residuals $\bm{\delta t}$, get_delay for any
deterministic delays, get_ndiag for the white-noise matrix $\bm{N}$, get_basis
and get_phi for Gaussian-process bases $\bm{T}$ and prior matrix $\bm{B}$. All
these methods take a Python dictionary of parameter values.

Each SignalCollection object consists of one Pulsar object and any number of
Signal objects. Pulsar uses the PINT or libstempo packages to read pulsar data
from par and tim files, standard formats for storing timing model parameters
and timing data (see, e.g., [57]), and it provides the vectors residuals (for
$\bm{\delta t}$), toaerrs (for $\bm{\sigma}$), and freqs (for the pulse radio
frequencies), as well as the timing-model design matrix Mmat (i.e., $\bm{M}$),
and several more pulsar properties. The Signal objects, which come in a
variety of subtypes, implement deterministic signals and noise components, and
provide the building blocks that SignalCollection and (downstream) PTA
assemble into a full likelihood. Specifically, deterministic signals define
get_delay, white-noise components define get_ndiag, and Gaussian processes
define get_basis and get_phi. _Common_ Gaussian processes, used most notably
for the HD-correlated background, have parallel bases for each pulsar (with
the same $n_{\mathrm{GP}}$ but different $n_{\mathrm{obs}}$), provide a
get_phicross method in addition to get_phi to fill off-diagonal prior-matrix
elements, and are meant to be used with common hyperparameters.

The package includes a number of optimizations for speed, memory, and ease of
configuration:

  * •

All objects keep track of which parameters affect the output of their methods,
and cache results if those parameters are not changed, or are set as
constants. This is useful, for instance, in analyses that freeze the white-
noise parameters, or in stochastic schemes that update parameters in blocks.
The SignalCollection and PTA objects have methods for intermediate matrix
combinations such as $\bm{T}^{T}\bm{N}^{-1}\bm{T}$, and these are also cached
on the basis of the relevant parameters.

  * •

The SignalCollection object combines delays from all deterministic signals
into a single vector, assembles white-noise matrices into a single matrix, and
stacks Gaussian-process bases and priors into two combined matrices. The
object has the ability of reusing basis vectors that appear identically in
multiple processes, such as Fourier vectors for red noise and the GWB.

  * •

The get_ndiag methods return custom “kernel” objects that represent the
constituents of the measurement-noise matrix $\bm{N}$. Whether the kernels are
represented internally as a vector for EFAC/EQUAD noise, a Sherman-Morrison
decomposition for ECORR [see Equation 28], or a sparse matrix for even more
general cases, they all provide “solve” methods that return combinations such
as $\bm{T}^{T}\bm{N}^{-1}\bm{y}$ and $\bm{T}^{T}\bm{N}^{-1}\bm{T}$ without
actually computing and storing $\bm{N}^{-1}$. Furthermore, the kernel objects
know how to combine themselves with other kernel objects, creating an
optimized object used in the likelihood calculation.

  * •

The PTA object defines an optimized get_phiinv method that can choose between
multiple inversion strategies to build the global $\bm{B}^{-1}$ depending on
the structure of $\bm{B}$.

Enterprise is configured by building a SignalCollection _template_ from a
sequence of Signal templates. Enterprise includes “factories” that build
templates for commonly used noise components such as MeasurementNoise (for
radiometer noise), EcorrKernelNoise (for jitter-like noise), TimingModel (for
the Gaussian process with improper prior used to marginalize over timing-model
corrections), FourierBasisGP (for Gaussian processes with a sine/cosine
Fourier basis), and FourierBasisCommonGP (for Fourier Gaussian processes with
correlations across pulsars), and more. Each template must be assigned one or
more Parameter objects that encode the priors chosen for the relevant
parameters. For instance, MeasurementNoise may be passed efac = Normal(1,
0.25) and log10_t2equad = Uniform(-8.5, -5) to indicate that in Equation 6
$F\sim\mathcal{N}(1,0.25)$ and $\log_{10}Q\sim U(-8,-5)$. Parameters may also
be passed as Constant if they will not be sampled stochastically.

The SignalCollection template is then applied to each Pulsar object in turn,
generating _instantiated_ SignalCollection and Signal objects, in which
parameters are specialized to the pulsar (e.g., efac becomes B1855+09_efac)
and the Pulsar quantities are made available locally. For most Signals, a
further kind of specialization is possible in which a _selection_ function is
passed to the Signal, effectively splitting it into multiple copies, each
applied with different parameters to a different section of the data. For
instance, the selection by_backend, when given to MeasurementNoise, would
split EFAC parameters for the NANOGrav B1855+09 pulsar into
B1855+09_430_ASP_efac, B1855+09_430_PUPPI_efac, B1855+09_L-wide_ASP_efac,
B1855+09_L-wide_PUPPI_efac, each applied to the subset of the residuals
obtained with that receiver backend.

Configuration is completed by collecting the set of instantiated
SignalCollection into a single PTA object. See the NANOGrav 12.5-yr and 15-yr
tutorials for examples of creating a standard PTA object that is ready to
compute log likelihoods and priors.

Enterprise includes a number of additional facilities and utilities that
simplify the task of creating a model. For instance, the prior for a Fourier
Gaussian processes can be written simply as Python functions that take a
vector of frequencies and return the diagonal components of $\bm{\phi}$. By
wrapping the prior in parameter.function, the arguments of the Python function
are automatically book-kept as model parameters and specialized to pulsar and
selections. Furthermore, if the Python function includes arguments that are
defined in the Pulsar object (such as residuals or toaerrs), these are passed
to the function automatically once the Signal that uses the function has been
specified.

Of course, Enterprise already defines a number of commonly used Gaussian-
process priors and ORFs. It also implements the deterministic model of solar-
system-ephemeris uncertainties [37] used for older NANOGrav data releases
(utils.FourierBasisCommonGP_physicalephem); it can draw random values of model
parameters according to their prior (parameter.sample); it can sample the
conditional distributions of Gaussian-process weights given their
hyperparameters (utils.ConditionalGP); it can simulate full realizations of
the data model (simulate); it can handle the faster likelihood of subsection
II.2 with a special noise kernel MarginalizingTimingModel that holds the
components of $\bm{D}$ internally and performs the Woodbury inversion of
Equation 34 transparently. See the Enterprise documentation for this and much
more.

##  III Bayesian Methods and Sampling the PTA Parameter Space

GW data analyses over the past decade have made use of both Bayesian and
frequentist statistical techniques. Bayesian methods rely on Bayes’ theorem
for parameter estimation and model selection,

$p(\bm{\theta}\mid\bm{y})=\frac{p(\bm{\theta})p(\bm{y}\mid\bm{\theta})}{p(\bm{y})}\,.$
(43)

Given a descriptive model of the relevant data as a likelihood
$p(\bm{y}\mid\bm{\theta})$ and a prior distribution $p(\bm{\theta})$ on the
parameters, we can sample the (unnormalized) posterior
$p(\bm{\theta}\mid\bm{y})$, thus reallocating credibility from the prior
across parameter values. The normalizing factor in the denominator is known as
the marginal likelihood or evidence, and it may be used in deciding which of a
set of models is preferred.

Numerical methods to estimate $p(\bm{\theta}\mid\bm{y})$ when the dimension of
$\bm{\theta}$ is large typically depend on a form of stochastic sampling to
perform high-dimensional integrals for both parameter estimation and model
selection. In NANOGrav, we use PTMCMCSampler [43], a parallel-tempering (PT)
enabled Metropolis-Hastings MCMC sampler. This has been the NANOGrav sampler
of choice for many years. It has also recently been compared to other options
and recommended for use in PTA data analysis [58]. To sample the high-
dimensional PTA parameter space, we use the default proposal distributions
that come with PTMCMCSampler, and some other custom proposals which will be
enumerated and described below. By default, PTMCMCSampler uses single-
component adaptive metropolis (SCAM), adaptive metropolis (AM), and
differential evolution (DE) proposal distributions, referred to in the code as
“jump proposals.” Additionally, the sampler supports adding custom
distributions. Details of this sampler were also discussed in [59], but we
reiterate them here with any changes that have since been made.

###  III.1 Metropolis-Hastings Algorithm

We use the Metropolis-Hastings algorithm [60, 61] to stochastically sample the
posterior distribution given by the likelihood and chosen priors. After enough
iterations, the sampler converges to a stable distribution that approximates
the posterior well regardless of where in the parameter space we start. We
will show tests of this in Section V. Under certain circumstances convergence
is guaranteed with an infinite number of samples, but in the high-dimensional
parameter space that we search over, the amount of time required for
convergence can be prohibitive without well-chosen proposal distributions.

First, we start with an initial set of parameters $\theta$ at iteration $t=0$.
Then, we draw a new point $\bm{\theta}_{*}$ from a proposal distribution for
each iteration $J_{t}\left(\bm{\theta}_{*},\bm{\theta}\right)$. In the
Metropolis algorithm [60], such proposal distributions are required to be
symmetric, but the Metropolis-Hastings (MH) algorithm [61] allows for
asymmetric distributions by the inclusion of the Hastings ratio. The
Metropolis-Hastings algorithm leads to the proposed point being accepted with
log probability [62]

$\ln A=\min\left(0,\ln R+\ln H\right)\,,$ (44)

where $R$ is the ratio of probabilities between the proposed point and the old
one,

$\ln R=\ln p(\bm{\theta}_{*}\mid\bm{y})-\ln p(\bm{\theta}\mid\bm{y})\,,$ (45)

and $H$ is a ratio which accounts for asymmetric proposal distributions

$\ln H=\ln J_{t}(\bm{\theta}\mid\bm{\theta}_{*})-\ln
J_{t}(\bm{\theta}_{*}\mid\bm{\theta})\,.$ (46)

Iterating for each $t$ returns a chain of samples from the posterior
distribution. However, the number of samples required to settle into a
stationary distribution that approximates the posterior well depends on
several factors including the auto-correlation length, which is related to the
number of parameters, and whether the distribution is multi-modal where the
chain might get stuck in a single mode. At each iteration, we draw a point
from the proposal distributions which is either accepted or rejected. If the
new point is accepted, it is added to the chain. Otherwise, if the new point
is rejected, the current point is added to the chain instead. Chains produced
in this way contain samples which are not completely independent, and they are
correlated with themselves. In the case of a MH sampler, the auto-correlation
length is typically of the same order as the number of parameters.

###  III.2 Parallel Tempering

To improve exploration and mixing of the sampler, we sample multiple chains
with different exponents, known as temperatures, and propose swaps between
them in a sampling scheme known as PT [63]. The posterior now contains the
likelihood raised to some power,

$p\left(\bm{\theta}\mid\bm{y},\beta\right)=p\left(\bm{\theta}\right)p\left(\bm{y}\mid\bm{\theta}\right)^{\beta}\,,$
(47)

where $\beta=1/T$ and $T$ is known as the chain’s temperature. Samples from
one chain are propagated to the next via swap proposals between chains of
different temperatures. The higher the temperature is, the more the posterior
becomes like the prior, enabling exploration of the parameter space and
reducing the auto-correlation length of all chains via swaps, increasing the
number of effectively independent samples. Though the cost of this scheme is
more evaluations of the likelihood, we often find that this is more efficient
because of the increased number of effectively independent samples returned.
It can also be used to find the evidence as will be discussed in a later
section. We set a temperature for each of the chains using a geometrically
spaced ladder,

$T_{i}=\left(\frac{T_{\text{max}}}{T_{\text{min}}}\right)\exp(i)\,,$ (48)

where $i=0,1,2,\ldots$, which should result in a $\sim 25\%$ temperature swap
acceptance rate between adjacent chains when sampling a multi-variate Gaussian
distribution [59]. In PTMCMCSampler, swaps are proposed between adjacent
chains, though in general swaps could be proposed between any two chains in
the ladder. Swaps are accepted between temperatures $T_{i}$ and $T_{j}$ with
log probability,

$\ln A_{ij}=\min\left(0,\left(\beta_{j}-\beta_{i}\right)\ln L_{ij}\right)\,,$
(49)

where

$\ln L_{ij}=\ln p(\bm{y}\mid\bm{\theta}_{i})-\ln
p(\bm{y}\mid\bm{\theta}_{j})\,,$ (50)

is the log likelihood ratio. Other temperature ladders are possible including
adaptive temperature spacing based on a constant acceptance rate [64]. This
improves the PT scheme when sampling a posterior which is not a multi-variate
Gaussian. While these different temperature spacings certainly have
advantages, they are not currently implemented in PTMCMCSampler.

To perform parallel tempering swaps, we use multiple cores to sample each
chain simultaneously through the use of MPI [65] and MPI4Py [66, 67].
Previously, the sampler used an asynchronous model for the temperature swaps,
so that chains could sample at their own paces and swap as soon as the next
one down reaches a specified interval. These processes are now synchronized
using blocking commands, which are necessary for the standard product space
sampling method that NANOGrav uses for model selection [68] (see Appendix A
for a discussion of the necessity of synchronizing the sampler).

###  III.3 Jump Proposals

The GWB analysis uses several proposal distributions. These distributions are
critical for timely convergence and exploration of the parameter space.
Ideally, jump proposals match the posterior closely to minimize the auto-
correlation length of the chain and thus reduce the number of samples that
need to be taken. The combination of all of these proposals has proven to work
well for the problem at hand, even with the $\mathcal{O}(100)$ parameters that
we work with in the 15-year data set. The proposal distributions discussed
here consist of the default in PTMCMCSampler along with empirical
distributions and prior draws.

####  III.3.1 Adaptive Metropolis

Upon initializing the sampler, PTMCMCSampler takes a list of “parameter
groups” as an argument. If we believe that multiple parameters will be
correlated, we can add them as a group, and the sampler will propose jumps in
this subspace. The full sample space of all parameters together is always a
group regardless of new groups. However, if no groups are given, the entire
sample space is considered the only group. Typically, power law amplitude and
spectral indices are grouped together with a sampling group for each
individual pulsar.

PTMCMCSampler also requires a sample covariance matrix at initialization.
Periodically, we compute a sample covariance matrix $\bm{C_{s}}$ from the
chain’s history for each sample group using an online algorithm (see Appendix
C) to avoid storing the entirety of the chain. $\bm{C_{s}}$ is then decomposed
via a singular-value decomposition (SVD),

$\bm{C_{s}}=\bm{U_{s}}\bm{\Sigma_{s}}\bm{V_{s}}^{T}\,.$ (51)

This provides a robust generalization of the eigenvalue decomposition and is
used to jump along correlated directions in the parameter space.

The adaptive Metropolis proposal distribution [69] uses the sample covariance
matrix to propose jumps in uncorrelated directions of the parameter space.
First, we move the usual parameters at the $i$th step of the chain
$\bm{\theta}_{i}$ into the parameter combinations using

$\bm{\zeta}_{i}=\bm{U_{s}}^{T}\bm{\theta}_{i}\,.$ (52)

Each jump proposed is given by a normal distribution,

$\bm{\zeta}_{i+1}=\bm{\zeta}_{i}+\sqrt{\bm{\Sigma_{g}}}\,\mathcal{N}(0,c_{d})\,,$
(53)

where $\bm{\Sigma_{g}}$ is the sample covariance matrix for the specific group
in which the jump is proposed, and

$c_{d}=\frac{2.4s}{\sqrt{2n_{\text{dim}}}}\,,$ (54)

where $s$ is a scale parameter and $n_{\text{dim}}$ depends on the number of
dimensions of the group that the jump is proposed in.. For each default jump,
$3\%$ of jumps have their scale multiplied by 10, $7\%$ of jumps have their
scale multiplied by 0.2, and the other $90\%$ have unmodified scale. In all
cases, the relative scale of the jump is adjusted based on the temperature of
the chain. The mix of small, medium, and large jumps helps the sampler find
the scale of the parameter space being explored. To project back into a jump
in the original parameters, we use

$\bm{\theta}_{i+1}=\bm{U}\bm{\zeta}_{i+1}\,.$ (55)

####  III.3.2 Single-Component Adaptive Metropolis

Similar to the adaptive Metropolis jump, the single-component adaptive
Metropolis (SCAM) jump [70] uses the sample covariance matrix, but only moves
along one uncorrelated direction in the parameter space. Once again, we start
by projecting onto the uncorrelated combinations of parameters
$\bm{\zeta}=\bm{U}^{T}\bm{\theta}$. We then propose a jump in a single
parameter direction as,

$\zeta_{i+1}^{j}=\zeta_{i}^{j}+\sigma_{s}^{j}\mathcal{N}\left(0,c_{d}\right)\,,$
(56)

where $\sigma_{s}^{j}$ is the $j$th diagonal element of the diagonal matrix
$\sqrt{\bm{\Sigma_{s}}}$, and $j$ labels the uncorrelated parameter for which
we are proposing a new point. Finally, we move back to the original set of
parameters using

$\bm{\theta}_{i+1}=\bm{U}\bm{\zeta}_{i+1}\,.$ (57)

####  III.3.3 Differential Evolution

The final default proposal distribution in PTMCMCSampler uses a simple genetic
algorithm known as differential evolution (DE) [71]. This algorithm takes two
samples from the history of the chain, subtracts them, and proposes a jump
along that direction. In the current version of PTMCMCMSampler, the full chain
is not stored, and it instead draws from a buffer formed over many unthinned
iterations of the sampler. By keeping this buffer much longer than the auto-
correlation length of the sampler, we ensure that the draws come from a
stationary distribution. The DE jump draws two samples $\bm{\theta}_{m}$ and
$\bm{\theta}_{n}$ and then proposes a jump

$\bm{\theta}_{i+1}=\bm{\theta}_{i}+s_{DE}\left(\bm{\theta}_{m}-\bm{\theta}_{n}\right)\,,$
(58)

where $s_{DE}=1$ or

$s_{DE}\sim\text{Uniform}\left[0,\frac{2.4}{\sqrt{2\beta
n_{\text{dim}}}}\right]\,,$ (59)

each with $50\%$ probability with $\beta=1/T$ and $n_{\text{dim}}$ the number
of dimensions in which the jump is proposed.

####  III.3.4 Empirical Distributions

Before the GWB analysis, we run a noise analysis on each pulsar individually.
This run includes white noise (EFAC, EQUAD, ECORR) and a power-law intrinsic
red noise (amplitude and spectral index). From each of these runs, we find a
posterior for the intrinsic red noise, and we create white noise dictionaries
with which to set the white noise values constant during the full analyses.
Out of these posteriors, we can make 1D or 2D histograms that we can then
sample from during the full GWB analysis to propose as new points. During the
creation of these histograms, all bins’ counts are incremented by one to allow
exploration of the entire prior. They are then checked to make sure they cover
the prior before starting any sampling that includes them so that they do not
bias parameter recovery. Such histograms are known as empirical distributions
(see Appendix A in [72]). Typically, we use 2D empirical distributions on the
intrinsic red noise parameters for each pulsar. This provides excellent
proposals if the empirical distributions are somewhat close to the posteriors
on the parameters which we propose jumps in. Empirical proposal distributions
reduce the number of samples required to achieve stationarity, sometimes
called the burn-in, significantly. Importantly, empirical distributions are
only a small part of our overall mix of proposal distributions which include
many proposals that suggest points across the entire prior space.

###  III.4 Parameter Estimation

Parameter estimation provides crucial information for astrophysical inference.
The full multi-dimensional posterior, can be used to find the marginalized
posteriors for each of the parameters. Convergence and exploration of the
parameter space is critical to finding the maxima of the likelihood function
and to sampling them effectively. GWB analyses typically involve millions of
likelihood evaluations due to the large auto-correlation lengths of the
chains. It can take many evaluations to get a reasonable number of effective
samples. As a test of the procedure, we have verified the different runs on a
single simulation return the same result regardless of where we start in the
parameter space.

####  III.4.1 Gelman-Rubin $\hat{R}$ Diagnostic

We use the Gelman-Rubin $\hat{R}$ diagnostic [73] to check the stationarity of
the chain. This does not necessarily mean that the chains have converged, but
it tells us that the samples are coming from one part of the parameter space.
The diagnostic splits the MCMC chain into multiple segments and checks the
within-chain and between-chain statistics to confirm that the chain is in a
stationary state. A threshold is required to tell whether the chain passes or
fails the diagnostic. As suggested by [74], we use 1.01 as the $\hat{R}$
threshold. In all tests performed in this work, we use the $\hat{R}$
diagnostic to make sure the chains are stationary before performing tests.
This diagnostic is implemented in la_forge for ease of use in PTA data
analysis.

###  III.5 Model Selection: Bayes Factor Calculations

To compare models, we need to compute the marginal likelihood or evidences for
each model. Division of these evidences return Bayes factors. The standard
method of calculation used in NANOGrav is a form of product space sampling
[75, 68]. We refer to this method as the HyperModel (henceforth hypermodel)
framework, and it is implemented in enterprise_extensions. Other methods for
finding the model evidence include reweighting the posterior [76],
thermodynamic integration [77, 78], and nested sampling [79].

####  III.5.1 HyperModel Framework

The hypermodel concatenates different models into a single model by combining
their joint parameters into a set of parameters that contains only the unique
parameters between the models. During sampling, a continuous “switch”
parameter called nmodel changes between the models turning on only the
parameters that belong to the “on” model. By sampling the models and this
switch between them, we can compute an odds ratio comparing how many times
each of the models were sampled. This corresponds directly to a Bayes factor
between the two models. The uncertainties are then computed using a standard
bootstrap in which we resample the thinned nmodel marginalized posterior with
replacement and recompute the odds ratio. The odds ratio is averaged and a
standard deviation calculated over a number of realizations to give the final
Bayes factor 333Note that the odds ratio is only equal to the Bayes factor if
the prior odds are the same for each model in consideration. with
uncertainties.

In current data, we often find the situation where Bayes factors for one model
or another are significantly disfavored. Adding a log weight to the log
likelihood can remedy such a situation. To accomplish this, we add a constant
value to the log likelihood, scaling the likelihood by the exponential of the
log weight. We can find a weight estimate by subtracting the maximum posterior
values between the two models being compared. This results in a more even
mixing between the two models, and the weight can be undone in post-processing
by multiplying the Bayes factor by the exponential of the log weight.

####  III.5.2 Reweighting

Reweighting is a simple technique that utilizes existing samples from a
probability distribution, the approximate, to obtain an estimate of some other
probability distribution, the target, which shares support with the
approximate. Each existing sample is “weighted” by the ratio of the target and
approximate probability densities. These weighted samples are an estimate of
the target distribution whereas the weights can be used to estimate Bayes
Factors and uncertainties. Since the samples of the distribution have already
been produced, each weight can be calculated in parallel, increasing the speed
at which the target space can be evaluated. Reweighting results in a reduction
in the number of effectively independent samples based on how disjoint the
posteriors are between the approximate and target, but it often is still much
faster than directly sampling the target posterior directly. This technique is
particularly effective in cases where the two distributions have similar
support on similar regions in the parameter space but one distribution is
significantly faster to evaluate. In the context of PTAs, reweighting has been
used to generate HD posteriors from the faster-to-evaluate CURN posteriors
[76].

####  III.5.3 Nested Sampling

Nested sampling [80], a staple for GW analyses in current ground-based
detectors, returns a Bayes factor and posterior samples. It computes the
evidence integral by turning the multi-dimensional integral into a one-
dimensional integral. $N$ “live” points are sampled from the priors and the
space outside of the lowest likelihood point is removed, thereby reducing the
volume considered by approximately $1/N$. In this way, nested sampling climbs
the likelihood distribution in a global way, and eventually reaches a stopping
criterion set by the error on the log evidence. Overall it has a reputation
for being easy to use, good at finding multi-modality, and has stopping
conditions that do not require much input from the user [81]. In PTA data, we
find it difficult to use nested sampling on the full GWB analysis due to the
high number of parameters. However, we use nested sampling with a reduced set
of the data below to check our Bayes factors with $30$ parameters. In this
case, using Ultranest [79] returns the evidence for each model, but with much
larger uncertainty over the same computation time span as the hypermodel.

####  III.5.4 Thermodynamic Integration

If we use enough chains over a broad enough set of temperatures, parallel
tempering also can be used to compute the evidence. By taking the average of
the log likelihood on each sampled temperature, we can integrate over this to
yield the model evidence. The log evidence is given by [77, 78],

$\ln p(\bm{\delta t}\mid M)=\int_{-\infty}^{0}\beta\mathrm{E}_{\beta}[\ln
p(\bm{\delta t}\mid\bm{\eta},M)]d\ln\beta\,,$ (60)

where $\beta=1/T$, $M$ is the model of interest, and we are integrating with
respect to $\ln\beta$. We use two separate methods to evaluate uncertainties
on evidence estimates. In one, we use a cubic spline which is fit using a
trans-dimensional algorithm known as reversible-jump MCMC [82], and in the
other we use a bootstrap on an interpolation between points in the integrand
of Equation 60. There are two types of uncertainty here: the discretization
error which is determined by the number of temperatures that we choose, and
the sampling error which is determined by the number of independent samples.

##  IV Frequentist Methods and the Optimal Statistic

Fully Bayesian methods, especially when used for model selection, can be
computationally expensive. In this section, we consider a detection statistic
and an estimator for the GWB that are built from directly cross-correlating
arrival times between pulsars. We first consider the statistic in the case
where the amplitude of the GWB is small compared to the noise, which is what
has traditionally been assumed. This statistic is still important for null
hypothesis testing, although it can be improved upon when used as an estimator
for the strength of the background. We then move on to consider the situation
where the background cannot be neglected, and present an estimator that
properly accounts for the background itself. We also discuss how to construct
a “binned” estimator across the sky to yield a HD reconstruction that takes
into account the size of the background. We finish by presenting a version of
the optimal statistic (OS) that simultaneously fits for multiple correlation
patterns. We also highlight the effect of the choice of noise parameters used
in constructing the OS.

###  IV.1 Traditional optimal statistic

We begin by considering the noise-weighted match between the correlation of
data in pulsar $a$ with data in pulsar $b$:

$\displaystyle\rho_{ab}=$ $\displaystyle\frac{\bm{\delta
t}_{a}^{T}\bm{\mathbf{P}}_{a}^{-1}\tilde{\bm{\Phi}}_{ab}\mathbf{P}_{b}^{-1}\bm{\delta
t}_{b}}{\operatorname{tr}\left(\mathbf{P}_{a}^{-1}\tilde{\bm{\Phi}}_{ab}\mathbf{P}_{b}^{-1}\tilde{\bm{\Phi}}_{ba}\right)}\equiv\bm{\delta
t}_{a}^{T}\bm{Q}_{ab}\,\bm{\delta t}_{b}\,,$ (61) $\displaystyle\bm{Q}_{ab}=$
$\displaystyle\frac{\mathbf{P}_{a}^{-1}\tilde{\bm{\Phi}}_{ab}\mathbf{P}_{b}^{-1}}{\operatorname{tr}\left(\mathbf{P}_{a}^{-1}\tilde{\bm{\Phi}}_{ab}\mathbf{P}_{b}^{-1}\tilde{\bm{\Phi}}_{ba}\right)},$
(62)

where for two different pulsars, $a$ and $b$, we have

$\displaystyle\mathbf{P}_{a}$ $\displaystyle=\left<\bm{\delta t}_{a}\bm{\delta
t}_{a}^{T}\right>=\bm{D}_{a}+\bm{F}_{a}\bm{\phi}_{aa}\bm{F}_{a}^{T}\,,$ (63)
$\displaystyle\bm{\tilde{\Phi}}_{ab}$
$\displaystyle=\frac{\bm{F}_{a}\bm{\phi}_{ab}\
\bm{F}_{b}^{T}}{\Gamma_{ab}A^{2}_{gw}}\,.$ (64)

The normalization of $\bm{\tilde{\Phi}}_{ab}$ is chosen such that
$\langle\rho_{ab}\rangle=\Gamma_{ab}A^{2}_{\textrm{gw}}.$ In the small signal
regime, the variance of this correlation is
$\sigma_{ab}^{2}=\left(\operatorname{tr}\mathbf{P}_{a}^{-1}\tilde{\bm{\Phi}}_{ab}\mathbf{P}_{b}^{-1}\tilde{\bm{\Phi}}_{ba}\right)^{-1}.$
If we assume that the cross-correlation for one pair of pulsars is not
correlated with the cross-correlation of another pair of pulsars, i.e.,
$\langle\rho_{ab}\rho_{cd}\rangle\propto\delta_{ac}\delta_{bd}$, then we can
perform a variance-weighted, HD-matched sum over these correlations to
estimate the amplitude of the GWB from all pairs,

$\displaystyle\hat{A}_{\textrm{gw}}^{2}$
$\displaystyle=\frac{\sum_{a}\sum_{b>a}\rho_{ab}\Gamma_{ab}\sigma_{ab}^{-2}}{\sum_{a}\sum_{b>a}\Gamma_{ab}^{2}\sigma_{ab}^{-2}}\,,$
(65) $\displaystyle\sigma_{\textrm{gw}}^{2}$
$\displaystyle=\left(\sum_{a}\sum_{b>a}\sigma_{ab}^{-2}\Gamma_{ab}^{2}\right)^{-1}.$
(66)

This statistic is often used for null hypothesis testing for GWB detection.
Under the (null) assumption that $A^{2}_{\textrm{gw}}=0$, the variance of this
estimator can be used to construct a signal-to-noise ratio (S/N),

S/N
$\displaystyle=\frac{\sum_{a}\sum_{b>a}\rho_{ab}\Gamma_{ab}\sigma_{ab}^{-2}}{\left(\sum_{a}\sum_{b>a}\Gamma_{ab}^{2}\sigma_{ab}^{-2}\right)^{1/2}}\,,$
(67)

which we calculate on the data and then compare to its expected distribution
under the null hypothesis. The distribution for this statistic is a
generalized $\chi^{2}$ distribution [83], but it is often estimated
empirically using methods that destroy correlations, but preserve potential
mismodelling. Two such methods are sky scrambles [84] and phase shifts [85].

Construction of the matrices in Equation 63 and Equation 64 requires a choice
of hyperparameters $\bm{\eta}$. Specifically, the red noise parameters for
each pulsar are used to construct $\bm{D}_{a}$, and the CURN amplitude and
spectral index are used to construct $\bm{\phi}_{ab}$. The natural choice for
these parameters are those taken from the Bayesian analysis. However, choosing
the maximum likelihood parameters for $\bm{\eta}$ from a fully Bayesian run,
or jointly maximizing the individual 1-D posteriors for each parameter, leads
to a bias in the recovered value of $A_{\textrm{gw}}^{2}$ [34]. Part of this
bias is due to making a single choice of noise parameters. By averaging the
statistic calculated over draws of $\bm{\eta}$ from a posterior chain,
resulting in what is referred to as the “noise marginalized optimal statistic”
(NMOS) this bias can be partially alleviated. Additionally, a single choice of
hyperparameters could result in a larger value of the S/N than is
representative of the dataset. In general, therefore, the S/N is averaged over
many draws from $p(\bm{\eta}|\bm{\delta t})$, and this average is used as a
detection statistic. Other approaches have also been proposed, e.g., averaging
the $p$-value associated with the S/N for each draw, instead of averaging the
S/N [35].

The OS defined in Equation 65 and Equation 66. It is implemented in the
compute_os method of the OptimalStatistic class, which is found in the
frequentist.optimal_statistic module of enterprise_extensions. The NMOS
(described in the previous paragraph) is implemented as the
compute_noise_marginalized_os of the same class, and takes a MCMC chain and a
list of parameter names as input.

###  IV.2 Optimal statistic with a non-negligible GWB

In the case where $A_{\textrm{gw}}^{2}$ is comparable to the red noise level
in some pulsars, the assumption that
$\langle\rho_{ab}\rho_{cd}\rangle\propto\delta_{ac}\delta_{bd}$ breaks down.
We must account for the covariance between correlations when constructing both
$\hat{A}^{2}_{\textrm{gw}}$ and especially its variance, which will be
dominated by the background itself. When accounting for the covariance between
correlations due to the GWB we find

$\displaystyle\bm{\Sigma_{ab,cd}}=$
$\displaystyle\langle\rho_{ab}\rho_{cd}\rangle-\langle\rho_{ab}\rangle\langle\rho_{cd}\rangle$
(68) $\displaystyle=$
$\displaystyle\operatorname{tr}\left(\bm{Q}_{ba}\mathbf{P}_{ac}\bm{Q}_{cd}\mathbf{P}_{db}\right)$
$\displaystyle+\operatorname{tr}\left(\bm{Q}_{ba}\mathbf{P}_{ad}\bm{Q}_{dc}\mathbf{P}_{cb}\right)\,,$
(69)

where

$\displaystyle\mathbf{P}_{ab}$ $\displaystyle=\left<\bm{\delta
t}_{a}\bm{\delta
t}_{b}^{T}\right>=\delta_{ab}\bm{D}_{a}+\bm{F}_{a}\bm{\phi}_{ab}\bm{F}_{b}^{T}\,,$
(70)

We can then construct a least-squares estimator for the background using this
covariance matrix,

$\displaystyle\hat{A}_{\textrm{gw}}^{2}$
$\displaystyle=\frac{\bm{\Gamma}^{T}\bm{\Sigma}^{-1}\bm{\rho}}{\bm{\Gamma}^{T}\bm{\Sigma}^{-1}\bm{\Gamma}}\,,$
(71) $\displaystyle\sigma_{A_{\textrm{gw}}^{2}}^{2}$
$\displaystyle=\left(\bm{\Gamma}^{T}\bm{\Sigma}^{-1}\bm{\Gamma}\right)^{-1}\,,$
(72)

where $\bm{\Sigma}$ is given by Equation 68, $\bm{\rho}$ is a vector of paired
correlations and $\bm{\Gamma}$ is a vector of the HD correlation coefficients
corresponding to each pair. It is important to note here that construction of
$\mathbf{P}_{ab}$ and therefore $\bm{\Sigma}$ requires a choice of
$\bm{\eta}$. That is, some choice of $A^{2}_{\textrm{gw}}$ is needed to
actually construct our estimator. One can take an iterative approach, where we
begin with a choice of $\bm{\eta}$ [e.g., the maximum a posteriori draw from
$p(\bm{\eta}|\bm{\delta t})$], evaluate Equation 71, and then use the
resulting $\hat{A}^{2}_{\textrm{gw}}$ to construct $\bm{\Sigma}$, iterating
until convergence. In practice, we have found this converges rapidly, and is
consistent with results that use a single iteration with an initial choice of
$A^{2}_{\textrm{gw}}$ estimated from the posterior $p(\bm{\eta}|\bm{\delta
t})$.

It is also common to estimate the strength of the observed background using
subsets of paired correlations that have similar separations on the sky. The
individual bins should then trace the HD curve. We collect pairs that fall
into a single bin, initially choosing the number of bins we would like, and
then assigning pulsar pairs to bins such that there are roughly the same
number of pairs in each bin. We label the average angular separation between
pulsars in the $i^{\textrm{th}}$ bin as $\xi_{i}$, and construct an estimator
for the correlated power in each bin, $\rho_{\textrm{opt,}i}$, whose
expectation is given by
$\langle\rho_{\textrm{opt,}i}\rangle=\Gamma_{\xi_{i}}A_{\textrm{gw}}^{2}$. In
this case $\Gamma_{\xi_{i}}$ is the HD correlation coefficient evaluated at
the average angular separation for pulsar pairs in the $i^{\textrm{th}}$ bin.
Other choices could also be made, and would slightly change the results [33].
Motivated by our choice of binning, we can imagine $\bm{\Sigma}$ taking a
block form where the $i,i$ block corresponds to correlations between pairs in
bin $\xi_{i}$, the $i,j$ block corresponds to correlations between pairs in
$\xi_{i}$ and pairs in $\xi_{j}$. The resulting estimator is given by [33]

$\displaystyle\rho_{\textrm{opt,}i}=\Gamma_{\xi_{i}}\frac{\bm{\Gamma}_{i}^{T}\bm{\Sigma}_{ii}^{-1}\bm{\rho}_{i}}{\bm{\Gamma}_{i}^{T}\bm{\Sigma}_{ii}^{-1}\bm{\Gamma}_{i}},$
(73)

where $\bm{\Gamma}_{i}$ is a vector of the overlap reduction function for all
pairs in bin $i$, $\bm{\rho}_{i}$ is the vector of paired correlations in bin
$i$. The non-zero GWB induces correlations between these bins as well, and
this covariance matrix is given by [33]

$\displaystyle B_{ij}$
$\displaystyle=\langle\rho_{\textrm{opt,}i}\rho_{\textrm{opt,}j}\rangle-\langle\rho_{\textrm{opt,}i}\rangle\langle\rho_{\textrm{opt,}j}\rangle,$
(74)
$\displaystyle=\Gamma_{\xi_{i}}\Gamma_{\xi_{j}}\frac{\bm{\Gamma}_{i}^{T}\bm{\Sigma}_{ii}^{-1}\bm{\Sigma}_{ij}\bm{\Sigma}_{jj}^{-1}\bm{\Gamma}_{j}}{\left(\bm{\Gamma}_{i}^{T}\bm{\Sigma}_{ii}^{-1}\bm{\Gamma}_{i}\right)\left(\bm{\Gamma}_{j}^{T}\bm{\Sigma}_{jj}^{-1}\bm{\Gamma}_{j}\right)}\,.$
(75)

###  IV.3 Optimal statistic for multiple correlation patterns

Reference [38] arrived at the multiple component optimal statistic (MCOS)
through a $\chi^{2}$ approach, with

$\displaystyle\chi^{2}=\sum_{ab}\left(\frac{\rho_{ab}-\sum_{\alpha}A^{2}_{\alpha}\Gamma_{ab}^{\alpha}}{\sigma_{ab}}\right)^{2},$
(76)

$\rho_{ab}$ given by Equation 61, and $\sigma_{ab}$ its associated
uncertainty. The label $\alpha$ now indexes different spatial correlation
patterns. For example, one can jointly minimize $\chi^{2}$ with respect to
$A^{2}$ for both a HD-correlated process and a monopole process, and calculate
the associated covariance matrix between those estimators.

By generalizing the optimal statistic to include more than one ORF
simultaneously, one arrives at the multiple-correlation optimal statistic
(MCOS),

$\hat{A}_{\alpha}^{2}=\sum_{\beta}B_{\alpha\beta}C^{\beta}\,,$ (77)

where $\alpha$ and $\beta$ are individual ORFs, and

$\displaystyle B^{\alpha\beta}$
$\displaystyle\equiv\sum_{a}\sum_{b>a}\frac{\Gamma_{ab}^{\alpha}\Gamma_{ab}^{\beta}}{\sigma_{ab}^{2}}\,,$
(78) $\displaystyle C^{\beta}$
$\displaystyle\equiv\sum_{a}\sum_{b>a}\frac{\rho_{ab}\Gamma_{ab}^{\beta}}{\sigma_{ab}^{2}}\,.$
(79)

When denoting matrices in this notation, upper and lower indices indicate
matrix inverses with respect to one another. The variance on the individual
estimators for each spatially-correlated process is given by

$\sigma_{\hat{A}_{\alpha}^{2}}^{2}=B_{\alpha\alpha}\,.$ (80)

The noise-marginalized version of the MCOS follows a similar structure as the
noise-marginalized version of the original optimal statistic. By drawing
parameters from the MCMC chains, we average over the noise and return a
distribution of values.

The MCOS is implemented in the compute_multiple_corr_os method of the
OptimalStatistic class. The noise marginalized version of the MCOS is
implemented as the compute_noise_marginalized_multiple_corr_os method.

##  V Tests of Bayesian Methods: PTMCMCSampler and Enterprise

Here, we perform tests to assess the code performance for robust Bayesian
inference. We begin these tests by checking prior recovery. This tests that
our proposal distributions satisfy detailed balance, a condition required for
the chain samples to reflect the posterior. Next, we create simulations based
on the 15-year NANOGrav data and check for unbiased posteriors. Simulations
are produced directly from Enterprise models using the TOAs from the 15-year
data. Therefore, good posterior recovery implies that the models we use are
self-consistent and that the recovered posteriors are in the correct place
with the right width. Finally, we use a reduced version of the data to check
that Bayes factors agree among different methods of computation. If the
different methods agree, we conclude that our calculations are working
properly.

###  V.1 Prior Recovery Tests

To test the proposal distributions incorporated in PTMCMCSampler and
enterprise_extensions, we sample a posterior that is equal to the priors by
setting the likelihood equal to the prior and the prior to a constant. PT only
tempers the likelihood in PTMCMCSampler, so it is necessary to sample the
prior as the likelihood to also test this part of the sampler. If the proposal
distributions satisfy detailed balance, then the recovered posterior equals
the input priors within sampling uncertainties. PTMCMCSampler contains three
default proposal distributions known inside the code as “jump proposals.”
Additional jump proposals come from enterprise_extensions, and which proposals
get used depends on the type of search being performed. The isotropic GWB
analyses includes AM, SCAM, and DE proposals. Additionally, this search
includes prior draws and two-dimensional empirical distributions on the
intrinsic red noise amplitudes and spectral indices for each pulsar.

To assess prior recovery, we use a quantile-quantile (Q-Q) plot. Q-Q plots
compare the quantiles of the distribution of our recovered prior with samples
drawn from a simulated distribution of the input. We subtract the mean of the
simulated distribution from every point in our plot so that the mean falls
along zero on the vertical axis. The expected result is a set of lines, one
for each parameter, falling within the given uncertainties with few venturing
outside the $3\sigma$ bounds. Bias could appear as a line remaining
significantly above or below the mean for the entire interval indicating that
more samples exist in one quantile of the distribution.

Figure 3:  Quantile-quantile plot showing the recovery of the prior. To
produce this plot, we sample the prior in place of the likelihood and set the
prior to a constant. This is not the same as setting the likelihood to a
constant, because we use a parallel tempering method that only tempers the
likelihood. Priors on the $\gamma$ parameters are Uniform$[0,7]$, and the
priors on the $\log_{10}A$ parameters are Uniform$[-18,-11]$. The chains
produced should be equal to the input prior distribution within the sampling
uncertainties after thinning. On the horizontal axis, we plot the quantiles of
the simulated uniform distribution. On the vertical axis, we plot the
quantiles of the distribution output from the sampling process minus the mean
of simulated prior quantiles so that the mean lies along the zero of the
vertical axis. The curved, solid lines show $1\sigma$, $2\sigma$, and
$3\sigma$ uncertainties. The uncertainty lines are created by taking the
average and standard deviation over 10,000 realizations of a uniform
distribution with the same number of samples as the observed distribution.
This plot shows that the priors were recovered correctly when using parallel
tempering, AM, SCAM, DE, empirical, and prior proposal distributions.

We use the same prior for each spectral index parameter $\gamma$ and for each
log amplitude parameter $\log_{10}A$,

$\displaystyle\gamma$ $\displaystyle\sim\text{Uniform}[0,7]\,,$ (81)
$\displaystyle\log_{10}A$ $\displaystyle\sim\text{Uniform}[-18,-11]\,.$ (82)

These priors are chosen specifically for this test and may be slightly
different in “production grade” analyses between the intrinsic red noise
amplitudes and the GWB amplitude, because we can rule out very large GWB
amplitudes since they have not been observed in previous data sets. Spectral
index priors are typically the same for all parameters. The NANOGrav 15-year
GWB analysis uses 67 pulsars. Assuming a power-law spectrum, each pulsar adds
a spectral index and an intrinsic red noise amplitude parameter. Along with
the common process spectral index and amplitude parameters, each of the plots
in Figure 3 contains 68 lines. The majority of the lines remain inside three
sigma uncertainty bounds and only occasionally do lines venture outside and
then back toward zero.

###  V.2 Injection Recovery Tests

Next, we create simulations to check that PTMCMCSampler returns unbiased
posteriors. To quantify whether our posteriors are unbiased, we use a P-P
plot444We could use Q-Q plots, which are quite general, for this too. However,
P-P plots are somewhat standard among the literature for testing injection
recovery.. We make these plots by creating 100 realizations of data, sampling
the posteriors, and finding the $p$-value at which the injected value falls.
These $p$-values should follow a uniform distribution. By taking the CDF of
these values and plotting them against the $p$-values, we expect each
parameter (shown as a line on the plot) to follow a diagonal line within some
confidence interval. Since we model the red noise as a power law and have 67
pulsars, there are 136 parameters searched over in each simulation: an
amplitude and spectral index for each pulsar’s intrinsic red noise and an
amplitude and spectral index for the common red process which is shared among
all pulsars. For this P-P plot, we use the CURN model, because the
computational expense is small compared to the HD correlated model, and the
posteriors between the two models are similar. Furthermore, Hourihane  _et
al._ [76] found that reweighting the CURN plots using the HD model’s
likelihood also resulted in diagonal P-P plots. Given the similarity between
the two models’ posteriors, PTMCMCSampler should have no extra problems with
exploring the HD parameter space given its effective exploration of the CURN
parameter space.

Figure 4: A P-P plot shows the recovery of injected parameters for 100
realizations simulated using the 15-year NANOGrav data with injections pulled
from the priors in Table 1. The model being simulated includes power law
intrinsic red noise parameters and a power law CURN. Each of the 136
parameters are presented as an individual line on the P-P plot. A diagonal
black line indicates perfect recovery, and the 68.27%, 95.45%, 99.73%
confidence intervals, found using an inverse CDF of a binomial distribution
[86], appear as curved black lines with 68.27% being the closest to the
diagonal and 99.73% being farthest away. The plot indicates no significant
bias in recovery of the posterior by PTMCMCSampler when the simulation and
recovery is performed via Enterprise.

We make simulations differently here than in other papers which have used
either libstempo [87] and TEMPO2 or PINT. The question of whether our models
match the simulations of realistic PTA data is outside the scope of this
study. Instead, we use the models in Enterprise to simulate pulsar residuals
with timing model uncertainties, specified white noise, intrinsic red noise,
and either a CURN or a HD correlated GWB. The priors set on these values are
as shown in Table 1. Crucially, the distributions which we draw values from
and the priors we search over must be the same. If they are not, the P-P test
will fail. We have reduced the prior space from the full production analyses
to limit ourselves to a detectable part of the parameter space. While this is
not required to make good P-P plots, our purpose is to test whether we get
unbiased results in the event of a set of parameters which have strong signals
visible.

Parameter | Prior | Values
---|---|---
IRN Amplitude | Uniform | $[-15,-12]$
IRN Spectral Index | Uniform | $[2,6]$
CURN Amplitude | Uniform | $[-16,-14]$
CURN Spectral Index | Uniform | $[2,6]$
Table 1: Priors sampled for the injected values of simulation realizations.
IRN indicates an intrinsic red noise parameter and CURN indicates a common
uncorrelated parameter which is common to all pulsars. These priors are also
used in the search of the simulations to recover posteriors for use in P-P
tests.

The P-P plot produced in Figure 4 indicates little bias in the recovery of the
injected values. Lines that remain above and below the diagonal could indicate
that the recovery is biased so that we find the injected values either too low
or too high consistently. A signature “S” shape in which the line goes above
(below) and then below (above) the diagonal indicates that the width of the
posterior distribution is over or under estimated. After checking each of
these 136 parameters individually, we find no evidence of either of these
issues.

###  V.3 Bayes Factor Recovery Tests

In the final segment of the Bayesian tests, we check how well our Bayes
factors are recovered. Here, we use a few different techniques of computing
the Bayes factors and compare between them. Due to computational limitations,
we reduce our data set to a set of 14 pulsars that have been timed for greater
than 15 years. This allows us to use nested sampling which does not converge
quickly for high dimensional spaces such as with the full 67 pulsar parameter
space. Simulations were made with varying noise realizations and GWB
injections across a broad range of log Bayes factors,

$\log_{10}\text{BF}=\log_{10}Z_{\text{HD}}-\log_{10}Z_{\text{CURN}}\,,$ (83)

from $-2.5$ to $17$, where $Z$ is the evidence. Bayes factors are computed for
the HD correlated model against the CURN model. We check the hypermodel using
PTMCMCSampler against nested sampling with UltraNest, and the hypermodel
results against Bayes factors returned with reweighting [76]. By running these
methods on the same 100 realizations with each method, we show that they
return consistent answers, although at different levels of uncertainty.

In the comparison between reweighting and the hypermodel, we take the
resulting chain from the hypermodel and reweight any of the uncorrelated model
samples to the HD correlated model. This gives us a set of weights which can
be used to compute a Bayes factor and uncertainties. We find that the two
methods give results on the log ratio that are consistent with zero in every
realization within $3\sigma$. In every case, reweighting gives a larger
uncertainty than the hypermodel and dominates in the subtraction of Bayes
factors in which errors are propagated in quadrature.

Nested sampling requires a stopping condition in terms of the uncertainty on
the log evidence. Unfortunately, we find that even with the reduced parameter
space, setting the stopping condition to $d\log Z<0.5$ led to week-long run
times. Once again, the uncertainties on the hypermodel are overwhelmed by the
uncertainties on nested sampling without requiring more computational
resources for nested sampling, and all realizations are consistent with zero
within $3\sigma$.

Figure 5: Logarithmic ratio of the Bayes factor computed for each of 100
simulations using the HyperModel framework and reweighting in the top panel
and the HyperModel framework and nested sampling in the bottom panel. Each
point indicates a mean value of the ratio and the $1\sigma$ uncertainties are
given as a vertical bar on each point. The red line indicates zero on the
vertical axis, where we expect these values to fall if the Bayes factors
returned from each method are consistent. Uncertainties are dominated by the
nested sampling and reweighting methods. The values are consistent with zero
within $3\sigma$ across all simulations.

As a final test of the Bayes factors, we run the hypermodel on a CURN model
against itself on the same 100 simulations that were used above. In this case,
we know that the Bayes factor must equal one, because a model should not be
preferred over itself. On top of checking whether we get an answer consistent
with the known value, this test shows whether the uncertainties are being
estimated properly. As shown in Figure 6, the Bayes factor of one is recovered
in every realization within $3\sigma$. This method represents an easy check
that can be performed for any situation to make sure that the hypermodel
calculation is working. However, this test is not sufficient to claim that the
method will work for all scenarios. The case of a model against itself does
not take into account the situation where the posteriors between two models
are very different. Therefore, the test of consistency against other samplers
remains necessary.

Figure 6: Bayes factors computed for a CURN model against itself for the same
100 simulations as used in the other tests of Bayes factors. The red line in
this figure indicates the true Bayes factor between the two models. All
realizations recover the true value within $3\sigma$. This represents a quick
check to see if the Bayes factor calculation returns correct answers in the
ideal scenario of posteriors that are exactly the same between the two models
being sampled.

##  VI Tests of frequentist methods: enterprise_extensions and the optimal
statistic

In this section we present a series of tests on the optimal statistic
presented in Section IV. We begin by using simulated data sets to compare the
“traditional” optimal statistic and the one that accounts for covariance
between pulsar pair correlations. We use those same data sets to evaluate how
the revised binned estimator performs as well. This is followed by a
discussion of the distinction between using these statistics as estimators for
the amplitude of the GWB vs. using them as detection statistics in a classical
null hypothesis testing scenario. We finish by summarizing recent work in the
literature on the multiple component optimal statistic and constructing
empirical and analytic distributions for the optimal statistic.

###  VI.1 The optimal statistic as an estimator

In the case where $A_{\textrm{gw}}^{2}$ is small compared to intrinsic red
noise in all pulsars, the distribution on $\hat{A}^{2}_{\textrm{gw}}$ in
Equation 65 is approximately Gaussian with a variance given by Equation 66
except in the tails [83]. In present analyses, $A^{2}_{\textrm{gw}}$ is not
smaller than the corresponding intrinsic red noise amplitude for at least a
few pulsars [19, 20]. Therefore, to test how well the small signal
approximation works, we perform 200 simulations with the length and cadence of
the data set given in [50] and with red noise drawn from
$p(\bm{\eta}|\bm{\delta t})$ from [8]. We draw $A_{\textrm{gw}}$ from a
uniform distribution, $\log_{10}A_{\textrm{gw, inj}}\in[-17,-13]$. On each
simulated dataset we calculate $\hat{A}^{2}_{\textrm{gw}}$ and its variance
using Equation 65 and Equation 66, and we use the method described in Appendix
B to construct PP-like plots, replacing the posterior samples with a Gaussian,
$\mathcal{N}(\hat{A}_{\textrm{gw}}^{2},\sigma^{2}_{\textrm{gw}})$. These are
not traditional PP-plots, as there is no well-defined prior from which we draw
our simulations and sample our posterior. However, performing the same
processing as in Appendix B does test how frequently the Gaussian distribution
centered on $\hat{A}^{2}_{\textrm{gw}}$ with variance
$\sigma_{\textrm{gw}}^{2}$ includes the injected value of
$A^{2}_{\textrm{gw}}$. We do the same thing for the estimator that includes
covariances between pulsar pairs, defined in Equation 71 and Equation 72.

The results of this test are shown in Figure 7, with the “traditional” optimal
statistic results shown in the blue, solid curve and the corrected results
shown in the orange, solid curve. The blue curve is consistent with the
traditional optimal statistic underestimating the error on the estimator by
not accounting for the GWB, and therefore not capturing the injected value in
its credible intervals as frequently as it should. The orange curve corrects
this, as it follows the expected line more closely. Therefore, if one is to
use the optimal statistic as an estimator for the GWB, the corrected statistic
performs significantly better.

Figure 7: P-P-like plot that characterizes how well the optimal statistic
functions as an estimator for the amplitude of the GWB. The blue, solid curve,
which uses the traditional OS does not follow the expected diagonal line,
indicating that it underestimates the variance on the estimator for the
background. The orange curve, meanwhile, follows the diagonal line with its
expected error bars (shaded region) because it properly estimates the variance
by including the contribution from the GWB.

We can perform the same test for the binned optimal statistic in Equation 73
and Equation 74 as well. For each simulation, we evaluate the cumulative
distribution function at the injected value, similar to the test in Appendix
B, but on each individual binned estimator in Equation 73 and Equation 74. We
plot the results of the P-P-like test in Figure 8, where we see similar
results to Figure 7. In this case, we show deviations from the predicted line,
and so expectation is zero, as opposed to $y=x$. The blue curves show excess
cases where the CDF is zero and close to one, indicating misestimation. The
orange curves, which includes the GWB in its variance, perform better,
especially near zero and one. The orange curves do reach the 3-sigma level
more than one might normally expect, and this is likely because we have
assumed that the distribution on the binned estimators is a Gaussian. In
practice, the distribution on the estimators is a generalized chi-squared
distribution, which can be well approximated by a Gaussian under certain
circumstances [83].

In these simulations we have used the injected amplitude of the background
when calculating the covariance between correlations. We do this for practical
purposes–the goal of these P-P-like tests is to show that our estimator is
unbiased and its error bars are correct when we take into account the
amplitude of the GWB. In practice, we do not have access to the GWB amplitude
a priori–we either use values drawn from a MCMC chain (e.g. the noise-
marginalized optimal statistic), or we could employ an iterative approach,
where we calculate the GWB estimator using the optimal statistic, and then use
the estimator in the covariance matrix to properly estimate the covariance
matrix from the correlations, and then repeat until convergence.

Figure 8: P-P-like plots (with diagonal subtracted off) that characterize how
well the angular-binned optimal statistic functions as an estimator for the
amplitude of the GWB, modulated by the Hellings-Downs correlations. We show
results estimating the amplitude in each of 11 individual bins. Top: the blue,
solid curve, uses the traditional binned OS does not follow the expected
horizontal line, indicating that it underestimates the variance on the
estimator for the background. This is especially obvious looking near zero and
near one, where we see the curves diverge from the expected range (shaded
regions) Bottom: the orange curve follows the horizontal line better than the
blue curves in the top, because it properly estimates the variance by
including the contribution from the GWB.

###  VI.2 The optimal statistic for detection

We have shown that the estimators which include the strength of the GWB in
their construction are better than the ones that do not, especially in the
case where the amplitude of the GWB is of a similar size to that of the
intrinsic red noise. However, the S/N calculated in Equation 67 is calculated
under the null hypothesis. Therefore, when making a detection, we construct a
distribution for this statistic under the null hypothesis (i.e., no correlated
power), and then we compare the same statistic calculated on the original data
to that distribution. Therefore, Equation 67 is the correct expression for a
detection statistic. It is common to use a noise-marginalized version of this
statistic, as discussed previously. That is, we take the average of this S/N
over many draws from $p(\bm{\eta}|\bm{\delta t})$, and compare this to a null
distribution.

Construction of a distribution for the null statistic takes a few forms. The
analytic distribution of this statistic was calculated in [83], and can be
used. However, pathologies in the data that are not correctly modeled would
not be accounted for in an analytic calculation. To preserve potential
mismodeling, but still approximate the null distribution of the detection
statistic, it is common to use sky scrambling [84] and phase shifting [85].
Sky scrambling involves assigning random sky positions to each pulsar, and
calculating the optimal statistic, while phase shifting applies random phases
to each frequency and each pulsar in the $\bm{F}$-matrix used to construct
Equation 64. Comparisons between the analytic distribution and distributions
constructed by sky scrambling and phase shifting are still in progress.
Additionally, some concerns have been raised about whether performing sky
scrambling and phase shifting without producing independent555A statistic one
could use to assess independence is presented in [88] scrambles or shifts
results in oversampling parts of the null distribution, and undersampling
other parts [88]. Some tests of this have been done to explore this in [19],
and found that more stringent “orthogonality” conditions do not produce
meaningfully different distributions than those calculated with less stringent
conditions. However, more tests are in progress.

###  VI.3 MCOS Tests

We do not present novel tests of the MCOS here, as tests have been presented
in separate work that more completely presents the statistic itself and
characterizes its behavior [38]. We summarize those results here, for
completeness.

In [38] and [89] it is shown that, due to the non-isotropic distribution of
pulsars on the sky, Hellings-Downs correlations are not orthogonal to
monopolar or dipolar correlations, as one might expect. This is what motivates
the construction of a statistic that simultaneously fits for multiple spatial
correlation patterns simultaneously. In [38], they perform injections of non-
Hellings-Downs correlation patterns, and show that it is possible using the
statistic in Equation 67 to find a spurious detection of Hellings-Downs
correlations. Using Equation 77 and Equation 80, remedies this situation, as
the new statistic correctly finds that the non-Hellings-Downs correlations are
preferred. Additionally, the authors consider the MCOS as an estimator for the
strength of the background associated with each spatial correlation pattern.
They find that, when multiple correlation patterns are injected into the data,
the MCOS estimator for the strength of each pattern performs better than
individually estimating the strength of each process separately. However,
there are issues with the estimation of signals and their uncertainties owing
to the fact that the MCOS is based on the traditional optimal statistic, which
does not take into account covariance between pulsar pairs (see subsection
VI.1).

##  VII Conclusions

Here we have checked both Bayesian and frequentist methods used in the 15-year
GWB analysis. These methods were outlined as they are used in the 15-year
analysis. We subjected the software implementations of each part of the
analysis to the most stringent tests of correctness performed so far, and we
find that each analysis returns unbiased, self-consistent results.

In the Bayesian tests, we find that the priors are recovered correctly while
including parallel tempering, AM, SCAM, DE, empirical distributions, and prior
proposal distributions that are used in the full analysis. This indicates that
all proposal distributions work properly and are not biasing our results. We
also find that the simulations made with injections pulled from a prior
distribution return diagonal P-P plots within acceptable uncertainties. By
creating realizations of the 15-year data with various amplitudes and spectral
indices of a GWB including HD correlations in them, we perform tests using
reweighting, nested sampling, and the hypermodel to return Bayes factors
between an HD correlated and CURN model. Each comparison between these methods
return log Bayes factor ratios consistent with zero. Finally, we perform model
comparison between a model and itself on these same simulations and find the
expected result of Gaussian distributed Bayes factors centered around unity.

In the frequentist tests, we show through simulating data sets that properly
including the GWB when constructing the optimal statistic is necessary in the
situation where the GWB is not small compared to the intrinsic red noise. The
estimator that accounts for the size of the GWB, however, is not consistent
with the null hypothesis of no correlated power in the timing residuals, and
therefore the “traditional” optimal statistic should be used as a detection
statistic. Such a detection statistic should be calculated on the observed
timing residuals and compared to a null distribution, which can be calculated
either using the analytic distribution [83] or a method that preserves
potential mismodeling but suppresses correlations [84, 85]. Finally, we
summarize the recently-proposed MCOS, which simultaneously fits for multiple
spatial correlation patterns, and helps prevent, e.g., monopolar correlations
from producing a spurious detection of Hellings-Downs correlations and vice-
versa.

Currently, the NANOGrav software collection allows one to perform inference
without a steep learning curve. As our data volume and parameter space
increases with each data set, we look for additional methods to improve
efficiency of memory management and the likelihood evaluation speed in
Enterprise. This includes possible data compression techniques to reduce the
number of points of data required to fit to perform inference (e.g, [90]).
Additionally, reducing the auto-correlation of our chains produced with
PTMCMCSampler could reduce computation time significantly. To this end, some
recent work has used JAX [91], a Python package that includes just-in-time
compilation (JIT) to speed up evaluation of loops and auto-differentiation to
take derivatives quickly and accurately, among other convenient features. With
this, Freedman  _et al._ [92] were able to use a No-U-Turn sampler (NUTs), a
Hamiltonian Monte Carlo (HMC) sampling technique, which reduced auto-
correlation significantly compared to PTMCMCSampler. In Bécsy  _et al._ [93],
JIT is used through the numba Python package alongside techniques to sample
some parameters more often than others to vastly increase the speed of
evaluation of continuous GW searches. While these papers have shown the
incredible speed gains that one may achieve via these new technologies, they
cannot be easily implemented in the current paradigm used by Enterprise. This
is primarily due to sub-classing of numpy arrays used in Enterprise, which is
not supported by these other programs.

Both the Bayesian and frequentist methods and their software representations
pass all the tests they were subjected to. As PTA data sets increase in
sensitivity, testing that our methods are reliable proves paramount. Here, we
validate the results of current and future analyses that use the methods
examined. Through rigorous testing of our software and methods, we form a
foundation on which astrophysical results may stand.

###### Acknowledgements.

The NANOGrav project receives support from National Science Foundation (NSF)
Physics Frontiers Center award number 1430284. A.D.J., K.C., and M.V.
acknowledge support from the Caltech and Jet Propulsion Laboratory President’s
and Director’s Research and Development Fund. A.D.J. and K.C. acknowledge
support from the Sloan Foundation. S.H. is supported by the National Science
Foundation Graduate Research Fellowship under Grant No. DGE-1745301. L.B.
acknowledges support from the National Science Foundation under award
AST-1909933 and from the Research Corporation for Science Advancement under
Cottrell Scholar Award No. 27553. P.R.B. is supported by the Science and
Technology Facilities Council, grant number ST/W000946/1. S.B. gratefully
acknowledges the support of a Sloan Fellowship, and the support of NSF under
award #1815664. M.C. and S.R.T. acknowledge support from NSF AST-2007993. M.C.
and N.S.P. were supported by the Vanderbilt Initiative in Data Intensive
Astrophysics (VIDA) Fellowship. Support for this work was provided by the NSF
through the Grote Reber Fellowship Program administered by Associated
Universities, Inc./National Radio Astronomy Observatory. Support for H.T.C. is
provided by NASA through the NASA Hubble Fellowship Program grant #HST-
HF2-51453.001 awarded by the Space Telescope Science Institute, which is
operated by the Association of Universities for Research in Astronomy, Inc.,
for NASA, under contract NAS5-26555. M.E.D. acknowledges support from the
Naval Research Laboratory by NASA under contract S-15633Y. T.D. and M.T.L. are
supported by an NSF Astronomy and Astrophysics Grant (AAG) award number
2009468. E.C.F. is supported by NASA under award number 80GSFC21M0002. G.E.F.,
S.C.S., and S.J.V. are supported by NSF award PHY-2011772. The Flatiron
Institute is supported by the Simons Foundation. The work of N.La. and X.S. is
partly supported by the George and Hannah Bolinger Memorial Fund in the
College of Science at Oregon State University. N.La. acknowledges the support
from Larry W. Martin and Joyce B. O’Neill Endowed Fellowship in the College of
Science at Oregon State University. Part of this research was carried out at
the Jet Propulsion Laboratory, California Institute of Technology, under a
contract with the National Aeronautics and Space Administration
(80NM0018D0004). M.A.M. is supported by NSF #1458952 and #2009425. C.M.F.M.
was supported in part by the National Science Foundation under Grants No. NSF
PHY-1748958 and AST-2106552. A.Mi. is supported by the Deutsche
Forschungsgemeinschaft under Germany’s Excellence Strategy - EXC 2121 Quantum
Universe - 390833306. K.D.O. was supported in part by NSF Grant No. 2207267.
T.T.P. acknowledges support from the Extragalactic Astrophysics Research Group
at Eötvös Loránd University, funded by the Eötvös Loránd Research Network
(ELKH), which was used during the development of this research. S.M.R. and
I.H.S. are CIFAR Fellows. Portions of this work performed at NRL were
supported by ONR 6.1 basic research funding. J.D.R. also acknowledges support
from start-up funds from Texas Tech University. J.S. is supported by an NSF
Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2202388,
and acknowledges previous support by the NSF under award 1847938. Pulsar
research at UBC is supported by an NSERC Discovery Grant and by CIFAR. S.R.T.
acknowledges support from an NSF CAREER award #2146016. C.U. acknowledges
support from BGU (Kreitman fellowship), and the Council for Higher Education
and Israel Academy of Sciences and Humanities (Excellence fellowship). C.A.W.
acknowledges support from CIERA, the Adler Planetarium, and the Brinson
Foundation through a CIERA-Adler postdoctoral fellowship. O.Y. is supported by
the National Science Foundation Graduate Research Fellowship under Grant No.
DGE-2139292. J.A.C.C. was supported in part by NASA CT Space Grant PTE Federal
Award Number 80NSSC20M0129, and also supportedin part by the National Science
Foundation’s NANOGrav Physics Frontier Center, Award Number 2020265.

##  Appendix A Asynchronous Parallel Tempering and the HyperModel Framework

PTMCMCSampler previously used an asynchronous parallelization model for
parallel tempering in which each chain sampled a set number of steps on its
own unless interrupted by another chain asking to propose a swap between
chains. During our tests, we noticed that this asynchronous parallel-tempered
MCMC biased Bayes factors computed with hypermodel in favor of the model that
took the longest evaluation time per iteration. As a simple example of this
issue, we simulate data with a sinusoidal signal $h(t)$ and noise $n(t)$,

$\textbf{d}(t)=\textbf{h}(t)+\textbf{n}(t)\,.$ (84)

The sinusoid signal is described by amplitude, angular frequency, and phase,

$h_{i}(t,A,\omega,\phi)=A\sin(\omega t_{i}+\phi)\,,$ (85)

and the noise $n(t)\sim\mathcal{N}(0,1)$. The log likelihood we use for this
signal model is

$p(d\mid
A,\omega,\phi)\propto-\frac{1}{2}\sum_{i}\left(d_{i}(t)-h_{i}(t,A,\omega,\phi)\right)^{2}\,.$
(86)

Priors are all chosen to be uniform with

$\displaystyle A$ $\displaystyle\sim\text{Uniform}[0,5]\,,$ (87)
$\displaystyle\omega$ $\displaystyle\sim\text{Uniform}[0,3]\,,$ (88)
$\displaystyle\phi$ $\displaystyle\sim\text{Uniform}[0,\pi]\,.$ (89)

Priors remain the same across each model used, and parameter estimation of the
injected values returns consistently with the posteriors.

The HyperModel framework consists of multiple models concatenated with a
continuous “switch” parameter between them. This framework chooses models
based on which of two bins the switch parameter falls into. In this case, the
first model is the sinusoidal model as described above, and the second model
is the same model, but with a time.sleep() call to increase the evaluation
time by a factor of a few.

Upon computing the odds ratios with the asynchronously updated sampler, we
find that the odds ratio is not consistent with the anticipated value of one,
as seen in Figure 9. Parallel tempering swap proposals occur every 100
samples. Thinning by multiples of 10 cause the samples to become increasingly
contaminated with swaps that pull the odds ratio away from one.

Figure 9: The Bayes factor for a model over itself should be one, which is
shown as a dashed, black line on this plot. Instead, we find that the Bayes
factor is inconsistent with a value of one. Parallel tempering swaps are
proposed every 100 samples, labeled “Tskip” in this plot (as it is in
PTMCMCSampler). Thinning increases the contamination from swaps by increasing
the number of samples that come from swaps out of the total number.
Uncertainties on the Bayes factors were computed via bootstrapping.

Syncing the model evaluation times gives odds ratios consistent with one,
contrary to what we see in Figure 9 Therefore, the problem appears to be
caused by the evaluation time of one of the models being much longer than the
other. The exact timeline of when these swaps are proposed is shown in Figure
10 and proves critical to figuring out what went wrong.

Figure 10: Two chains, one at $T_{1}=1$ and the other at $T_{1}=2$, progress
through iterations. Iterations start at the white circles and last until the
next circle. Two models are being considered, one of which takes 10 times
longer to evaluate. After each iteration, the hot chain, shown in red on top,
runs a non-blocking probe to check if anything has been sent from the cold
chain, shown in black on bottom. Once the cold chain reaches 100 samples, a
blocking send pushes data required for a swap to the hot chain. This data gets
accepted by the hot chain after the current iteration finishes and a swap is
proposed. Because of this disparity in evaluation times, we find that the hot
chain is 10 times more likely to be in the model that takes longer to
evaluate. This results in a bias in the model that goes into the swap
proposals. Fundamentally, this means that the swaps proposed are more likely
to happen in certain parts of the HyperModel parameter space, breaking
detailed balance.

In asynchronous parallel tempering, swaps are proposed whenever a chain gets
to a set number of samples. The chain sends data to the next chain up in the
temperature ladder, halting the chain until the other chain is done with its
current sample evaluation. While waiting for this signal, the hotter chain
continues collecting samples, probing the lower chain for data after each
sample. Once the signal has been received, the chains swap with their most
recent sample iteration.

Upon proposing a swap, the cold chain finds the hot chain in the model that
has the longer evaluation time more often than not. This means that the swaps
are proposed with a dependence on where in the switch parameter space we are,
violating detailed balance for parallel tempering swap proposals.

One option to fix this involves synchronizing the swap proposals so that all
chains propose swaps at the same iterations. The other option is to weight the
proposals with weights proportional to the ratio of the evaluation times. In
some situations this option is not possible due to dependence of evaluation
time on where we are in the parameter space. We opted to synchronize the
sampler for simplicity and to keep the sampler as generic as possible.

##  Appendix B Understanding Confidence Intervals of P-P Plots

In Figure 4, the confidence intervals are found as in Ashton and Talbot [86].
Unlike in the Q-Q plots where we have tens of thousands of samples, the
standard bootstrap approximated confidence intervals (which estimates the
confidence intervals as Gaussian) overshoot the exact bounds near p-values of
0 and 1. For a given realization of the data, we have an injection for a
single parameter $\theta_{\text{inj}}$, and we can compute a CDF of this
value,

$F(\theta_{\text{inj}})=\int_{-\infty}^{\theta_{\text{inj}}}p(x\mid\bm{\delta
t})dx.$ (90)

If we take many new realizations, their associated CDFs should be uniformly
distributed between 0 and 1. Picking a particular p-value on the horizontal
axis will split the uniform distribution into two segments. Let us define a
success as $F(\theta_{inj})\leq s$ and a failure as $F(\theta_{inj})>s$ where
$s$ is the chosen horizontal axis value. From $N$ draws, the probability of
$k$ successes can be found with a binomial distribution,

$p(k\text{ successes})=\binom{N}{k}s^{k}(1-s)^{N-k}.$ (91)

We would like to have coverage $(1-\alpha)/2$ on both sides of the mean of the
distribution with

$\alpha=(0.6827,0.9545,0.9973)$ (92)

for the three bounds that we show. Using the quantile function for the
binomial distribution, $F^{-1}(s)$ (binom.ppf in scipy.stats) on each of the
horizontal axis values, the offset from the mean is

$\sigma=F^{-1}(s)/N.$ (93)

##  Appendix C Batch Updates of the Sample Mean and Covariance

Let $\bm{x_{0}}$ be an $N\times n_{\text{param}}$ matrix of samples from the
history of an MCMC chain and let $\bm{x_{1}}$ be the next $M$ samples of the
$n_{\text{param}}$ parameters. Finally, let $\bm{x}$ be the concatenation of
$\bm{x_{0}}$ and $\bm{x_{1}}$ with total length $N+M$. The mean for $\bm{x}$
along the $N+M$ axis can be computed as

$\bar{x}^{j}=\bar{x}_{0}^{j}+\frac{1}{M+N}\sum_{i=1}^{N}\left(x_{1}^{ij}-\bar{x}_{0}^{j}\right),$
(94)

where an over bar denotes an average. Next, we can compute the sample
covariance through a batch update as

$\Sigma=\frac{(N-1)\Sigma_{0}+\Sigma_{1}}{N+M-1},$ (95)

where

$\Sigma_{1}^{ij}=\left(x_{1}^{ji}-\bar{x}^{j}\right)\left(x_{1}^{ij}-\bar{x}_{0}^{j}\right).$
(96)

In PTMCMCSampler, these methods use $M=1$ and iterate over the whole chain.

## References

  * Abbott  _et al._ [2016] R. Abbott _et al._ , Observation of gravitational waves from a binary black hole merger, Physical Review Letters 116, 061102 (2016).