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-06-16T15:29:56.733011
# PINT_version: 0+untagged.351.g83be0a2
# User: docs
# Host: build-33166772-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:07.267595
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM3
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566809954
FINISH 56985.0000000464353819
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1951.584733541817
CHI2R 1.0002997096575177
TRES 0.9790470088140139426
RAJ 4:59:59.99976877 1 0.00012290824797668560
DECJ 15:00:00.00502328 1 0.01237636593940938737
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999890516 1 5.655111173044867538e-13
F1 -1.0000117911287594679e-15 1 2.0031634978341586121e-19
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 14.99999676768375974 1 4.7673521682351223763e-06
WXEPOCH 55000.0000000000000000
WXFREQ_0001 0.00025100401605860213
WXSIN_0001 -1.3208778359357957e-05 1 6.346781293497607e-07
WXCOS_0001 -2.283770015983773e-06 1 1.2096647326138968e-05
WXFREQ_0002 0.0005020080321172043
WXSIN_0002 1.65714671200821e-06 1 3.20675138397766e-07
WXCOS_0002 -2.9266561467908325e-06 1 3.072107254843646e-06
WXFREQ_0003 0.0007530120481758064
WXSIN_0003 -1.1637066168258049e-06 1 2.240664764541234e-07
WXCOS_0003 9.878580909681147e-07 1 1.404725446888021e-06
WXFREQ_0004 0.0010040160642344085
WXSIN_0004 9.343014366915392e-07 1 1.7831855201555246e-07
WXCOS_0004 -1.3686819408741409e-07 1 8.256106690606549e-07
WXFREQ_0005 0.0012550200802930107
WXSIN_0005 1.0595839474520305e-08 1 1.5495359636243588e-07
WXCOS_0005 4.5294791101782764e-07 1 5.621852151372626e-07
WXFREQ_0006 0.0015060240963516128
WXSIN_0006 3.972574658788551e-07 1 1.4505290718738397e-07
WXCOS_0006 -8.350127200896726e-08 1 4.271870589634109e-07
WXFREQ_0007 0.0017570281124102147
WXSIN_0007 -5.795887173075903e-07 1 1.4465077168237787e-07
WXCOS_0007 -4.899848230142127e-08 1 3.554636531179947e-07
WXFREQ_0008 0.002008032128468817
WXSIN_0008 8.868962186947401e-08 1 1.5813400367602905e-07
WXCOS_0008 2.4097051339166314e-07 1 3.2847502494889654e-07
WXFREQ_0009 0.002259036144527419
WXSIN_0009 -5.142973177415298e-07 1 1.982737841105345e-07
WXCOS_0009 2.0816322078920427e-07 1 3.536193302112138e-07
WXFREQ_0010 0.0025100401605860213
WXSIN_0010 7.777715191633866e-07 1 3.4537099264683965e-07
WXCOS_0010 1.7319377167798162e-08 1 5.335498456415592e-07
WXFREQ_0011 0.0027610441766446232
WXSIN_0011 6.779051843385078e-06 1 2.86585994494125e-06
WXCOS_0011 -4.6961555016139953e-07 1 3.806068310972492e-06
WXFREQ_0012 0.0030120481927032256
WXSIN_0012 -6.033629157598381e-07 1 2.1017080508434213e-07
WXCOS_0012 -6.730873248825573e-09 1 2.4137895036359155e-07
WXFREQ_0013 0.0032630522087618275
WXSIN_0013 2.8013223679470676e-07 1 9.86794272096636e-08
WXCOS_0013 -6.000297549738066e-08 1 1.0066036058583734e-07
WXFREQ_0014 0.0035140562248204294
WXSIN_0014 3.3598990609374996e-08 1 6.26566581632472e-08
WXCOS_0014 -3.647691544293056e-08 1 5.932698833959353e-08
WXFREQ_0015 0.0037650602408790318
WXSIN_0015 5.459297804466251e-08 1 4.749098801931556e-08
WXCOS_0015 7.928872576197998e-08 1 4.517299969431721e-08
WXFREQ_0016 0.004016064256937634
WXSIN_0016 -9.805914255199913e-08 1 4.0162355261782676e-08
WXCOS_0016 1.3021144024650625e-09 1 3.966715493904549e-08
WXFREQ_0017 0.004267068272996236
WXSIN_0017 -3.9779510652563426e-08 1 3.664516291783137e-08
WXCOS_0017 -1.0643862087359721e-07 1 3.672049715198773e-08
WXFREQ_0018 0.004518072289054838
WXSIN_0018 -6.900856633784327e-08 1 3.451758682410306e-08
WXCOS_0018 1.965917271345588e-07 1 3.537802466844712e-08
WXFREQ_0019 0.00476907630511344
WXSIN_0019 2.4123446148180242e-08 1 3.338876309013701e-08
WXCOS_0019 -5.648524071266121e-09 1 3.4735993006996573e-08
WXFREQ_0020 0.005020080321172043
WXSIN_0020 -9.102638009777712e-08 1 3.266655229953428e-08
WXCOS_0020 6.619383196128747e-08 1 3.443605511952207e-08
WXFREQ_0021 0.005271084337230644
WXSIN_0021 8.908257281473297e-08 1 3.233339825682494e-08
WXCOS_0021 3.036139006960431e-08 1 3.46710488597315e-08
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.0004023028099447161 1 0.001000026270911382
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-06-16T15:29:56.783648
# PINT_version: 0+untagged.351.g83be0a2
# User: docs
# Host: build-33166772-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:07.267595
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM3
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566809954
FINISH 56985.0000000464353819
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1951.584733541817
CHI2R 1.0002997096575177
TRES 0.9790470088140139426
RAJ 4:59:59.99976877 1 0.00012290824797668560
DECJ 15:00:00.00502328 1 0.01237636593940938737
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999890516 1 5.655111173044867538e-13
F1 -1.0000117911287594679e-15 1 2.0031634978341586121e-19
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 14.99999676768375974 1 4.7673521682351223763e-06
TNREDAMP -12.72195348774838 0 0.07259657553335375
TNREDGAM 3.2116791198723833 0 0.29738214897124854
TNREDC 21
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.0004023028099447161 1 0.001000026270911382
[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 0x7cf8578db750>
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 = 28
[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:30:55.367947
# PINT_version: 0+untagged.351.g83be0a2
# User: docs
# Host: build-33166772-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:57.585086
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM4
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566998032
FINISH 56985.0000000464155671
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1933.9475280542529
CHI2R 0.9984241239309515
TRES 0.9914951298652044531
RAJ 5:00:00.00000085 1 0.00000186584718916673
DECJ 14:59:59.99998385 1 0.00016225945899174527
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999996574 1 3.6043675695916718273e-14
F1 -9.999977280390689541e-16 1 8.2584617722719201526e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.000007916902145865 1 4.8766877304889739086e-06
DMWXEPOCH 55000.0000000000000000
DMWXFREQ_0001 0.00025100401605860484
DMWXSIN_0001 0.003944851450022517 1 5.9969390873376055e-06
DMWXCOS_0001 -0.0011508949521212297 1 6.6315144517716605e-06
DMWXFREQ_0002 0.0005020080321172097
DMWXSIN_0002 0.00028192285924450366 1 4.758738258912527e-06
DMWXCOS_0002 -0.0004409386690095899 1 4.4786190291747786e-06
DMWXFREQ_0003 0.0007530120481758146
DMWXSIN_0003 -0.00010411660824808581 1 4.4735883354710236e-06
DMWXCOS_0003 0.0006586044631343503 1 4.2861620625498755e-06
DMWXFREQ_0004 0.0010040160642344194
DMWXSIN_0004 8.718103657410036e-05 1 4.371288193998216e-06
DMWXCOS_0004 -4.5240183897777344e-05 1 4.282425565693022e-06
DMWXFREQ_0005 0.0012550200802930243
DMWXSIN_0005 6.152517357532041e-05 1 4.360127604071423e-06
DMWXCOS_0005 -5.104361914622239e-05 1 4.2314142607487054e-06
DMWXFREQ_0006 0.0015060240963516293
DMWXSIN_0006 8.449220389770276e-05 1 4.3219400549233395e-06
DMWXCOS_0006 -6.354952896931756e-05 1 4.244614309168858e-06
DMWXFREQ_0007 0.001757028112410234
DMWXSIN_0007 -9.331860500814688e-05 1 4.2757035939982175e-06
DMWXCOS_0007 0.00010085259529356617 1 4.278598450221941e-06
DMWXFREQ_0008 0.0020080321284688387
DMWXSIN_0008 2.6570777627852083e-05 1 4.313918201916955e-06
DMWXCOS_0008 -0.0001054799977991553 1 4.230459774194885e-06
DMWXFREQ_0009 0.0022590361445274437
DMWXSIN_0009 -5.4768508652516836e-05 1 4.279559043272739e-06
DMWXCOS_0009 5.4307996524547384e-05 1 4.294273655352419e-06
DMWXFREQ_0010 0.0025100401605860486
DMWXSIN_0010 0.00011301192389737522 1 4.316568677954768e-06
DMWXCOS_0010 1.484929697755365e-05 1 4.3026593939583745e-06
DMWXFREQ_0011 0.0027610441766446536
DMWXSIN_0011 7.5794319679485814e-06 1 6.7547859171800925e-06
DMWXCOS_0011 -1.4266517542903668e-05 1 6.930842023974544e-06
DMWXFREQ_0012 0.0030120481927032585
DMWXSIN_0012 -1.4691141110190099e-05 1 4.24220794433152e-06
DMWXCOS_0012 6.268841313613307e-05 1 4.367198761883052e-06
DMWXFREQ_0013 0.003263052208761863
DMWXSIN_0013 1.8187292461372242e-05 1 4.217998271207673e-06
DMWXCOS_0013 6.9151727033496375e-06 1 4.334414915288688e-06
DMWXFREQ_0014 0.003514056224820468
DMWXSIN_0014 8.46944642466001e-07 1 4.352075802126238e-06
DMWXCOS_0014 -1.3110108688312778e-05 1 4.189448468433252e-06
DMWXFREQ_0015 0.003765060240879073
DMWXSIN_0015 -2.427662955382817e-05 1 4.24944858335635e-06
DMWXCOS_0015 1.3624048812311805e-05 1 4.296231389289082e-06
DMWXFREQ_0016 0.0040160642569376775
DMWXSIN_0016 -2.7399316017184847e-06 1 4.295959475641541e-06
DMWXCOS_0016 -7.1975551750687054e-06 1 4.251702798600754e-06
DMWXFREQ_0017 0.004267068272996282
DMWXSIN_0017 6.400770853156522e-06 1 4.284849597626741e-06
DMWXCOS_0017 -2.3877074096900198e-05 1 4.261605296842162e-06
DMWXFREQ_0018 0.004518072289054887
DMWXSIN_0018 1.6258412133636523e-05 1 4.304421057215738e-06
DMWXCOS_0018 3.733456869322243e-06 1 4.246824929982875e-06
DMWXFREQ_0019 0.004769076305113492
DMWXSIN_0019 -2.82332165404682e-05 1 4.368194773361173e-06
DMWXCOS_0019 4.674784550847678e-07 1 4.172590718971024e-06
DMWXFREQ_0020 0.005020080321172097
DMWXSIN_0020 -1.3585683863118855e-05 1 4.274795700131409e-06
DMWXCOS_0020 -2.194662847940208e-05 1 4.266935400167186e-06
DMWXFREQ_0021 0.005271084337230702
DMWXSIN_0021 6.708884837625576e-06 1 4.280996256929388e-06
DMWXCOS_0021 -5.67435993504744e-06 1 4.267563450646817e-06
DMWXFREQ_0022 0.005522088353289307
DMWXSIN_0022 1.794147409123793e-06 1 4.248384384770392e-06
DMWXCOS_0022 1.432089847977288e-05 1 4.301762209406312e-06
DMWXFREQ_0023 0.005773092369347912
DMWXSIN_0023 -4.147701445232653e-06 1 4.237627012664696e-06
DMWXCOS_0023 2.0162651040834566e-05 1 4.312334804706647e-06
DMWXFREQ_0024 0.006024096385406517
DMWXSIN_0024 1.2278693909423492e-06 1 4.26469103366907e-06
DMWXCOS_0024 -7.318237107650921e-06 1 4.292040152742763e-06
DMWXFREQ_0025 0.006275100401465122
DMWXSIN_0025 -9.481377120308524e-06 1 4.342551298640034e-06
DMWXCOS_0025 1.652233926515226e-05 1 4.214320935125864e-06
DMWXFREQ_0026 0.006526104417523726
DMWXSIN_0026 -1.6794427653943042e-05 1 4.360968791238735e-06
DMWXCOS_0026 3.40607843967227e-06 1 4.201650724934321e-06
DMWXFREQ_0027 0.006777108433582331
DMWXSIN_0027 7.531120462994882e-06 1 4.312545277894153e-06
DMWXCOS_0027 -8.936894325951539e-06 1 4.251164478994344e-06
DMWXFREQ_0028 0.007028112449640936
DMWXSIN_0028 1.2322592864142744e-05 1 4.319524588354342e-06
DMWXCOS_0028 -8.587651770846054e-06 1 4.234210102559955e-06
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.00020805363836361095 1 5.520289782651061e-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)
28
[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:30:55.427768
# PINT_version: 0+untagged.351.g83be0a2
# User: docs
# Host: build-33166772-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:57.585086
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM4
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566998032
FINISH 56985.0000000464155671
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1933.9475280542529
CHI2R 0.9984241239309515
TRES 0.9914951298652044531
RAJ 5:00:00.00000085 1 0.00000186584718916673
DECJ 14:59:59.99998385 1 0.00016225945899174527
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999996574 1 3.6043675695916718273e-14
F1 -9.999977280390689541e-16 1 8.2584617722719201526e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.000007916902145865 1 4.8766877304889739086e-06
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.00020805363836361095 1 5.520289782651061e-06
TNDMAMP -13.037576171429485 0 0.04418788463048408
TNDMGAM 3.418341263535047 0 0.22859000979200736
TNDMC 28
[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()
28 28
[18]:
<matplotlib.legend.Legend at 0x7cf8558d6f90>
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:05.649697
# PINT_version: 0+untagged.351.g83be0a2
# User: docs
# Host: build-33166772-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:56.098951
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM5
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566747801
FINISH 56985.0000000437054398
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 2082.2372385015005
CHI2R 1.077204986291516
TRES 1.0113495951876170086
RAJ 5:00:00.00000010 1 0.00000140239870330976
DECJ 15:00:00.00007314 1 0.00012120361705024044
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999993615 1 2.7908091376039748657e-14
F1 -1.0000002869016239469e-15 1 6.268905856420818565e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.0
CM 1.1552543165493932638 1 0.04928779845880076199
TNCHROMIDX 3.5
CMWXEPOCH 55000.0000000000000000
CMWXFREQ_0001 0.00025100401605877354
CMWXSIN_0001 24.2528583323045 1 0.06601844102615165
CMWXCOS_0001 63.64895961178953 1 0.06664163239159845
CMWXFREQ_0002 0.0005020080321175471
CMWXSIN_0002 -29.452666182658415 1 0.05795624293358263
CMWXCOS_0002 -78.0782555241342 1 0.05676005993901208
CMWXFREQ_0003 0.0007530120481763206
CMWXSIN_0003 5.297815740858238 1 0.05603499321903633
CMWXCOS_0003 -1.658215950214373 1 0.056671793443033205
CMWXFREQ_0004 0.0010040160642350942
CMWXSIN_0004 8.431095455837257 1 0.05862440116113204
CMWXCOS_0004 -4.569180279471016 1 0.053581570737380256
CMWXFREQ_0005 0.0012550200802938678
CMWXSIN_0005 -12.198701191257134 1 0.055285743921599566
CMWXCOS_0005 -3.1994387535170516 1 0.056686205828646016
CMWXFREQ_0006 0.0015060240963526413
CMWXSIN_0006 -1.238018601456843 1 0.05677385884649177
CMWXCOS_0006 0.35188758250670055 1 0.05513812190865324
CMWXFREQ_0007 0.001757028112411415
CMWXSIN_0007 1.5349890075730517 1 0.0556709771674127
CMWXCOS_0007 5.21854628848217 1 0.05627430052826464
CMWXFREQ_0008 0.0020080321284701884
CMWXSIN_0008 -5.483089242724355 1 0.054770463465851744
CMWXCOS_0008 -0.6089108961704854 1 0.05718545094141403
CMWXFREQ_0009 0.002259036144528962
CMWXSIN_0009 6.7498047583712895 1 0.05715143334500793
CMWXCOS_0009 4.592167687690279 1 0.05455931576314862
CMWXFREQ_0010 0.0025100401605877357
CMWXSIN_0010 0.3259412308906067 1 0.05527035708535987
CMWXCOS_0010 -5.2826370700815595 1 0.05678140365451549
CMWXFREQ_0011 0.0027610441766465093
CMWXSIN_0011 3.3410007608916144 1 0.06727860979720093
CMWXCOS_0011 -0.8314999327717051 1 0.07098383153758876
CMWXFREQ_0012 0.0030120481927052825
CMWXSIN_0012 0.35240944193751766 1 0.05662995346029895
CMWXCOS_0012 -4.603849297126597 1 0.0550003733506768
CMWXFREQ_0013 0.003263052208764056
CMWXSIN_0013 1.9688031013426548 1 0.05599283944623908
CMWXCOS_0013 -0.10994310452251944 1 0.05535471319322811
CMWXFREQ_0014 0.00351405622482283
CMWXSIN_0014 0.5037075645225052 1 0.057530943037966936
CMWXCOS_0014 0.4615470307965615 1 0.053717706117246924
CMWXFREQ_0015 0.0037650602408816035
CMWXSIN_0015 0.008541193728753017 1 0.055575515849292205
CMWXCOS_0015 0.7099093626745742 1 0.05568589321636468
CMWXFREQ_0016 0.004016064256940377
CMWXSIN_0016 1.6232379093946698 1 0.05646125018111799
CMWXCOS_0016 -1.362062634581622 1 0.054947614043596354
CMWXFREQ_0017 0.004267068272999151
CMWXSIN_0017 -0.2961821681503682 1 0.05657422296017559
CMWXCOS_0017 0.8094412044079063 1 0.05482160849799616
CMWXFREQ_0018 0.004518072289057924
CMWXSIN_0018 1.4057588192419872 1 0.05536810863670256
CMWXCOS_0018 -1.7298814159542122 1 0.055875990409268314
CMWXFREQ_0019 0.004769076305116697
CMWXSIN_0019 1.3202165732574291 1 0.056240646492222025
CMWXCOS_0019 -0.18828145703270355 1 0.055011928411504124
CMWXFREQ_0020 0.005020080321175471
CMWXSIN_0020 0.7932057747712411 1 0.05643760266266435
CMWXCOS_0020 1.232407902592499 1 0.054854815705462655
CMWXFREQ_0021 0.0052710843372342445
CMWXSIN_0021 -0.0923580897165119 1 0.05512205022971997
CMWXCOS_0021 0.4925247625342274 1 0.05587636811502565
CMWXFREQ_0022 0.005522088353293019
CMWXSIN_0022 -0.0451921901026705 1 0.054707816093424524
CMWXCOS_0022 0.08632987987032902 1 0.05627606119440757
CMWXFREQ_0023 0.005773092369351792
CMWXSIN_0023 -0.40826765470578513 1 0.05405140987575455
CMWXCOS_0023 -0.5149663485832364 1 0.05664136292597702
CMWXFREQ_0024 0.006024096385410565
CMWXSIN_0024 0.12439273647812511 1 0.05474852228326578
CMWXCOS_0024 0.14465332865330596 1 0.05589497331199401
CMWXFREQ_0025 0.006275100401469339
CMWXSIN_0025 -0.5529580012733823 1 0.05425523778341763
CMWXCOS_0025 -0.418659698911582 1 0.05633724984478436
CMWXFREQ_0026 0.006526104417528112
CMWXSIN_0026 -0.08580838429663527 1 0.055757925404464934
CMWXCOS_0026 -0.18053869781082277 1 0.054846964945206585
CMWXFREQ_0027 0.0067771084335868865
CMWXSIN_0027 -0.09651067938986281 1 0.05601295656779403
CMWXCOS_0027 -0.4223089462309307 1 0.05492514501951058
CMWXFREQ_0028 0.00702811244964566
CMWXSIN_0028 -0.0567434548717356 1 0.055300709478640124
CMWXCOS_0028 0.5354998117876048 1 0.05557654341131873
CMWXFREQ_0029 0.007279116465704433
CMWXSIN_0029 0.5682789775200966 1 0.05480640847292998
CMWXCOS_0029 0.0012723055244423282 1 0.05598652043455457
CMWXFREQ_0030 0.007530120481763207
CMWXSIN_0030 0.08295058402457779 1 0.05569367157706305
CMWXCOS_0030 0.057714957393351064 1 0.05504845682215121
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -9.647095965909566e-05 1 4.192052616838693e-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:05.709960
# PINT_version: 0+untagged.351.g83be0a2
# User: docs
# Host: build-33166772-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:56.098951
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM5
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566747801
FINISH 56985.0000000437054398
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 2082.2372385015005
CHI2R 1.077204986291516
TRES 1.0113495951876170086
RAJ 5:00:00.00000010 1 0.00000140239870330976
DECJ 15:00:00.00007314 1 0.00012120361705024044
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999993615 1 2.7908091376039748657e-14
F1 -1.0000002869016239469e-15 1 6.268905856420818565e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.0
CM 1.1552543165493932638 1 0.04928779845880076199
TNCHROMIDX 3.5
TNCHROMAMP -13.004209507724797 0 0.04020851691493547
TNCHROMGAM 3.5483570745337563 0 0.2793595537815093
TNCHROMC 30
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -9.647095965909566e-05 1 4.192052616838693e-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 0x7cf857570390>
[ ]: