This Jupyter notebook can be downloaded from rednoise-fit-example.ipynb, or viewed as a python script at rednoise-fit-example.py.
Red noise, DM noise, and chromatic noise fitting examples
This notebook provides an example on how to fit for red noise and DM noise using PINT using simulated datasets.
We will use the PLRedNoise and PLDMNoise models to generate noise realizations (these models provide Fourier Gaussian process descriptions of achromatic red noise and DM noise respectively).
We will fit the generated datasets using the WaveX and DMWaveX models, which provide deterministic Fourier representations of achromatic red noise and DM noise respectively.
Finally, we will convert the WaveX/DMWaveX amplitudes into spectral parameters and compare them with the injected values.
[1]:
from pint import DMconst
from pint.models import get_model
from pint.simulation import make_fake_toas_uniform
from pint.logging import setup as setup_log
from pint.fitter import WLSFitter
from pint.utils import (
cmwavex_setup,
dmwavex_setup,
find_optimal_nharms,
plchromnoise_from_cmwavex,
wavex_setup,
plrednoise_from_wavex,
pldmnoise_from_dmwavex,
)
from io import StringIO
import numpy as np
import astropy.units as u
from matplotlib import pyplot as plt
from copy import deepcopy
setup_log(level="WARNING")
[1]:
1
Red noise fitting
Simulation
The first step is to generate a simulated dataset for demonstration. Note that we are adding PHOFF as a free parameter. This is required for the fit to work properly.
[2]:
par_sim = """
PSR SIM3
RAJ 05:00:00 1
DECJ 15:00:00 1
PEPOCH 55000
F0 100 1
F1 -1e-15 1
PHOFF 0 1
DM 15 1
TNREDAMP -13
TNREDGAM 3.5
TNREDC 30
TZRMJD 55000
TZRFRQ 1400
TZRSITE gbt
UNITS TDB
EPHEM DE440
CLOCK TT(BIPM2019)
"""
m = get_model(StringIO(par_sim))
[3]:
# Now generate the simulated TOAs.
ntoas = 2000
toaerrs = np.random.uniform(0.5, 2.0, ntoas) * u.us
freqs = np.linspace(500, 1500, 8) * u.MHz
t = make_fake_toas_uniform(
startMJD=53001,
endMJD=57001,
ntoas=ntoas,
model=m,
freq=freqs,
obs="gbt",
error=toaerrs,
add_noise=True,
add_correlated_noise=True,
name="fake",
include_bipm=True,
multi_freqs_in_epoch=True,
)
Optimal number of harmonics
The optimal number of harmonics can be estimated by minimizing the Akaike Information Criterion (AIC). This is implemented in the pint.utils.find_optimal_nharms function.
[4]:
m1 = deepcopy(m)
m1.remove_component("PLRedNoise")
nharm_opt, d_aics = find_optimal_nharms(m1, t, "WaveX", 30)
print("Optimum no of harmonics = ", nharm_opt)
Optimum no of harmonics = 18
[5]:
print(np.argmin(d_aics))
18
[6]:
# The Y axis is plotted in log scale only for better visibility.
plt.scatter(list(range(len(d_aics))), d_aics + 1)
plt.axvline(nharm_opt, color="red", label="Optimum number of harmonics")
plt.axvline(
int(m.TNREDC.value), color="black", ls="--", label="Injected number of harmonics"
)
plt.xlabel("Number of harmonics")
plt.ylabel("AIC - AIC$_\\min{} + 1$")
plt.legend()
plt.yscale("log")
# plt.savefig("sim3-aic.pdf")
[7]:
# Now create a new model with the optimum number of harmonics
m2 = deepcopy(m1)
Tspan = t.get_mjds().max() - t.get_mjds().min()
wavex_setup(m2, T_span=Tspan, n_freqs=nharm_opt, freeze_params=False)
ftr = WLSFitter(t, m2)
ftr.fit_toas(maxiter=10)
m2 = ftr.model
print(m2)
# Created: 2026-06-16T15:30:48.848415
# PINT_version: 1.1.5+169.ge4c9571
# User: docs
# Host: build-33166788-project-85767-nanograv-pint
# OS: Linux-7.0.0-1004-aws-x86_64-with-glibc2.35
# Python: 3.11.14 (main, Apr 27 2026, 17:28:30) [GCC 11.4.0]
# Format: pint
# read_time: 2026-06-16T15:29:59.873849
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM3
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566498843
FINISH 56985.0000000464006366
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1923.9784834863046
CHI2R 0.9831264606470642
TRES 0.99161851118135482043
RAJ 4:59:59.99998969 1 0.00011766275985031068
DECJ 14:59:59.99093491 1 0.01196695415329324629
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999967901 1 5.71130312176685461e-13
F1 -9.9984421525715518555e-16 1 1.8763155423666869642e-19
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 14.99999542229031424 1 4.7818388965941336467e-06
WXEPOCH 55000.0000000000000000
WXFREQ_0001 0.00025100401605860257
WXSIN_0001 7.87630556749117e-06 1 6.395015583920052e-07
WXCOS_0001 -4.8313082227793976e-06 1 1.1334393962608818e-05
WXFREQ_0002 0.0005020080321172051
WXSIN_0002 4.446539271089234e-06 1 3.228766993754119e-07
WXCOS_0002 1.7077979616939418e-06 1 2.879560871355152e-06
WXFREQ_0003 0.0007530120481758077
WXSIN_0003 1.7770112498406102e-07 1 2.252346510942049e-07
WXCOS_0003 -1.2928260259488196e-06 1 1.3182172427722592e-06
WXFREQ_0004 0.0010040160642344103
WXSIN_0004 -1.0171136469853734e-07 1 1.7902516708664308e-07
WXCOS_0004 5.308400452612013e-07 1 7.755855549969352e-07
WXFREQ_0005 0.0012550200802930128
WXSIN_0005 -2.518130603083194e-07 1 1.56135754956151e-07
WXCOS_0005 -5.409773712695599e-07 1 5.297234117200345e-07
WXFREQ_0006 0.0015060240963516154
WXSIN_0006 8.625095536785448e-08 1 1.4505112475429255e-07
WXCOS_0006 7.329399661204629e-07 1 4.0288276288236733e-07
WXFREQ_0007 0.001757028112410218
WXSIN_0007 -9.402844238732055e-09 1 1.4491603784484472e-07
WXCOS_0007 -3.6015801922273094e-07 1 3.374973577233556e-07
WXFREQ_0008 0.0020080321284688205
WXSIN_0008 4.3530605907453027e-07 1 1.5632543539876843e-07
WXCOS_0008 3.1965112788241467e-07 1 3.124960757158423e-07
WXFREQ_0009 0.0022590361445274233
WXSIN_0009 -1.5907229664778663e-07 1 1.9671850448928628e-07
WXCOS_0009 -5.274669149296134e-07 1 3.3801688477128653e-07
WXFREQ_0010 0.0025100401605860257
WXSIN_0010 4.0976098097240536e-07 1 3.3951941702306133e-07
WXCOS_0010 4.491245513980958e-07 1 5.144165152890652e-07
WXFREQ_0011 0.0027610441766446284
WXSIN_0011 1.96657968682622e-06 1 2.7932237427074086e-06
WXCOS_0011 3.0788657961974026e-06 1 3.7025194919351555e-06
WXFREQ_0012 0.003012048192703231
WXSIN_0012 -1.5650982542769485e-07 1 2.039685783714564e-07
WXCOS_0012 -1.3133173689001158e-07 1 2.3655394309063525e-07
WXFREQ_0013 0.0032630522087618336
WXSIN_0013 2.1430417666281822e-07 1 9.584541770356817e-08
WXCOS_0013 6.11156768577608e-08 1 9.871441387099666e-08
WXFREQ_0014 0.003514056224820436
WXSIN_0014 -5.605581824932753e-08 1 6.133951846763445e-08
WXCOS_0014 -1.1351191995712074e-07 1 5.927287593503186e-08
WXFREQ_0015 0.0037650602408790387
WXSIN_0015 5.408446861509041e-08 1 4.703833366109986e-08
WXCOS_0015 1.4062989170927806e-07 1 4.426338467396346e-08
WXFREQ_0016 0.004016064256937641
WXSIN_0016 -9.222551201617452e-08 1 3.997369503084622e-08
WXCOS_0016 8.611987038331605e-08 1 3.813201166403714e-08
WXFREQ_0017 0.004267068272996243
WXSIN_0017 7.62054906960942e-08 1 3.630992236453892e-08
WXCOS_0017 -3.9524255161295953e-08 1 3.66180454601658e-08
WXFREQ_0018 0.004518072289054847
WXSIN_0018 -8.351371079657899e-08 1 3.5230684050886065e-08
WXCOS_0018 -1.1163604002476528e-07 1 3.466297355340797e-08
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.00035838892606131103 1 0.0009372785015328656
Estimating the spectral parameters from the WaveX fit.
[8]:
# Get the Fourier amplitudes and powers and their uncertainties.
idxs = np.array(m2.components["WaveX"].get_indices())
a = np.array([m2[f"WXSIN_{idx:04d}"].quantity.to_value("s") for idx in idxs])
da = np.array([m2[f"WXSIN_{idx:04d}"].uncertainty.to_value("s") for idx in idxs])
b = np.array([m2[f"WXCOS_{idx:04d}"].quantity.to_value("s") for idx in idxs])
db = np.array([m2[f"WXCOS_{idx:04d}"].uncertainty.to_value("s") for idx in idxs])
print(len(idxs))
P = (a**2 + b**2) / 2
dP = ((a * da) ** 2 + (b * db) ** 2) ** 0.5
f0 = (1 / Tspan).to_value(u.Hz)
fyr = (1 / u.year).to_value(u.Hz)
18
[9]:
# We can create a `PLRedNoise` model from the `WaveX` model.
# This will estimate the spectral parameters from the `WaveX`
# amplitudes.
m3 = plrednoise_from_wavex(m2)
print(m3)
# Created: 2026-06-16T15:30:48.900577
# PINT_version: 1.1.5+169.ge4c9571
# User: docs
# Host: build-33166788-project-85767-nanograv-pint
# OS: Linux-7.0.0-1004-aws-x86_64-with-glibc2.35
# Python: 3.11.14 (main, Apr 27 2026, 17:28:30) [GCC 11.4.0]
# Format: pint
# read_time: 2026-06-16T15:29:59.873849
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM3
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566498843
FINISH 56985.0000000464006366
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1923.9784834863046
CHI2R 0.9831264606470642
TRES 0.99161851118135482043
RAJ 4:59:59.99998969 1 0.00011766275985031068
DECJ 14:59:59.99093491 1 0.01196695415329324629
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999967901 1 5.71130312176685461e-13
F1 -9.9984421525715518555e-16 1 1.8763155423666869642e-19
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 14.99999542229031424 1 4.7818388965941336467e-06
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.00035838892606131103 1 0.0009372785015328656
TNREDAMP -12.805785754475949 0 0.08984350944880343
TNREDGAM 3.278850409118224 0 0.3479611181230549
TNREDC 18
[10]:
# Now let us plot the estimated spectrum with the injected
# spectrum.
plt.subplot(211)
plt.errorbar(
idxs * f0,
b * 1e6,
db * 1e6,
ls="",
marker="o",
label="$\\hat{a}_j$ (WXCOS)",
color="red",
)
plt.errorbar(
idxs * f0,
a * 1e6,
da * 1e6,
ls="",
marker="o",
label="$\\hat{b}_j$ (WXSIN)",
color="blue",
)
plt.axvline(fyr, color="black", ls="dotted")
plt.axhline(0, color="grey", ls="--")
plt.ylabel("Fourier coeffs ($\mu$s)")
plt.xscale("log")
plt.legend(fontsize=8)
plt.subplot(212)
plt.errorbar(
idxs * f0, P, dP, ls="", marker="o", label="Spectral power (PINT)", color="k"
)
P_inj = m.components["PLRedNoise"].get_noise_weights(t)[::2][:nharm_opt]
plt.plot(idxs * f0, P_inj, label="Injected Spectrum", color="r")
P_est = m3.components["PLRedNoise"].get_noise_weights(t)[::2][:nharm_opt]
print(len(idxs), len(P_est))
plt.plot(idxs * f0, P_est, label="Estimated Spectrum", color="b")
plt.xscale("log")
plt.yscale("log")
plt.ylabel("Spectral power (s$^2$)")
plt.xlabel("Frequency (Hz)")
plt.axvline(fyr, color="black", ls="dotted", label="1 yr$^{-1}$")
plt.legend()
18 18
[10]:
<matplotlib.legend.Legend at 0x7c48d4497390>
Note the outlier in the 1 year^-1 bin. This is caused by the covariance with RA and DEC, which introduce a delay with the same frequency.
DM noise fitting
Let us now do a similar kind of analysis for DM noise.
[11]:
par_sim = """
PSR SIM4
RAJ 05:00:00 1
DECJ 15:00:00 1
PEPOCH 55000
F0 100 1
F1 -1e-15 1
PHOFF 0 1
DM 15 1
TNDMAMP -13
TNDMGAM 3.5
TNDMC 30
TZRMJD 55000
TZRFRQ 1400
TZRSITE gbt
UNITS TDB
EPHEM DE440
CLOCK TT(BIPM2019)
"""
m = get_model(StringIO(par_sim))
[12]:
# Generate the simulated TOAs.
ntoas = 2000
toaerrs = np.random.uniform(0.5, 2.0, ntoas) * u.us
freqs = np.linspace(500, 1500, 8) * u.MHz
t = make_fake_toas_uniform(
startMJD=53001,
endMJD=57001,
ntoas=ntoas,
model=m,
freq=freqs,
obs="gbt",
error=toaerrs,
add_noise=True,
add_correlated_noise=True,
name="fake",
include_bipm=True,
multi_freqs_in_epoch=True,
)
[13]:
# Find the optimum number of harmonics by minimizing AIC.
m1 = deepcopy(m)
m1.remove_component("PLDMNoise")
m2 = deepcopy(m1)
nharm_opt, d_aics = find_optimal_nharms(m2, t, "DMWaveX", 30)
print("Optimum no of harmonics = ", nharm_opt)
Optimum no of harmonics = 29
[14]:
# The Y axis is plotted in log scale only for better visibility.
plt.scatter(list(range(len(d_aics))), d_aics + 1)
plt.axvline(nharm_opt, color="red", label="Optimum number of harmonics")
plt.axvline(
int(m.TNDMC.value), color="black", ls="--", label="Injected number of harmonics"
)
plt.xlabel("Number of harmonics")
plt.ylabel("AIC - AIC$_\\min{} + 1$")
plt.legend()
plt.yscale("log")
# plt.savefig("sim3-aic.pdf")
[15]:
# Now create a new model with the optimum number of
# harmonics
m2 = deepcopy(m1)
Tspan = t.get_mjds().max() - t.get_mjds().min()
dmwavex_setup(m2, T_span=Tspan, n_freqs=nharm_opt, freeze_params=False)
ftr = WLSFitter(t, m2)
ftr.fit_toas(maxiter=10)
m2 = ftr.model
print(m2)
# Created: 2026-06-16T15:31:47.195513
# PINT_version: 1.1.5+169.ge4c9571
# User: docs
# Host: build-33166788-project-85767-nanograv-pint
# OS: Linux-7.0.0-1004-aws-x86_64-with-glibc2.35
# Python: 3.11.14 (main, Apr 27 2026, 17:28:30) [GCC 11.4.0]
# Format: pint
# read_time: 2026-06-16T15:30:49.550593
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM4
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566265741
FINISH 56985.0000000458077662
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1997.094858291999
CHI2R 1.0320903660423766
TRES 1.0043218436103544272
RAJ 4:59:59.99999456 1 0.00000191632881776798
DECJ 15:00:00.00027756 1 0.00016158622528123995
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999999348 1 3.6880733219803857885e-14
F1 -1.0000009966135068385e-15 1 8.4011048442534754603e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 14.999991231380702635 1 4.916714383096330572e-06
DMWXEPOCH 55000.0000000000000000
DMWXFREQ_0001 0.00025100401605863834
DMWXSIN_0001 -0.004673438815593799 1 5.937931393007601e-06
DMWXCOS_0001 0.004162625664029142 1 6.783586918030218e-06
DMWXFREQ_0002 0.0005020080321172767
DMWXSIN_0002 5.737023470751262e-05 1 4.710568018929444e-06
DMWXCOS_0002 0.0011317224339904266 1 4.4170594319774205e-06
DMWXFREQ_0003 0.000753012048175915
DMWXSIN_0003 -0.0004290383487264485 1 4.431597151607152e-06
DMWXCOS_0003 0.00036770094395078195 1 4.2793764055045436e-06
DMWXFREQ_0004 0.0010040160642345534
DMWXSIN_0004 0.00020329135826252607 1 4.312648386483265e-06
DMWXCOS_0004 -3.533769681517434e-05 1 4.2620330263794345e-06
DMWXFREQ_0005 0.0012550200802931917
DMWXSIN_0005 -3.111871922953532e-05 1 4.3731223765629275e-06
DMWXCOS_0005 8.390837392550991e-05 1 4.173953145146236e-06
DMWXFREQ_0006 0.00150602409635183
DMWXSIN_0006 8.621203291518582e-05 1 4.335450947374546e-06
DMWXCOS_0006 0.00011572465472972941 1 4.187427368481524e-06
DMWXFREQ_0007 0.0017570281124104684
DMWXSIN_0007 -6.227800373772157e-06 1 4.185629587601254e-06
DMWXCOS_0007 -0.0001357061673552174 1 4.327411674345253e-06
DMWXFREQ_0008 0.0020080321284691068
DMWXSIN_0008 2.6546823309771e-05 1 4.253527442453624e-06
DMWXCOS_0008 1.0480710076258976e-05 1 4.272239946728468e-06
DMWXFREQ_0009 0.002259036144527745
DMWXSIN_0009 2.7077213784615694e-05 1 4.2975962930038194e-06
DMWXCOS_0009 1.396523335827454e-06 1 4.245835018522734e-06
DMWXFREQ_0010 0.0025100401605863834
DMWXSIN_0010 -1.725104633771263e-05 1 4.280935066963676e-06
DMWXCOS_0010 2.470835754829421e-05 1 4.321836035341292e-06
DMWXFREQ_0011 0.002761044176645022
DMWXSIN_0011 6.328467604167578e-05 1 7.093595924801682e-06
DMWXCOS_0011 1.3120354942562288e-05 1 6.740135327420645e-06
DMWXFREQ_0012 0.00301204819270366
DMWXSIN_0012 2.1433177716624473e-05 1 4.3147212172285275e-06
DMWXCOS_0012 -9.823089243969911e-06 1 4.239564635766524e-06
DMWXFREQ_0013 0.0032630522087622985
DMWXSIN_0013 2.2902284513509346e-05 1 4.270974597641839e-06
DMWXCOS_0013 3.941542225910975e-05 1 4.240534608339072e-06
DMWXFREQ_0014 0.003514056224820937
DMWXSIN_0014 -3.0252477446788357e-05 1 4.267597870506938e-06
DMWXCOS_0014 -2.5488454438297174e-05 1 4.227337039318785e-06
DMWXFREQ_0015 0.003765060240879575
DMWXSIN_0015 3.0815479857198805e-05 1 4.221059801942717e-06
DMWXCOS_0015 4.9710737405938e-06 1 4.28192144330884e-06
DMWXFREQ_0016 0.0040160642569382135
DMWXSIN_0016 -4.3706242488209015e-05 1 4.322644599862427e-06
DMWXCOS_0016 -4.157507790937137e-05 1 4.163882926360797e-06
DMWXFREQ_0017 0.004267068272996851
DMWXSIN_0017 -4.078917527991374e-05 1 4.1529818669144804e-06
DMWXCOS_0017 4.072755794441562e-07 1 4.3347552479160055e-06
DMWXFREQ_0018 0.00451807228905549
DMWXSIN_0018 2.076752746130703e-05 1 4.250477891089431e-06
DMWXCOS_0018 -1.701446912236452e-05 1 4.246327383829583e-06
DMWXFREQ_0019 0.004769076305114128
DMWXSIN_0019 -3.1701026058245616e-06 1 4.3148640080051684e-06
DMWXCOS_0019 -1.8441617734253816e-06 1 4.163974442004694e-06
DMWXFREQ_0020 0.005020080321172767
DMWXSIN_0020 1.4007723501690241e-05 1 4.2106038746245126e-06
DMWXCOS_0020 2.4394755365587845e-05 1 4.23989033342807e-06
DMWXFREQ_0021 0.005271084337231405
DMWXSIN_0021 4.077771165596948e-06 1 4.205055652523034e-06
DMWXCOS_0021 9.99918583529193e-07 1 4.255949696489197e-06
DMWXFREQ_0022 0.005522088353290044
DMWXSIN_0022 -1.380771970358398e-05 1 4.236822162651564e-06
DMWXCOS_0022 -1.6671415643071027e-06 1 4.2237084512418865e-06
DMWXFREQ_0023 0.0057730923693486815
DMWXSIN_0023 -4.457548084875104e-06 1 4.186018711102732e-06
DMWXCOS_0023 2.5208013869700126e-06 1 4.251969527670012e-06
DMWXFREQ_0024 0.00602409638540732
DMWXSIN_0024 9.107416817200981e-06 1 4.196930802999487e-06
DMWXCOS_0024 1.3985190283652823e-06 1 4.245663766172168e-06
DMWXFREQ_0025 0.006275100401465958
DMWXSIN_0025 7.626331539142819e-06 1 4.109473079143571e-06
DMWXCOS_0025 1.6571324357938042e-05 1 4.330191540014328e-06
DMWXFREQ_0026 0.006526104417524597
DMWXSIN_0026 3.4259960279564395e-06 1 4.213606386277551e-06
DMWXCOS_0026 3.957915071913053e-06 1 4.228607128745629e-06
DMWXFREQ_0027 0.006777108433583235
DMWXSIN_0027 1.1549554697443532e-05 1 4.194137861665957e-06
DMWXCOS_0027 -6.068077237099529e-06 1 4.2500757751283866e-06
DMWXFREQ_0028 0.007028112449641874
DMWXSIN_0028 8.928723579735839e-06 1 4.209308679359159e-06
DMWXCOS_0028 5.758201465044699e-06 1 4.236034390918995e-06
DMWXFREQ_0029 0.0072791164657005115
DMWXSIN_0029 1.04356732242846e-05 1 4.240786991793992e-06
DMWXCOS_0029 -2.020493319842389e-06 1 4.201299933880773e-06
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.0012249280178032544 1 5.580901321283589e-06
Estimating the spectral parameters from the DMWaveX fit.
[16]:
# Get the Fourier amplitudes and powers and their uncertainties.
# Note that the `DMWaveX` amplitudes have the units of DM.
# We multiply them by a constant factor to convert them to dimensions
# of time so that they are consistent with `PLDMNoise`.
scale = DMconst / (1400 * u.MHz) ** 2
idxs = np.array(m2.components["DMWaveX"].get_indices())
a = np.array(
[(scale * m2[f"DMWXSIN_{idx:04d}"].quantity).to_value("s") for idx in idxs]
)
da = np.array(
[(scale * m2[f"DMWXSIN_{idx:04d}"].uncertainty).to_value("s") for idx in idxs]
)
b = np.array(
[(scale * m2[f"DMWXCOS_{idx:04d}"].quantity).to_value("s") for idx in idxs]
)
db = np.array(
[(scale * m2[f"DMWXCOS_{idx:04d}"].uncertainty).to_value("s") for idx in idxs]
)
print(len(idxs))
P = (a**2 + b**2) / 2
dP = ((a * da) ** 2 + (b * db) ** 2) ** 0.5
f0 = (1 / Tspan).to_value(u.Hz)
fyr = (1 / u.year).to_value(u.Hz)
29
[17]:
# We can create a `PLDMNoise` model from the `DMWaveX` model.
# This will estimate the spectral parameters from the `DMWaveX`
# amplitudes.
m3 = pldmnoise_from_dmwavex(m2)
print(m3)
# Created: 2026-06-16T15:31:47.254288
# PINT_version: 1.1.5+169.ge4c9571
# User: docs
# Host: build-33166788-project-85767-nanograv-pint
# OS: Linux-7.0.0-1004-aws-x86_64-with-glibc2.35
# Python: 3.11.14 (main, Apr 27 2026, 17:28:30) [GCC 11.4.0]
# Format: pint
# read_time: 2026-06-16T15:30:49.550593
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM4
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566265741
FINISH 56985.0000000458077662
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1997.094858291999
CHI2R 1.0320903660423766
TRES 1.0043218436103544272
RAJ 4:59:59.99999456 1 0.00000191632881776798
DECJ 15:00:00.00027756 1 0.00016158622528123995
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999999348 1 3.6880733219803857885e-14
F1 -1.0000009966135068385e-15 1 8.4011048442534754603e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 14.999991231380702635 1 4.916714383096330572e-06
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.0012249280178032544 1 5.580901321283589e-06
TNDMAMP -13.026328536325252 0 0.04346738636666348
TNDMGAM 3.6434267837432 0 0.21157276694374783
TNDMC 29
[18]:
# Now let us plot the estimated spectrum with the injected
# spectrum.
plt.subplot(211)
plt.errorbar(
idxs * f0,
b * 1e6,
db * 1e6,
ls="",
marker="o",
label="$\\hat{a}_j$ (DMWXCOS)",
color="red",
)
plt.errorbar(
idxs * f0,
a * 1e6,
da * 1e6,
ls="",
marker="o",
label="$\\hat{b}_j$ (DMWXSIN)",
color="blue",
)
plt.axvline(fyr, color="black", ls="dotted")
plt.axhline(0, color="grey", ls="--")
plt.ylabel("Fourier coeffs ($\mu$s)")
plt.xscale("log")
plt.legend(fontsize=8)
plt.subplot(212)
plt.errorbar(
idxs * f0, P, dP, ls="", marker="o", label="Spectral power (PINT)", color="k"
)
P_inj = m.components["PLDMNoise"].get_noise_weights(t)[::2][:nharm_opt]
plt.plot(idxs * f0, P_inj, label="Injected Spectrum", color="r")
P_est = m3.components["PLDMNoise"].get_noise_weights(t)[::2][:nharm_opt]
print(len(idxs), len(P_est))
plt.plot(idxs * f0, P_est, label="Estimated Spectrum", color="b")
plt.xscale("log")
plt.yscale("log")
plt.ylabel("Spectral power (s$^2$)")
plt.xlabel("Frequency (Hz)")
plt.axvline(fyr, color="black", ls="dotted", label="1 yr$^{-1}$")
plt.legend()
29 29
[18]:
<matplotlib.legend.Legend at 0x7c48d09c6d90>
Chromatic noise fitting
Let us now do a similar kind of analysis for chromatic noise.
[19]:
par_sim = """
PSR SIM5
RAJ 05:00:00 1
DECJ 15:00:00 1
PEPOCH 55000
F0 100 1
F1 -1e-15 1
PHOFF 0 1
DM 15
CM 1.2 1
TNCHROMIDX 3.5
TNCHROMAMP -13
TNCHROMGAM 3.5
TNCHROMC 30
TZRMJD 55000
TZRFRQ 1400
TZRSITE gbt
UNITS TDB
EPHEM DE440
CLOCK TT(BIPM2019)
"""
m = get_model(StringIO(par_sim))
[20]:
# Generate the simulated TOAs.
ntoas = 2000
toaerrs = np.random.uniform(0.5, 2.0, ntoas) * u.us
freqs = np.linspace(500, 1500, 8) * u.MHz
t = make_fake_toas_uniform(
startMJD=53001,
endMJD=57001,
ntoas=ntoas,
model=m,
freq=freqs,
obs="gbt",
error=toaerrs,
add_noise=True,
add_correlated_noise=True,
name="fake",
include_bipm=True,
multi_freqs_in_epoch=True,
)
[21]:
# Find the optimum number of harmonics by minimizing AIC.
m1 = deepcopy(m)
m1.remove_component("PLChromNoise")
m2 = deepcopy(m1)
nharm_opt = m.TNCHROMC.value
[22]:
# Now create a new model with the optimum number of
# harmonics
m2 = deepcopy(m1)
Tspan = t.get_mjds().max() - t.get_mjds().min()
cmwavex_setup(m2, T_span=Tspan, n_freqs=nharm_opt, freeze_params=False)
ftr = WLSFitter(t, m2)
ftr.fit_toas(maxiter=10)
m2 = ftr.model
print(m2)
# Created: 2026-06-16T15:31:57.123331
# PINT_version: 1.1.5+169.ge4c9571
# User: docs
# Host: build-33166788-project-85767-nanograv-pint
# OS: Linux-7.0.0-1004-aws-x86_64-with-glibc2.35
# Python: 3.11.14 (main, Apr 27 2026, 17:28:30) [GCC 11.4.0]
# Format: pint
# read_time: 2026-06-16T15:31:47.919521
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM5
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567779977
FINISH 56985.0000000491110185
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1933.5821732208865
CHI2R 1.0003011760066665
TRES 0.9731343536272374595
RAJ 4:59:59.99999888 1 0.00000141582133384566
DECJ 15:00:00.00019916 1 0.00012299125695230159
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.00000000000000988 1 2.842218334601495247e-14
F1 -9.999996318535180839e-16 1 6.3778723022545978236e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.0
CM 1.2275655284862623061 1 0.05033914643101146391
TNCHROMIDX 3.5
CMWXEPOCH 55000.0000000000000000
CMWXFREQ_0001 0.0002510040160584398
CMWXSIN_0001 -156.47881162298444 1 0.06784235519959746
CMWXCOS_0001 -209.2730293264276 1 0.06892134880046212
CMWXFREQ_0002 0.0005020080321168797
CMWXSIN_0002 -92.97425377479061 1 0.061483302822058356
CMWXCOS_0002 -8.794130655808134 1 0.05650175556551221
CMWXFREQ_0003 0.0007530120481753195
CMWXSIN_0003 -38.37968152001221 1 0.05805873749430889
CMWXCOS_0003 -11.000204782712641 1 0.05805292685014816
CMWXFREQ_0004 0.0010040160642337593
CMWXSIN_0004 -1.7827190783400626 1 0.059169304567122785
CMWXCOS_0004 -34.8134911448648 1 0.05648990240558088
CMWXFREQ_0005 0.0012550200802921992
CMWXSIN_0005 7.412224572417315 1 0.05811005213789977
CMWXCOS_0005 -12.220696571872764 1 0.05721208161077406
CMWXFREQ_0006 0.001506024096350639
CMWXSIN_0006 0.08132739258565867 1 0.05740933994353317
CMWXCOS_0006 -8.928788618854632 1 0.0580328924624548
CMWXFREQ_0007 0.001757028112409079
CMWXSIN_0007 3.215197747786588 1 0.05724526771704209
CMWXCOS_0007 3.669145894745538 1 0.05802175546094752
CMWXFREQ_0008 0.0020080321284675186
CMWXSIN_0008 1.4164326834944465 1 0.058113657746037696
CMWXCOS_0008 -8.000349399667199 1 0.0571653591485541
CMWXFREQ_0009 0.0022590361445259588
CMWXSIN_0009 -0.8117958671157368 1 0.057904485277933966
CMWXCOS_0009 2.65699064009519 1 0.05742448749260342
CMWXFREQ_0010 0.0025100401605843985
CMWXSIN_0010 -3.943746428855493 1 0.05733617364787144
CMWXCOS_0010 -3.2027716742945347 1 0.05839129470160487
CMWXFREQ_0011 0.002761044176642838
CMWXSIN_0011 -3.8881489973977823 1 0.06942188653089774
CMWXCOS_0011 2.539247837203363 1 0.07339374463795488
CMWXFREQ_0012 0.003012048192701278
CMWXSIN_0012 -0.13121826064744335 1 0.05960426470491224
CMWXCOS_0012 -3.7137116977896922 1 0.05568081426281196
CMWXFREQ_0013 0.003263052208759718
CMWXSIN_0013 -3.020096120247841 1 0.05766021026651007
CMWXCOS_0013 0.43398348261873 1 0.05772456512426281
CMWXFREQ_0014 0.003514056224818158
CMWXSIN_0014 0.548344778443811 1 0.05614613473492575
CMWXCOS_0014 1.1293692637182606 1 0.05920218756183869
CMWXFREQ_0015 0.0037650602408765975
CMWXSIN_0015 1.010405689338949 1 0.057239061292953405
CMWXCOS_0015 -1.5063958501169494 1 0.05719679427756715
CMWXFREQ_0016 0.004016064256935037
CMWXSIN_0016 -2.991601590388431 1 0.05936361001927001
CMWXCOS_0016 -1.0565150548783475 1 0.055021387869380926
CMWXFREQ_0017 0.004267068272993477
CMWXSIN_0017 -0.8662084714802362 1 0.05766984650555683
CMWXCOS_0017 0.02180045009922303 1 0.05679223117103709
CMWXFREQ_0018 0.0045180722890519175
CMWXSIN_0018 1.7089221883598649 1 0.05523059018386513
CMWXCOS_0018 -0.09202013728433725 1 0.05897161394151831
CMWXFREQ_0019 0.004769076305110357
CMWXSIN_0019 0.1777490638726145 1 0.057945137940994346
CMWXCOS_0019 -0.6848623534829872 1 0.05617302443119135
CMWXFREQ_0020 0.005020080321168797
CMWXSIN_0020 0.36475072174082973 1 0.0583401258952611
CMWXCOS_0020 0.013806313372791777 1 0.0559011857410082
CMWXFREQ_0021 0.005271084337227236
CMWXSIN_0021 -0.5865496769839305 1 0.05771374903297803
CMWXCOS_0021 1.014351909118882 1 0.05642683873666322
CMWXFREQ_0022 0.005522088353285676
CMWXSIN_0022 -1.4573324109824202 1 0.05705018392544127
CMWXCOS_0022 0.1513029979858186 1 0.05707869748874025
CMWXFREQ_0023 0.0057730923693441166
CMWXSIN_0023 0.6351027200005497 1 0.05689709104164802
CMWXCOS_0023 0.45238690451339203 1 0.057239922428756175
CMWXFREQ_0024 0.006024096385402556
CMWXSIN_0024 0.37596538246353095 1 0.05716629909089215
CMWXCOS_0024 -1.039383926725432 1 0.05689790950944832
CMWXFREQ_0025 0.006275100401460996
CMWXSIN_0025 -0.5528388454076073 1 0.0578499535845621
CMWXCOS_0025 1.6419337605593427 1 0.056070343698456586
CMWXFREQ_0026 0.006526104417519436
CMWXSIN_0026 0.42049369206271187 1 0.057441282197353484
CMWXCOS_0026 -0.327966211416106 1 0.05644061516823244
CMWXFREQ_0027 0.006777108433577875
CMWXSIN_0027 0.1556049735163054 1 0.05753060410909398
CMWXCOS_0027 -0.4043050730262996 1 0.056325769515768495
CMWXFREQ_0028 0.007028112449636316
CMWXSIN_0028 0.6409345796036154 1 0.05731111572812722
CMWXCOS_0028 0.4204223251072528 1 0.05667493488733421
CMWXFREQ_0029 0.007279116465694756
CMWXSIN_0029 -0.17363995806569427 1 0.05798260158321058
CMWXCOS_0029 -0.11558631018383 1 0.055751614523670605
CMWXFREQ_0030 0.007530120481753195
CMWXSIN_0030 0.1815994260094562 1 0.05593687229020186
CMWXCOS_0030 0.022656871606959447 1 0.057901798593078306
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.0011675865366946173 1 4.141406912581703e-06
Estimating the spectral parameters from the CMWaveX fit.
[23]:
# Get the Fourier amplitudes and powers and their uncertainties.
# Note that the `CMWaveX` amplitudes have the units of pc/cm^3/MHz^2.
# We multiply them by a constant factor to convert them to dimensions
# of time so that they are consistent with `PLChromNoise`.
scale = DMconst / 1400**m.TNCHROMIDX.value
idxs = np.array(m2.components["CMWaveX"].get_indices())
a = np.array(
[(scale * m2[f"CMWXSIN_{idx:04d}"].quantity).to_value("s") for idx in idxs]
)
da = np.array(
[(scale * m2[f"CMWXSIN_{idx:04d}"].uncertainty).to_value("s") for idx in idxs]
)
b = np.array(
[(scale * m2[f"CMWXCOS_{idx:04d}"].quantity).to_value("s") for idx in idxs]
)
db = np.array(
[(scale * m2[f"CMWXCOS_{idx:04d}"].uncertainty).to_value("s") for idx in idxs]
)
print(len(idxs))
P = (a**2 + b**2) / 2
dP = ((a * da) ** 2 + (b * db) ** 2) ** 0.5
f0 = (1 / Tspan).to_value(u.Hz)
fyr = (1 / u.year).to_value(u.Hz)
30
[24]:
# We can create a `PLChromNoise` model from the `CMWaveX` model.
# This will estimate the spectral parameters from the `CMWaveX`
# amplitudes.
m3 = plchromnoise_from_cmwavex(m2)
print(m3)
# Created: 2026-06-16T15:31:57.183490
# PINT_version: 1.1.5+169.ge4c9571
# User: docs
# Host: build-33166788-project-85767-nanograv-pint
# OS: Linux-7.0.0-1004-aws-x86_64-with-glibc2.35
# Python: 3.11.14 (main, Apr 27 2026, 17:28:30) [GCC 11.4.0]
# Format: pint
# read_time: 2026-06-16T15:31:47.919521
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM5
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567779977
FINISH 56985.0000000491110185
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1933.5821732208865
CHI2R 1.0003011760066665
TRES 0.9731343536272374595
RAJ 4:59:59.99999888 1 0.00000141582133384566
DECJ 15:00:00.00019916 1 0.00012299125695230159
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.00000000000000988 1 2.842218334601495247e-14
F1 -9.999996318535180839e-16 1 6.3778723022545978236e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.0
CM 1.2275655284862623061 1 0.05033914643101146391
TNCHROMIDX 3.5
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.0011675865366946173 1 4.141406912581703e-06
TNCHROMAMP -12.897628510889847 0 0.040055750428270226
TNCHROMGAM 3.66359423787523 0 0.22487030143936632
TNCHROMC 30
[25]:
# Now let us plot the estimated spectrum with the injected
# spectrum.
plt.subplot(211)
plt.errorbar(
idxs * f0,
b * 1e6,
db * 1e6,
ls="",
marker="o",
label="$\\hat{a}_j$ (CMWXCOS)",
color="red",
)
plt.errorbar(
idxs * f0,
a * 1e6,
da * 1e6,
ls="",
marker="o",
label="$\\hat{b}_j$ (CMWXSIN)",
color="blue",
)
plt.axvline(fyr, color="black", ls="dotted")
plt.axhline(0, color="grey", ls="--")
plt.ylabel("Fourier coeffs ($\mu$s)")
plt.xscale("log")
plt.legend(fontsize=8)
plt.subplot(212)
plt.errorbar(
idxs * f0, P, dP, ls="", marker="o", label="Spectral power (PINT)", color="k"
)
P_inj = m.components["PLChromNoise"].get_noise_weights(t)[::2]
plt.plot(idxs * f0, P_inj, label="Injected Spectrum", color="r")
P_est = m3.components["PLChromNoise"].get_noise_weights(t)[::2]
print(len(idxs), len(P_est))
plt.plot(idxs * f0, P_est, label="Estimated Spectrum", color="b")
plt.xscale("log")
plt.yscale("log")
plt.ylabel("Spectral power (s$^2$)")
plt.xlabel("Frequency (Hz)")
plt.axvline(fyr, color="black", ls="dotted", label="1 yr$^{-1}$")
plt.legend()
30 30
[25]:
<matplotlib.legend.Legend at 0x7c48c4586e90>
[ ]: