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 = 21
[5]:
print(np.argmin(d_aics))
21
[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-08-10T10:10:48.126942
# PINT_version: 1.1.6
# User: docs
# Host: build-33990685-project-85767-nanograv-pint
# OS: Linux-7.0.0-1004-aws-x86_64-with-glibc2.35
# Python: 3.11.15 (main, Jun 25 2026, 19:09:59) [GCC 11.4.0]
# Format: pint
# read_time: 2026-08-10T10:10:00.608061
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
# ell1h_shapiro: full
PSR SIM3
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566839353
FINISH 56985.0000000464228241
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1842.306875317326
CHI2R 0.944288506057061
TRES 0.9633547039608704
RAJ 5:00:00.00004011 1 0.00012933186343048518
DECJ 15:00:00.00883405 1 0.01234923138653131955
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.000000000000335516 1 6.0219670820867707056e-13
F1 -1.0002865970237783904e-15 1 2.0037728494403568477e-19
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 14.999997816322307586 1 4.9499092287333770575e-06
WXEPOCH 55000.0000000000000000
WXFREQ_0001 0.0002510040160586035
WXSIN_0001 -1.0125079012386734e-05 1 6.747327635173275e-07
WXCOS_0001 1.7716877736487565e-05 1 1.2097361619979654e-05
WXFREQ_0002 0.000502008032117207
WXSIN_0002 -3.1305271364755854e-06 1 3.4199132866948384e-07
WXCOS_0002 -3.822213154436823e-06 1 3.0712454342797366e-06
WXFREQ_0003 0.0007530120481758105
WXSIN_0003 1.2470340555131879e-06 1 2.3765661163362963e-07
WXCOS_0003 3.325256989925778e-06 1 1.4044176148770118e-06
WXFREQ_0004 0.001004016064234414
WXSIN_0004 -2.3052463809270375e-07 1 1.9004462700352514e-07
WXCOS_0004 -8.754165013126249e-07 1 8.2309099906726e-07
WXFREQ_0005 0.0012550200802930174
WXSIN_0005 3.115364293546183e-07 1 1.6508194399934906e-07
WXCOS_0005 1.1565885639500766e-06 1 5.606449785344647e-07
WXFREQ_0006 0.001506024096351621
WXSIN_0006 -1.647213419572226e-07 1 1.5399183296291006e-07
WXCOS_0006 -4.0092959666786525e-07 1 4.250053420482762e-07
WXFREQ_0007 0.0017570281124102245
WXSIN_0007 -1.4681660371050598e-07 1 1.5381665999081215e-07
WXCOS_0007 6.340285747012862e-07 1 3.5325722564633857e-07
WXFREQ_0008 0.002008032128468828
WXSIN_0008 -1.0150637044903032e-07 1 1.6740252246330775e-07
WXCOS_0008 -4.800433559253762e-07 1 3.2511941432318265e-07
WXFREQ_0009 0.0022590361445274315
WXSIN_0009 1.8904765476417068e-07 1 2.0894607794529595e-07
WXCOS_0009 4.057207379357608e-07 1 3.495195263942994e-07
WXFREQ_0010 0.0025100401605860348
WXSIN_0010 -1.8832179126502472e-07 1 3.649635289479369e-07
WXCOS_0010 -6.115456532815659e-07 1 5.275680569209263e-07
WXFREQ_0011 0.0027610441766446384
WXSIN_0011 -2.858489991340295e-06 1 3.009761786727167e-06
WXCOS_0011 -3.054541392482574e-06 1 3.7627645828047955e-06
WXFREQ_0012 0.003012048192703242
WXSIN_0012 1.5100839612473117e-07 1 2.195611496409546e-07
WXCOS_0012 1.596701960356686e-07 1 2.3973904635791444e-07
WXFREQ_0013 0.0032630522087618453
WXSIN_0013 1.2400770705250664e-08 1 1.0256200407642977e-07
WXCOS_0013 -7.864310285597518e-08 1 9.947580408181371e-08
WXFREQ_0014 0.003514056224820449
WXSIN_0014 3.5789500943003796e-08 1 6.451600752182331e-08
WXCOS_0014 4.130241423790594e-08 1 5.954857758113654e-08
WXFREQ_0015 0.0037650602408790526
WXSIN_0015 1.0203523220235725e-08 1 4.9871146177864125e-08
WXCOS_0015 -1.927659269008461e-08 1 4.569931399607159e-08
WXFREQ_0016 0.004016064256937656
WXSIN_0016 -1.0380882892057015e-07 1 4.172370565617514e-08
WXCOS_0016 -2.9418699466232256e-08 1 3.978163026420516e-08
WXFREQ_0017 0.004267068272996259
WXSIN_0017 -4.827082438755104e-08 1 3.8045868568631256e-08
WXCOS_0017 8.032782115931235e-08 1 3.682178490199846e-08
WXFREQ_0018 0.004518072289054863
WXSIN_0018 5.558334109585222e-09 1 3.559125318693464e-08
WXCOS_0018 3.1165135543025857e-09 1 3.5896319032417385e-08
WXFREQ_0019 0.004769076305113466
WXSIN_0019 -4.671956247473819e-08 1 3.406751682054381e-08
WXCOS_0019 1.0311268196426344e-07 1 3.53556773990748e-08
WXFREQ_0020 0.0050200803211720695
WXSIN_0020 -6.262638838107317e-09 1 3.323083285196919e-08
WXCOS_0020 -1.1314176279299631e-07 1 3.480333512617898e-08
WXFREQ_0021 0.005271084337230674
WXSIN_0021 -1.3638323069847597e-07 1 3.31201245187844e-08
WXCOS_0021 4.98047840814985e-09 1 3.491533001211159e-08
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.00033639049386909573 1 0.000999886940607425
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)
21
[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-08-10T10:10:48.171777
# PINT_version: 1.1.6
# User: docs
# Host: build-33990685-project-85767-nanograv-pint
# OS: Linux-7.0.0-1004-aws-x86_64-with-glibc2.35
# Python: 3.11.15 (main, Jun 25 2026, 19:09:59) [GCC 11.4.0]
# Format: pint
# read_time: 2026-08-10T10:10:00.608061
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
# ell1h_shapiro: full
PSR SIM3
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566839353
FINISH 56985.0000000464228241
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1842.306875317326
CHI2R 0.944288506057061
TRES 0.9633547039608704
RAJ 5:00:00.00004011 1 0.00012933186343048518
DECJ 15:00:00.00883405 1 0.01234923138653131955
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.000000000000335516 1 6.0219670820867707056e-13
F1 -1.0002865970237783904e-15 1 2.0037728494403568477e-19
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 14.999997816322307586 1 4.9499092287333770575e-06
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.00033639049386909573 1 0.000999886940607425
TNREDAMP -12.848295384013536 0 0.08896577400964638
TNREDGAM 3.5226650209790593 0 0.36106570183671055
TNREDC 21
[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()
21 21
[10]:
<matplotlib.legend.Legend at 0x7cb65b02d3d0>
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 = 30
[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-08-10T10:11:44.259027
# PINT_version: 1.1.6
# User: docs
# Host: build-33990685-project-85767-nanograv-pint
# OS: Linux-7.0.0-1004-aws-x86_64-with-glibc2.35
# Python: 3.11.15 (main, Jun 25 2026, 19:09:59) [GCC 11.4.0]
# Format: pint
# read_time: 2026-08-10T10:10:48.700700
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
# ell1h_shapiro: full
PSR SIM4
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567161227
FINISH 56985.0000000468601273
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1773.0158156244734
CHI2R 0.9172352900281807
TRES 0.9482653554749005
RAJ 5:00:00.00000088 1 0.00000188711617330221
DECJ 14:59:59.99987556 1 0.00016302553344418170
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999996768 1 3.568705882408330807e-14
F1 -1.000000148808350431e-15 1 8.37975029792795233e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.000005342538324451 1 4.8420015269153692395e-06
DMWXEPOCH 55000.0000000000000000
DMWXFREQ_0001 0.00025100401605857785
DMWXSIN_0001 -0.0038079236773060406 1 5.810743262895173e-06
DMWXCOS_0001 -0.0030565162174329068 1 6.6941051345300455e-06
DMWXFREQ_0002 0.0005020080321171557
DMWXSIN_0002 0.0003246801438623925 1 4.657117935767851e-06
DMWXCOS_0002 -0.0002098390406724936 1 4.43826381771098e-06
DMWXFREQ_0003 0.0007530120481757334
DMWXSIN_0003 -0.00030103521010844866 1 4.416650318796112e-06
DMWXCOS_0003 0.0003190830315197773 1 4.2357061835729424e-06
DMWXFREQ_0004 0.0010040160642343114
DMWXSIN_0004 -0.0001904543248709882 1 4.243884698086943e-06
DMWXCOS_0004 -2.8675022290491754e-05 1 4.319843451886819e-06
DMWXFREQ_0005 0.001255020080292889
DMWXSIN_0005 7.439938236255072e-05 1 4.334514837055442e-06
DMWXCOS_0005 0.00010517698594713213 1 4.156714226277986e-06
DMWXFREQ_0006 0.0015060240963514669
DMWXSIN_0006 0.000242814269543571 1 4.348316808316206e-06
DMWXCOS_0006 -0.00022154429402135596 1 4.140921516913525e-06
DMWXFREQ_0007 0.0017570281124100447
DMWXSIN_0007 4.311994134669624e-05 1 4.192584228653275e-06
DMWXCOS_0007 7.011258512808586e-05 1 4.3023256959127915e-06
DMWXFREQ_0008 0.0020080321284686228
DMWXSIN_0008 6.375090715720695e-07 1 4.213925729706536e-06
DMWXCOS_0008 0.00013397449998138182 1 4.268843547479331e-06
DMWXFREQ_0009 0.0022590361445272004
DMWXSIN_0009 -4.1027263154724674e-05 1 4.290838754415032e-06
DMWXCOS_0009 1.962330435549146e-05 1 4.231372418360385e-06
DMWXFREQ_0010 0.002510040160585778
DMWXSIN_0010 -4.367634440876183e-05 1 4.281982463991482e-06
DMWXCOS_0010 -3.2884476276930484e-05 1 4.274793170026374e-06
DMWXFREQ_0011 0.002761044176644356
DMWXSIN_0011 0.00013320024378035778 1 6.817287237801698e-06
DMWXCOS_0011 5.4247213645385526e-05 1 6.806970825438207e-06
DMWXFREQ_0012 0.0030120481927029337
DMWXSIN_0012 9.813346009663466e-06 1 4.278512760145588e-06
DMWXCOS_0012 6.766807947313797e-05 1 4.257134582927263e-06
DMWXFREQ_0013 0.0032630522087615118
DMWXSIN_0013 -3.6448178243252915e-05 1 4.3061448816498985e-06
DMWXCOS_0013 3.997282358251552e-05 1 4.170051099924139e-06
DMWXFREQ_0014 0.0035140562248200894
DMWXSIN_0014 1.5462706335345484e-05 1 4.279303550759608e-06
DMWXCOS_0014 3.0493688076823184e-05 1 4.20107346354003e-06
DMWXFREQ_0015 0.0037650602408786675
DMWXSIN_0015 -2.555304092561154e-05 1 4.210476899061622e-06
DMWXCOS_0015 -1.0755109907846927e-05 1 4.262738598387476e-06
DMWXFREQ_0016 0.0040160642569372455
DMWXSIN_0016 1.2103242947259025e-05 1 4.218419490979351e-06
DMWXCOS_0016 -3.294008823733775e-05 1 4.251619772716002e-06
DMWXFREQ_0017 0.004267068272995823
DMWXSIN_0017 -2.182230758862097e-05 1 4.261074060256383e-06
DMWXCOS_0017 1.623111478321896e-05 1 4.215298952104221e-06
DMWXFREQ_0018 0.004518072289054401
DMWXSIN_0018 1.6607733074296317e-06 1 4.179838474165949e-06
DMWXCOS_0018 -3.4457086679132695e-05 1 4.304295994506785e-06
DMWXFREQ_0019 0.004769076305112979
DMWXSIN_0019 4.324412912061833e-06 1 4.102921252963899e-06
DMWXCOS_0019 1.939551653035899e-05 1 4.36509509106343e-06
DMWXFREQ_0020 0.005020080321171556
DMWXSIN_0020 -2.969564079702465e-06 1 4.20191105448612e-06
DMWXCOS_0020 2.18879359971371e-05 1 4.2801576669947365e-06
DMWXFREQ_0021 0.005271084337230134
DMWXSIN_0021 2.178119387130424e-05 1 4.251675698906013e-06
DMWXCOS_0021 1.6060616988125625e-05 1 4.2369234152579935e-06
DMWXFREQ_0022 0.005522088353288712
DMWXSIN_0022 -2.224382519774949e-05 1 4.240308793929227e-06
DMWXCOS_0022 3.7841529813021664e-06 1 4.247296090400192e-06
DMWXFREQ_0023 0.00577309236934729
DMWXSIN_0023 9.722817566722268e-06 1 4.251876260637681e-06
DMWXCOS_0023 -2.752779725001777e-06 1 4.221844006863506e-06
DMWXFREQ_0024 0.006024096385405867
DMWXSIN_0024 4.687790022099295e-06 1 4.370406735538767e-06
DMWXCOS_0024 1.6912115810043765e-05 1 4.099190791729044e-06
DMWXFREQ_0025 0.0062751004014644455
DMWXSIN_0025 -3.5230146337844144e-06 1 4.1638836822189055e-06
DMWXCOS_0025 9.699683314658086e-06 1 4.314867828147191e-06
DMWXFREQ_0026 0.0065261044175230236
DMWXSIN_0026 1.3270317514195725e-05 1 4.148233346629196e-06
DMWXCOS_0026 8.139650172604014e-06 1 4.329728837983825e-06
DMWXFREQ_0027 0.006777108433581601
DMWXSIN_0027 -1.1986682687120947e-05 1 4.197086533020083e-06
DMWXCOS_0027 9.814547596057298e-06 1 4.2838526499666385e-06
DMWXFREQ_0028 0.007028112449640179
DMWXSIN_0028 5.754920450954809e-06 1 4.306403002923543e-06
DMWXCOS_0028 6.909695192667511e-06 1 4.183321289835546e-06
DMWXFREQ_0029 0.007279116465698757
DMWXSIN_0029 1.3879323104929954e-05 1 4.2490989190128636e-06
DMWXCOS_0029 -6.476069870633514e-06 1 4.250153910103605e-06
DMWXFREQ_0030 0.007530120481757335
DMWXSIN_0030 7.665015106301059e-06 1 4.2313755568746815e-06
DMWXCOS_0030 1.4434284486107814e-05 1 4.26622806814973e-06
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.0005653409228365521 1 5.676054692702931e-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)
30
[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-08-10T10:11:44.317483
# PINT_version: 1.1.6
# User: docs
# Host: build-33990685-project-85767-nanograv-pint
# OS: Linux-7.0.0-1004-aws-x86_64-with-glibc2.35
# Python: 3.11.15 (main, Jun 25 2026, 19:09:59) [GCC 11.4.0]
# Format: pint
# read_time: 2026-08-10T10:10:48.700700
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
# ell1h_shapiro: full
PSR SIM4
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567161227
FINISH 56985.0000000468601273
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1773.0158156244734
CHI2R 0.9172352900281807
TRES 0.9482653554749005
RAJ 5:00:00.00000088 1 0.00000188711617330221
DECJ 14:59:59.99987556 1 0.00016302553344418170
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999996768 1 3.568705882408330807e-14
F1 -1.000000148808350431e-15 1 8.37975029792795233e-22
PEPOCH 55000.0000000000000000
TNDMAMP -12.940750353397455 0 0.04172638398317115
TNDMGAM 3.403261645948323 0 0.23734454795829635
TNDMC 30
PLANET_SHAPIRO N
DM 15.000005342538324451 1 4.8420015269153692395e-06
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.0005653409228365521 1 5.676054692702931e-06
[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()
30 30
[18]:
<matplotlib.legend.Legend at 0x7cb65b003590>
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-08-10T10:11:53.810640
# PINT_version: 1.1.6
# User: docs
# Host: build-33990685-project-85767-nanograv-pint
# OS: Linux-7.0.0-1004-aws-x86_64-with-glibc2.35
# Python: 3.11.15 (main, Jun 25 2026, 19:09:59) [GCC 11.4.0]
# Format: pint
# read_time: 2026-08-10T10:11:44.763689
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
# ell1h_shapiro: full
PSR SIM5
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566785186
FINISH 56985.0000000454472917
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1861.2912037530832
CHI2R 0.962902847259743
TRES 0.9641518902005014
RAJ 5:00:00.00000115 1 0.00000142870787399426
DECJ 14:59:59.99991140 1 0.00012395731002708311
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.999999999999999244 1 2.7847054506167165792e-14
F1 -1.000000740207259087e-15 1 6.337177503176894825e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.0
CM 1.1581582454906943621 1 0.052297256504414343925
TNCHROMIDX 3.5
CMWXEPOCH 55000.0000000000000000
CMWXFREQ_0001 0.0002510040160586645
CMWXSIN_0001 -70.1823303257482 1 0.06963314547575705
CMWXCOS_0001 46.70336380361093 1 0.07157592456305598
CMWXFREQ_0002 0.000502008032117329
CMWXSIN_0002 5.114977223671063 1 0.06318392116265154
CMWXCOS_0002 17.27793318847422 1 0.05924038611434245
CMWXFREQ_0003 0.0007530120481759934
CMWXSIN_0003 -4.896030356923332 1 0.06040329974811901
CMWXCOS_0003 18.19568730037808 1 0.0597548413000732
CMWXFREQ_0004 0.001004016064234658
CMWXSIN_0004 -3.2300610245261367 1 0.06000076088458665
CMWXCOS_0004 -12.942627449546837 1 0.05973810663940792
CMWXFREQ_0005 0.0012550200802933223
CMWXSIN_0005 -1.973100281394335 1 0.059522427902953236
CMWXCOS_0005 -21.85876906600081 1 0.05963507024957328
CMWXFREQ_0006 0.0015060240963519868
CMWXSIN_0006 11.552080078873958 1 0.05885933005124603
CMWXCOS_0006 -6.005726038229253 1 0.060312299114850886
CMWXFREQ_0007 0.0017570281124106512
CMWXSIN_0007 -2.395113006847731 1 0.059347358938853456
CMWXCOS_0007 -1.6546168382490765 1 0.059817815693410485
CMWXFREQ_0008 0.002008032128469316
CMWXSIN_0008 -0.9973616362755124 1 0.06092519757761954
CMWXCOS_0008 -5.131916658216853 1 0.058044959755987884
CMWXFREQ_0009 0.00225903614452798
CMWXSIN_0009 -3.5050075805845817 1 0.06172106037629442
CMWXCOS_0009 2.773175783134212 1 0.05711304984311804
CMWXFREQ_0010 0.0025100401605866445
CMWXSIN_0010 1.0569036027100898 1 0.060420470634067164
CMWXCOS_0010 -4.207698441270572 1 0.05851218402000446
CMWXFREQ_0011 0.002761044176645309
CMWXSIN_0011 -2.010778516910985 1 0.07172059708908303
CMWXCOS_0011 -1.10242546933938 1 0.07534934553474795
CMWXFREQ_0012 0.0030120481927039737
CMWXSIN_0012 -2.60788868783269 1 0.060693541904580046
CMWXCOS_0012 -2.221002059983309 1 0.058181126176079824
CMWXFREQ_0013 0.003263052208762638
CMWXSIN_0013 -0.9419047765362768 1 0.05975065703704542
CMWXCOS_0013 0.06187202542082926 1 0.05858531785935004
CMWXFREQ_0014 0.0035140562248213024
CMWXSIN_0014 0.9484741522326854 1 0.059796632134138064
CMWXCOS_0014 -1.8051822550719738 1 0.058472296139283494
CMWXFREQ_0015 0.0037650602408799668
CMWXSIN_0015 -0.805451838208905 1 0.06003179352147423
CMWXCOS_0015 1.7210852877989404 1 0.058444949265424706
CMWXFREQ_0016 0.004016064256938632
CMWXSIN_0016 -0.14226836820992292 1 0.05852379683967763
CMWXCOS_0016 -0.6985352393008305 1 0.05996714127389343
CMWXFREQ_0017 0.0042670682729972955
CMWXSIN_0017 0.15527256471590645 1 0.05854134560056935
CMWXCOS_0017 -0.8174532454668375 1 0.059973877657177675
CMWXFREQ_0018 0.00451807228905596
CMWXSIN_0018 -0.18315621486485198 1 0.05784336865798854
CMWXCOS_0018 0.5924181736580001 1 0.06057868392312835
CMWXFREQ_0019 0.004769076305114625
CMWXSIN_0019 1.759763471894801 1 0.05968889526302939
CMWXCOS_0019 -0.8607902650388735 1 0.058848836415688284
CMWXFREQ_0020 0.005020080321173289
CMWXSIN_0020 0.8842182215073995 1 0.059147376140564464
CMWXCOS_0020 0.1652492644278259 1 0.059467909894479844
CMWXFREQ_0021 0.005271084337231954
CMWXSIN_0021 -0.39647730795831515 1 0.05753467845556319
CMWXCOS_0021 -0.8687119453136597 1 0.06111469279836891
CMWXFREQ_0022 0.005522088353290618
CMWXSIN_0022 0.21158581798352719 1 0.05976180242129498
CMWXCOS_0022 0.09682579677499832 1 0.05876059758909464
CMWXFREQ_0023 0.0057730923693492826
CMWXSIN_0023 -1.0817696986482588 1 0.057954665430151354
CMWXCOS_0023 -0.46704928955869646 1 0.06056340295915835
CMWXFREQ_0024 0.006024096385407947
CMWXSIN_0024 -0.5080869098371633 1 0.059778133445736
CMWXCOS_0024 0.5593380018649725 1 0.05868559941884537
CMWXFREQ_0025 0.006275100401466611
CMWXSIN_0025 -0.05719750255544224 1 0.06000513826147847
CMWXCOS_0025 0.7152212106328439 1 0.05864196147932282
CMWXFREQ_0026 0.006526104417525276
CMWXSIN_0026 -0.2008377406942896 1 0.05940644349790113
CMWXCOS_0026 0.3190371440618515 1 0.05935646444181074
CMWXFREQ_0027 0.00677710843358394
CMWXSIN_0027 -0.15385206776929303 1 0.059032404949413586
CMWXCOS_0027 -0.08365855078242622 1 0.05965444014429133
CMWXFREQ_0028 0.007028112449642605
CMWXSIN_0028 -0.5046448620240755 1 0.061670617776188645
CMWXCOS_0028 -1.1028082560197965 1 0.056840377991122626
CMWXFREQ_0029 0.00727911646570127
CMWXSIN_0029 -0.8669039528572352 1 0.0594994467336283
CMWXCOS_0029 -0.5402104735916821 1 0.05920990080552273
CMWXFREQ_0030 0.0075301204817599336
CMWXSIN_0030 -0.49978301136348335 1 0.058202471423731535
CMWXCOS_0030 0.28827390952041115 1 0.06076822141675765
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.0001033021121693765 1 4.273918973644669e-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-08-10T10:11:53.871096
# PINT_version: 1.1.6
# User: docs
# Host: build-33990685-project-85767-nanograv-pint
# OS: Linux-7.0.0-1004-aws-x86_64-with-glibc2.35
# Python: 3.11.15 (main, Jun 25 2026, 19:09:59) [GCC 11.4.0]
# Format: pint
# read_time: 2026-08-10T10:11:44.763689
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
# ell1h_shapiro: full
PSR SIM5
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566785186
FINISH 56985.0000000454472917
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1861.2912037530832
CHI2R 0.962902847259743
TRES 0.9641518902005014
RAJ 5:00:00.00000115 1 0.00000142870787399426
DECJ 14:59:59.99991140 1 0.00012395731002708311
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.999999999999999244 1 2.7847054506167165792e-14
F1 -1.000000740207259087e-15 1 6.337177503176894825e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.0
CM 1.1581582454906943621 1 0.052297256504414343925
TNCHROMIDX 3.5
TNCHROMAMP -12.990224296660577 0 0.040129860937228215
TNCHROMGAM 3.1522694089859997 0 0.2450868749712295
TNCHROMC 30
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.0001033021121693765 1 4.273918973644669e-06
[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 0x7cb658628250>
[ ]: