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 = 17
[5]:
print(np.argmin(d_aics))
17
[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-11T13:46:51.618767
# PINT_version: 0+untagged.358.g3534316
# User: docs
# Host: build-33094749-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-11T13:46:01.760322
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM3
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566199421
FINISH 56985.0000000463431250
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 2137.309933985424
CHI2R 1.0910208953473322
TRES 1.0198563575311198425
RAJ 5:00:00.00027243 1 0.00010863394973928482
DECJ 14:59:59.98326786 1 0.01147777773965706280
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.000000000001199145 1 5.2199504969153217947e-13
F1 -9.997224121729288738e-16 1 1.7908484068095687825e-19
PEPOCH 55000.0000000000000000
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.000724703363393334 1 0.0008949538911816551
PLANET_SHAPIRO N
DM 14.999996681118964236 1 4.838938253647725367e-06
WXEPOCH 55000.0000000000000000
WXFREQ_0001 0.0002510040160586044
WXSIN_0001 1.988845985244696e-06 1 5.845261521303025e-07
WXCOS_0001 -8.656395761718925e-06 1 1.0818352928929756e-05
WXFREQ_0002 0.0005020080321172088
WXSIN_0002 5.010078014297111e-07 1 2.9578608392829814e-07
WXCOS_0002 3.1188902004808038e-06 1 2.748801125285131e-06
WXFREQ_0003 0.0007530120481758132
WXSIN_0003 1.3365465505747103e-07 1 2.0633265970772213e-07
WXCOS_0003 -1.4248950033534733e-06 1 1.26005799261847e-06
WXFREQ_0004 0.0010040160642344176
WXSIN_0004 2.2104283511014212e-07 1 1.641983002704933e-07
WXCOS_0004 6.461636713214089e-07 1 7.415200316456581e-07
WXFREQ_0005 0.001255020080293022
WXSIN_0005 7.45188169904801e-07 1 1.4356195007177694e-07
WXCOS_0005 -1.1694910132766195e-06 1 5.070822267999735e-07
WXFREQ_0006 0.0015060240963516265
WXSIN_0006 -1.861691304015269e-07 1 1.336122930824571e-07
WXCOS_0006 6.906271146374542e-07 1 3.848000854186677e-07
WXFREQ_0007 0.0017570281124102308
WXSIN_0007 3.176665844297749e-07 1 1.3260385314823517e-07
WXCOS_0007 -4.767446665653848e-07 1 3.234161051169419e-07
WXFREQ_0008 0.0020080321284688353
WXSIN_0008 -3.3867055428821445e-07 1 1.4448664292631638e-07
WXCOS_0008 4.270605135671799e-07 1 3.0159594051378095e-07
WXFREQ_0009 0.00225903614452744
WXSIN_0009 4.736671777060554e-07 1 1.78707853275465e-07
WXCOS_0009 -3.546388132822399e-07 1 3.2591827968071506e-07
WXFREQ_0010 0.002510040160586044
WXSIN_0010 -6.767528126938449e-07 1 3.077343904492864e-07
WXCOS_0010 7.515176261980837e-07 1 4.942679161610601e-07
WXFREQ_0011 0.0027610441766446484
WXSIN_0011 -5.995781895396237e-06 1 2.5285953445832744e-06
WXCOS_0011 4.159272007662487e-06 1 3.558537606887618e-06
WXFREQ_0012 0.003012048192703253
WXSIN_0012 3.5300819009555507e-07 1 1.8506347809556867e-07
WXCOS_0012 -2.0726476007427142e-07 1 2.2659662113179767e-07
WXFREQ_0013 0.0032630522087618574
WXSIN_0013 -3.5966856888969666e-07 1 8.63838840920329e-08
WXCOS_0013 9.077027986545232e-08 1 9.407905584109688e-08
WXFREQ_0014 0.0035140562248204615
WXSIN_0014 1.1491152256803388e-07 1 5.7266459521624464e-08
WXCOS_0014 6.229741355654915e-08 1 5.67486051428223e-08
WXFREQ_0015 0.003765060240879066
WXSIN_0015 -1.5926227334837086e-07 1 4.37399427388115e-08
WXCOS_0015 8.40156356516682e-08 1 4.324998019218788e-08
WXFREQ_0016 0.0040160642569376705
WXSIN_0016 2.2601163520247756e-08 1 3.785061689042232e-08
WXCOS_0016 2.390170162048386e-08 1 3.719423620162383e-08
WXFREQ_0017 0.004267068272996275
WXSIN_0017 -2.7249371539861618e-08 1 3.470223027368066e-08
WXCOS_0017 1.0988989179171135e-07 1 3.540766934383782e-08
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)
17
[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-11T13:46:51.665854
# PINT_version: 0+untagged.358.g3534316
# User: docs
# Host: build-33094749-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-11T13:46:01.760322
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM3
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566199421
FINISH 56985.0000000463431250
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 2137.309933985424
CHI2R 1.0910208953473322
TRES 1.0198563575311198425
RAJ 5:00:00.00027243 1 0.00010863394973928482
DECJ 14:59:59.98326786 1 0.01147777773965706280
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.000000000001199145 1 5.2199504969153217947e-13
F1 -9.997224121729288738e-16 1 1.7908484068095687825e-19
PEPOCH 55000.0000000000000000
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.000724703363393334 1 0.0008949538911816551
TNREDAMP -12.679373707226073 0 0.09873787982826318
TNREDGAM 2.0138551240046967 0 0.5540771690302679
TNREDC 17
PLANET_SHAPIRO N
DM 14.999996681118964236 1 4.838938253647725367e-06
[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()
17 17
[10]:
<matplotlib.legend.Legend at 0x713539ec4090>
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-06-11T13:47:50.287440
# PINT_version: 0+untagged.358.g3534316
# User: docs
# Host: build-33094749-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-11T13:46:52.386764
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM4
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566857755
FINISH 56985.0000000463691088
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1915.3788043255029
CHI2R 0.990884016722971
TRES 0.986608999272166773
RAJ 5:00:00.00000268 1 0.00000190066703634089
DECJ 15:00:00.00005030 1 0.00016116959119829902
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999998855 1 3.6214478305399477814e-14
F1 -1.000000295184137167e-15 1 8.470114338831499001e-22
PEPOCH 55000.0000000000000000
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.0002580513947990552 1 5.644264446900314e-06
PLANET_SHAPIRO N
DM 14.9999955975508036075 1 5.011542322785246341e-06
DMWXEPOCH 55000.0000000000000000
DMWXFREQ_0001 0.0002510040160586067
DMWXSIN_0001 -0.0010367744431310747 1 6.0680070038837885e-06
DMWXCOS_0001 0.0010379105790988348 1 6.9647295687457745e-06
DMWXFREQ_0002 0.0005020080321172134
DMWXSIN_0002 -0.0005930877276584442 1 4.793643664182902e-06
DMWXCOS_0002 4.4675563803074385e-05 1 4.6524988941670945e-06
DMWXFREQ_0003 0.0007530120481758201
DMWXSIN_0003 -0.0002357782876967481 1 4.666656732682045e-06
DMWXCOS_0003 0.0002047297140146733 1 4.349869843048785e-06
DMWXFREQ_0004 0.0010040160642344267
DMWXSIN_0004 -8.409196003227821e-05 1 4.6517308229716095e-06
DMWXCOS_0004 0.00028867775590641633 1 4.263369070017021e-06
DMWXFREQ_0005 0.0012550200802930334
DMWXSIN_0005 3.103956886620585e-05 1 4.4651398661883355e-06
DMWXCOS_0005 -0.00042999164517669485 1 4.396167303522242e-06
DMWXFREQ_0006 0.0015060240963516401
DMWXSIN_0006 8.109678489196654e-05 1 4.465633949016279e-06
DMWXCOS_0006 2.0874191591450272e-05 1 4.347894567698819e-06
DMWXFREQ_0007 0.0017570281124102468
DMWXSIN_0007 -5.3869150624798836e-05 1 4.38240317978187e-06
DMWXCOS_0007 3.364069095808586e-05 1 4.4113638077040245e-06
DMWXFREQ_0008 0.0020080321284688535
DMWXSIN_0008 -4.888537870476354e-05 1 4.477916754172091e-06
DMWXCOS_0008 -4.960075148077053e-06 1 4.273263892366744e-06
DMWXFREQ_0009 0.00225903614452746
DMWXSIN_0009 3.686527347016429e-05 1 4.254247362118613e-06
DMWXCOS_0009 1.2702059713231734e-05 1 4.485974732865493e-06
DMWXFREQ_0010 0.002510040160586067
DMWXSIN_0010 -3.9984144768808895e-05 1 4.324888065781149e-06
DMWXCOS_0010 6.942934665112012e-05 1 4.468379122797948e-06
DMWXFREQ_0011 0.0027610441766446735
DMWXSIN_0011 7.816633745149227e-05 1 7.148942165666428e-06
DMWXCOS_0011 3.351332002224217e-08 1 6.883847799007546e-06
DMWXFREQ_0012 0.0030120481927032802
DMWXSIN_0012 2.7859067637117046e-05 1 4.402247359908801e-06
DMWXCOS_0012 -2.9134057953379664e-05 1 4.3555415856058924e-06
DMWXFREQ_0013 0.003263052208761887
DMWXSIN_0013 -5.906885442786023e-06 1 4.355694286843811e-06
DMWXCOS_0013 -1.5484175366817797e-05 1 4.347532825818152e-06
DMWXFREQ_0014 0.0035140562248204936
DMWXSIN_0014 1.4010366570279612e-05 1 4.368538100623782e-06
DMWXCOS_0014 -3.367543809067054e-05 1 4.35045263608568e-06
DMWXFREQ_0015 0.0037650602408791007
DMWXSIN_0015 3.179896921219616e-05 1 4.335547694757826e-06
DMWXCOS_0015 -1.344104764043655e-05 1 4.376535511667348e-06
DMWXFREQ_0016 0.004016064256937707
DMWXSIN_0016 -2.3743630552506286e-05 1 4.271751619693612e-06
DMWXCOS_0016 5.740611836492479e-06 1 4.420878217187466e-06
DMWXFREQ_0017 0.004267068272996314
DMWXSIN_0017 1.8768678371242444e-05 1 4.3636488281161415e-06
DMWXCOS_0017 9.88921399388446e-06 1 4.310622420525412e-06
DMWXFREQ_0018 0.00451807228905492
DMWXSIN_0018 -2.55486075219855e-05 1 4.3507913604932865e-06
DMWXCOS_0018 -1.3810396472367071e-05 1 4.337237847049603e-06
DMWXFREQ_0019 0.004769076305113527
DMWXSIN_0019 3.0296942349694885e-05 1 4.303529148160076e-06
DMWXCOS_0019 3.0690467766398084e-05 1 4.387117016111635e-06
DMWXFREQ_0020 0.005020080321172134
DMWXSIN_0020 -9.820394107027999e-06 1 4.431488025878916e-06
DMWXCOS_0020 -8.89274256563913e-06 1 4.246882907071746e-06
DMWXFREQ_0021 0.00527108433723074
DMWXSIN_0021 8.479219229877216e-06 1 4.309511064668012e-06
DMWXCOS_0021 -1.054290113860323e-05 1 4.364632832582027e-06
DMWXFREQ_0022 0.005522088353289347
DMWXSIN_0022 -3.0674486197396964e-06 1 4.411953973307959e-06
DMWXCOS_0022 -7.867863208667077e-06 1 4.260250378102626e-06
DMWXFREQ_0023 0.005773092369347954
DMWXSIN_0023 -3.311917661959573e-06 1 4.350012156150728e-06
DMWXCOS_0023 1.7781018011992893e-05 1 4.3073477481268595e-06
DMWXFREQ_0024 0.0060240963854065604
DMWXSIN_0024 1.2514676274487918e-05 1 4.31406687410843e-06
DMWXCOS_0024 -4.104511221771724e-06 1 4.344027307997093e-06
DMWXFREQ_0025 0.006275100401465167
DMWXSIN_0025 1.6407150114088566e-05 1 4.3896390425630705e-06
DMWXCOS_0025 -8.906402002771867e-06 1 4.271051049810413e-06
DMWXFREQ_0026 0.006526104417523774
DMWXSIN_0026 -6.257209718058854e-06 1 4.324497862947431e-06
DMWXCOS_0026 1.0439759484342997e-06 1 4.338555403884152e-06
DMWXFREQ_0027 0.0067771084335823805
DMWXSIN_0027 -5.185087240484576e-06 1 4.322086318652115e-06
DMWXCOS_0027 -1.0559670675920138e-05 1 4.3356709966292515e-06
DMWXFREQ_0028 0.007028112449640987
DMWXSIN_0028 -1.0467813075940182e-05 1 4.369971947457168e-06
DMWXCOS_0028 5.651156421536296e-06 1 4.2743381785520386e-06
DMWXFREQ_0029 0.007279116465699594
DMWXSIN_0029 3.117054939462295e-06 1 4.338869225619787e-06
DMWXCOS_0029 1.5885114797136547e-05 1 4.3252023311356445e-06
DMWXFREQ_0030 0.007530120481758201
DMWXSIN_0030 6.6410749764222295e-06 1 4.3485573174131024e-06
DMWXCOS_0030 -1.0361995729533963e-05 1 4.304293288075911e-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-06-11T13:47:50.349706
# PINT_version: 0+untagged.358.g3534316
# User: docs
# Host: build-33094749-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-11T13:46:52.386764
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM4
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566857755
FINISH 56985.0000000463691088
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1915.3788043255029
CHI2R 0.990884016722971
TRES 0.986608999272166773
RAJ 5:00:00.00000268 1 0.00000190066703634089
DECJ 15:00:00.00005030 1 0.00016116959119829902
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999998855 1 3.6214478305399477814e-14
F1 -1.000000295184137167e-15 1 8.470114338831499001e-22
PEPOCH 55000.0000000000000000
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.0002580513947990552 1 5.644264446900314e-06
PLANET_SHAPIRO N
DM 14.9999955975508036075 1 5.011542322785246341e-06
TNDMAMP -13.031238700524591 0 0.0423643999694402
TNDMGAM 3.2132482418826775 0 0.2615049186032542
TNDMC 30
[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 0x71353863c090>
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-11T13:48:00.278210
# PINT_version: 0+untagged.358.g3534316
# User: docs
# Host: build-33094749-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-11T13:47:51.099642
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM5
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567547223
FINISH 56985.0000000477440856
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 2019.783302701264
CHI2R 1.0448956558206228
TRES 1.0030578908965874359
RAJ 4:59:59.99999901 1 0.00000143122630949969
DECJ 15:00:00.00004099 1 0.00012359840628712365
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.999999999999959296 1 2.8292859287771445585e-14
F1 -9.999999518312910092e-16 1 6.3689966058646552728e-22
PEPOCH 55000.0000000000000000
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.00015324138027010466 1 4.256674982962686e-06
PLANET_SHAPIRO N
DM 15.0
CM 1.1750715217272361218 1 0.050362716065102501595
TNCHROMIDX 3.5
CMWXEPOCH 55000.0000000000000000
CMWXFREQ_0001 0.0002510040160585242
CMWXSIN_0001 -54.97586832640586 1 0.06614722212403307
CMWXCOS_0001 -8.354705173509526 1 0.06946030352961517
CMWXFREQ_0002 0.0005020080321170484
CMWXSIN_0002 40.43861742419806 1 0.060012162076794835
CMWXCOS_0002 52.469240432200955 1 0.05625676214680594
CMWXFREQ_0003 0.0007530120481755725
CMWXSIN_0003 -6.192501095049541 1 0.05776322863744177
CMWXCOS_0003 -4.459700100650175 1 0.05668986363734605
CMWXFREQ_0004 0.0010040160642340967
CMWXSIN_0004 1.046557520423805 1 0.05851986081060668
CMWXCOS_0004 6.820923129682449 1 0.055110710120421645
CMWXFREQ_0005 0.001255020080292621
CMWXSIN_0005 -1.8925281677493375 1 0.05703650426776246
CMWXCOS_0005 -2.8426362891298926 1 0.0566545270841416
CMWXFREQ_0006 0.001506024096351145
CMWXSIN_0006 2.5110852941143564 1 0.05561966287329346
CMWXCOS_0006 -5.1700342132566295 1 0.057613595381067695
CMWXFREQ_0007 0.0017570281124096694
CMWXSIN_0007 6.386749890508412 1 0.05623021617735105
CMWXCOS_0007 -2.4423325760638637 1 0.05710356242710705
CMWXFREQ_0008 0.0020080321284681934
CMWXSIN_0008 -2.941896875943095 1 0.056031778342793215
CMWXCOS_0008 3.8191475492501428 1 0.05729398686465582
CMWXFREQ_0009 0.0022590361445267177
CMWXSIN_0009 -0.48494653314455505 1 0.05761224825298585
CMWXCOS_0009 1.7540213784457621 1 0.05564048087424295
CMWXFREQ_0010 0.002510040160585242
CMWXSIN_0010 7.570237386435061 1 0.05621805061999051
CMWXCOS_0010 1.823692575031655 1 0.05747623135796964
CMWXFREQ_0011 0.002761044176643766
CMWXSIN_0011 -3.7870516062086423 1 0.07125150063519776
CMWXCOS_0011 -3.5724271887959764 1 0.0704582899805713
CMWXFREQ_0012 0.00301204819270229
CMWXSIN_0012 -0.9539402189481808 1 0.05594515274826638
CMWXCOS_0012 1.5537214670316353 1 0.057554814203294766
CMWXFREQ_0013 0.0032630522087608144
CMWXSIN_0013 3.2764582985126744 1 0.05703759804736903
CMWXCOS_0013 -0.7925310426000259 1 0.056494971377661324
CMWXFREQ_0014 0.0035140562248193387
CMWXSIN_0014 -0.6411952777336041 1 0.05645655816983062
CMWXCOS_0014 0.9403776505545767 1 0.05696133410793995
CMWXFREQ_0015 0.003765060240877863
CMWXSIN_0015 -0.3219137973490672 1 0.056720286753132274
CMWXCOS_0015 0.762561106891925 1 0.05661698638875554
CMWXFREQ_0016 0.004016064256936387
CMWXSIN_0016 0.04077810335367234 1 0.0574380504710672
CMWXCOS_0016 -1.5356586009581568 1 0.05585544185199848
CMWXFREQ_0017 0.004267068272994911
CMWXSIN_0017 -0.37134455558918833 1 0.05711274002014052
CMWXCOS_0017 0.8214212664901883 1 0.056214605808491885
CMWXFREQ_0018 0.004518072289053435
CMWXSIN_0018 0.8371959926224066 1 0.05638364607283639
CMWXCOS_0018 -1.3601984896072723 1 0.05694445270847065
CMWXFREQ_0019 0.00476907630511196
CMWXSIN_0019 0.37768009070130454 1 0.0570767956950141
CMWXCOS_0019 -2.114002375328593 1 0.05644404967551848
CMWXFREQ_0020 0.005020080321170484
CMWXSIN_0020 -0.3273211765765581 1 0.05613771220510653
CMWXCOS_0020 0.07605828321904541 1 0.05708174393420795
CMWXFREQ_0021 0.005271084337229008
CMWXSIN_0021 -0.22580739597521363 1 0.05636638088529479
CMWXCOS_0021 -0.020570125992400757 1 0.056826207041904934
CMWXFREQ_0022 0.005522088353287532
CMWXSIN_0022 -0.32702609007541655 1 0.05648730357798839
CMWXCOS_0022 -1.0634297457074187 1 0.05674281721840948
CMWXFREQ_0023 0.005773092369346056
CMWXSIN_0023 -0.05448830937205539 1 0.05625530016016467
CMWXCOS_0023 0.5798284200838386 1 0.05669596962366156
CMWXFREQ_0024 0.00602409638540458
CMWXSIN_0024 0.2259842543782285 1 0.05616627669291044
CMWXCOS_0024 -0.061370284413552044 1 0.05679191419508297
CMWXFREQ_0025 0.0062751004014631045
CMWXSIN_0025 0.5938369123325025 1 0.05609652298894384
CMWXCOS_0025 0.7214699367635364 1 0.05693143683985742
CMWXFREQ_0026 0.006526104417521629
CMWXSIN_0026 -0.6601417910028969 1 0.05764090927178216
CMWXCOS_0026 -0.6926187997927375 1 0.055438691107442914
CMWXFREQ_0027 0.006777108433580153
CMWXSIN_0027 -0.255085124242085 1 0.056918918834580084
CMWXCOS_0027 0.23141402875732894 1 0.05610374989162022
CMWXFREQ_0028 0.007028112449638677
CMWXSIN_0028 0.16832363384083845 1 0.05725917653558905
CMWXCOS_0028 0.16139106154205402 1 0.055922830535948716
CMWXFREQ_0029 0.007279116465697202
CMWXSIN_0029 -0.08805380086332278 1 0.057322728021979755
CMWXCOS_0029 -0.7234335249734771 1 0.056099941561064605
CMWXFREQ_0030 0.007530120481755726
CMWXSIN_0030 0.4780697260885422 1 0.0575328636837912
CMWXCOS_0030 -0.23413068921042404 1 0.05562178593271261
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-11T13:48:00.339408
# PINT_version: 0+untagged.358.g3534316
# User: docs
# Host: build-33094749-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-11T13:47:51.099642
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM5
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567547223
FINISH 56985.0000000477440856
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 2019.783302701264
CHI2R 1.0448956558206228
TRES 1.0030578908965874359
RAJ 4:59:59.99999901 1 0.00000143122630949969
DECJ 15:00:00.00004099 1 0.00012359840628712365
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.999999999999959296 1 2.8292859287771445585e-14
F1 -9.999999518312910092e-16 1 6.3689966058646552728e-22
PEPOCH 55000.0000000000000000
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.00015324138027010466 1 4.256674982962686e-06
PLANET_SHAPIRO N
DM 15.0
CM 1.1750715217272361218 1 0.050362716065102501595
TNCHROMIDX 3.5
TNCHROMAMP -13.010147523343033 0 0.04013578987524881
TNCHROMGAM 3.151978158691859 0 0.2561443328247285
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 0x71353863c210>
[ ]: