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-05-06T20:20:09.562489
# PINT_version: 0+untagged.322.g9210334
# User: docs
# Host: build-32572216-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-05-06T20:19:33.384102
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM3
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567508797
FINISH 56985.0000000464876042
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1986.895240926832
CHI2R 1.014239530845754
TRES 1.0095838234754384235
RAJ 5:00:00.00025578 1 0.00011549587010344783
DECJ 14:59:59.97346803 1 0.01183068961438333601
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.00000000000100272 1 5.5681539872505851913e-13
F1 -9.995580985906916916e-16 1 1.858747550040015298e-19
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.000000741057712476 1 4.94834771908840113e-06
WXEPOCH 55000.0000000000000000
WXFREQ_0001 0.0002510040160586035
WXSIN_0001 7.188182392562512e-07 1 6.237462163809919e-07
WXCOS_0001 -3.264872994252502e-05 1 1.122866950227744e-05
WXFREQ_0002 0.000502008032117207
WXSIN_0002 -2.608350578235678e-06 1 3.158222079444693e-07
WXCOS_0002 5.596258575502815e-06 1 2.8530135068135727e-06
WXFREQ_0003 0.0007530120481758105
WXSIN_0003 -7.76343809897526e-07 1 2.1958257146374852e-07
WXCOS_0003 -2.2701873140178604e-06 1 1.3073334724037228e-06
WXFREQ_0004 0.001004016064234414
WXSIN_0004 -3.9162421097126566e-07 1 1.7474273513733433e-07
WXCOS_0004 2.255188498346629e-06 1 7.694599642413411e-07
WXFREQ_0005 0.0012550200802930174
WXSIN_0005 6.437230791169769e-07 1 1.5153448647850143e-07
WXCOS_0005 -9.638789105211389e-07 1 5.254331499529449e-07
WXFREQ_0006 0.001506024096351621
WXSIN_0006 -2.623402662337335e-07 1 1.4204135822872187e-07
WXCOS_0006 1.3563466843461294e-06 1 4.004850893153201e-07
WXFREQ_0007 0.0017570281124102245
WXSIN_0007 1.928570896587876e-07 1 1.418470688783745e-07
WXCOS_0007 -5.444760777075522e-07 1 3.34555939609001e-07
WXFREQ_0008 0.002008032128468828
WXSIN_0008 -4.702606051365586e-08 1 1.530703629867642e-07
WXCOS_0008 6.313086936526267e-07 1 3.108636311191118e-07
WXFREQ_0009 0.0022590361445274315
WXSIN_0009 1.9234283126400893e-07 1 1.9099031093582474e-07
WXCOS_0009 -6.581895459477355e-07 1 3.349942351120584e-07
WXFREQ_0010 0.0025100401605860348
WXSIN_0010 -5.968883553333772e-07 1 3.312398434030435e-07
WXCOS_0010 9.681007355606008e-07 1 5.087403036900302e-07
WXFREQ_0011 0.0027610441766446384
WXSIN_0011 -3.765925892453741e-06 1 2.70926286193917e-06
WXCOS_0011 7.618769102725033e-06 1 3.655579202664158e-06
WXFREQ_0012 0.003012048192703242
WXSIN_0012 2.214344356692541e-07 1 1.9618068501400423e-07
WXCOS_0012 -3.7637527619979226e-07 1 2.32963938856772e-07
WXFREQ_0013 0.0032630522087618453
WXSIN_0013 -1.538071436459271e-07 1 9.184776101451053e-08
WXCOS_0013 2.1432176543970581e-07 1 9.669666610311194e-08
WXFREQ_0014 0.003514056224820449
WXSIN_0014 7.13305342228305e-08 1 5.895314296387463e-08
WXCOS_0014 1.942796822397363e-08 1 5.921205651617226e-08
WXFREQ_0015 0.0037650602408790526
WXSIN_0015 -6.394766579917266e-08 1 4.522503885402573e-08
WXCOS_0015 1.2448692541512665e-07 1 4.3576889843248306e-08
WXFREQ_0016 0.004016064256937656
WXSIN_0016 -6.55292384738508e-08 1 3.9707618660326575e-08
WXCOS_0016 1.054007477893173e-08 1 3.833929522068669e-08
WXFREQ_0017 0.004267068272996259
WXSIN_0017 1.220668625830455e-07 1 3.6256305565684997e-08
WXCOS_0017 -2.1491109400306244e-08 1 3.5385465120704896e-08
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.0004796225559498829 1 0.0009285617977242727
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-05-06T20:20:09.592549
# PINT_version: 0+untagged.322.g9210334
# User: docs
# Host: build-32572216-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-05-06T20:19:33.384102
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM3
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567508797
FINISH 56985.0000000464876042
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1986.895240926832
CHI2R 1.014239530845754
TRES 1.0095838234754384235
RAJ 5:00:00.00025578 1 0.00011549587010344783
DECJ 14:59:59.97346803 1 0.01183068961438333601
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.00000000000100272 1 5.5681539872505851913e-13
F1 -9.995580985906916916e-16 1 1.858747550040015298e-19
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.000000741057712476 1 4.94834771908840113e-06
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.0004796225559498829 1 0.0009285617977242727
TNREDAMP -12.746594381311958 0 0.10651000743117006
TNREDGAM 3.8550644840520913 0 0.4886715128769787
TNREDC 17
[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 0x7ff048153d90>
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 = 27
[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-05-06T20:20:51.871076
# PINT_version: 0+untagged.322.g9210334
# User: docs
# Host: build-32572216-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-05-06T20:20:10.002750
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM4
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567227316
FINISH 56985.0000000462840625
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1982.660469777562
CHI2R 1.0225170034953903
TRES 0.9930337181814973046
RAJ 4:59:59.99999863 1 0.00000186390448041385
DECJ 15:00:00.00003938 1 0.00015943900067729035
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.000000000000018124 1 3.5957472339125036196e-14
F1 -1.0000001721746152676e-15 1 8.1632490509007485924e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.000001816907212079 1 4.931235386074550482e-06
DMWXEPOCH 55000.0000000000000000
DMWXFREQ_0001 0.0002510040160586145
DMWXSIN_0001 0.0008946003096600588 1 5.938439672410723e-06
DMWXCOS_0001 -0.0008935861726063368 1 6.767389669648368e-06
DMWXFREQ_0002 0.000502008032117229
DMWXSIN_0002 0.0014288817694810376 1 4.716347466562437e-06
DMWXCOS_0002 -0.000857302885176323 1 4.488770832137023e-06
DMWXFREQ_0003 0.0007530120481758435
DMWXSIN_0003 0.00020303069131812177 1 4.487810267292325e-06
DMWXCOS_0003 0.00033427511908857877 1 4.326514076989947e-06
DMWXFREQ_0004 0.001004016064234458
DMWXSIN_0004 -0.0002239804775318629 1 4.430686035506522e-06
DMWXCOS_0004 -0.00025485793661216905 1 4.285112426354928e-06
DMWXFREQ_0005 0.0012550200802930725
DMWXSIN_0005 -0.000226418962488874 1 4.319866507504622e-06
DMWXCOS_0005 0.00012838950617389844 1 4.3519589599466964e-06
DMWXFREQ_0006 0.001506024096351687
DMWXSIN_0006 -0.0001343999614218248 1 4.417441071878044e-06
DMWXCOS_0006 -0.00015675534352189326 1 4.220500683310905e-06
DMWXFREQ_0007 0.0017570281124103014
DMWXSIN_0007 0.00011599339803708177 1 4.380485248422713e-06
DMWXCOS_0007 2.10326312186497e-05 1 4.240339945975684e-06
DMWXFREQ_0008 0.002008032128468916
DMWXSIN_0008 7.1503820093445e-05 1 4.302101270011841e-06
DMWXCOS_0008 5.027307747523312e-06 1 4.313564679317319e-06
DMWXFREQ_0009 0.0022590361445275304
DMWXSIN_0009 -6.777167312946779e-05 1 4.325733283695729e-06
DMWXCOS_0009 -9.492001226651565e-07 1 4.299746752535506e-06
DMWXFREQ_0010 0.002510040160586145
DMWXSIN_0010 3.977592687172212e-05 1 4.350394677397452e-06
DMWXCOS_0010 -0.00010402397216093929 1 4.3495585858025856e-06
DMWXFREQ_0011 0.0027610441766447594
DMWXSIN_0011 5.066250158153932e-05 1 6.810428957473435e-06
DMWXCOS_0011 -5.5943966715807864e-05 1 6.996016822004043e-06
DMWXFREQ_0012 0.003012048192703374
DMWXSIN_0012 -1.564824317206421e-05 1 4.316135751549212e-06
DMWXCOS_0012 8.529142018964935e-06 1 4.32514868933664e-06
DMWXFREQ_0013 0.0032630522087619884
DMWXSIN_0013 2.9215309794474335e-05 1 4.338725581864244e-06
DMWXCOS_0013 -2.3081287621522036e-05 1 4.263271872857371e-06
DMWXFREQ_0014 0.003514056224820603
DMWXSIN_0014 -1.4886118254284707e-05 1 4.278808950716149e-06
DMWXCOS_0014 2.6510992725988513e-05 1 4.316352439330853e-06
DMWXFREQ_0015 0.0037650602408792174
DMWXSIN_0015 -4.0551078896910996e-05 1 4.384244853962622e-06
DMWXCOS_0015 -1.3952467734372883e-05 1 4.199409537705386e-06
DMWXFREQ_0016 0.004016064256937832
DMWXSIN_0016 -1.6849035883705655e-05 1 4.2231954079889334e-06
DMWXCOS_0016 1.1034030348324043e-05 1 4.359365730645938e-06
DMWXFREQ_0017 0.004267068272996446
DMWXSIN_0017 1.4402663805571534e-05 1 4.2609502359278775e-06
DMWXCOS_0017 -6.969530589994863e-06 1 4.324451874846207e-06
DMWXFREQ_0018 0.004518072289055061
DMWXSIN_0018 1.4962594776748035e-05 1 4.353155493049104e-06
DMWXCOS_0018 2.8811538196651236e-06 1 4.218598238538541e-06
DMWXFREQ_0019 0.004769076305113675
DMWXSIN_0019 1.1013299077129405e-05 1 4.271928630367011e-06
DMWXCOS_0019 -1.2347716541960371e-05 1 4.30669872114703e-06
DMWXFREQ_0020 0.00502008032117229
DMWXSIN_0020 1.3072702174647734e-05 1 4.293016553286559e-06
DMWXCOS_0020 -1.0434403033821246e-05 1 4.292819768361764e-06
DMWXFREQ_0021 0.005271084337230904
DMWXSIN_0021 -1.501965394171976e-07 1 4.217234337733182e-06
DMWXCOS_0021 -2.2872812351160824e-05 1 4.367723313225131e-06
DMWXFREQ_0022 0.005522088353289519
DMWXSIN_0022 2.2886531851164675e-07 1 4.3111114823601745e-06
DMWXCOS_0022 1.2802331107290144e-06 1 4.27633732485273e-06
DMWXFREQ_0023 0.005773092369348133
DMWXSIN_0023 9.278396020182633e-06 1 4.3424093964478025e-06
DMWXCOS_0023 1.8717659286732253e-05 1 4.239500644624376e-06
DMWXFREQ_0024 0.006024096385406748
DMWXSIN_0024 -4.898567901648147e-06 1 4.278595048470354e-06
DMWXCOS_0024 7.062595726284702e-06 1 4.315007960119053e-06
DMWXFREQ_0025 0.006275100401465362
DMWXSIN_0025 -7.639796640451749e-06 1 4.2440555470425926e-06
DMWXCOS_0025 -1.4050330610079399e-05 1 4.343622067869429e-06
DMWXFREQ_0026 0.006526104417523977
DMWXSIN_0026 -2.0900289881828715e-06 1 4.315363255757676e-06
DMWXCOS_0026 -4.024531389434772e-06 1 4.277029841889079e-06
DMWXFREQ_0027 0.006777108433582591
DMWXSIN_0027 2.3587741594770703e-06 1 4.293682968200411e-06
DMWXCOS_0027 -1.1437810506120345e-05 1 4.300620863118548e-06
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.00040212877089185667 1 5.462679129400423e-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)
27
[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-05-06T20:20:51.909093
# PINT_version: 0+untagged.322.g9210334
# User: docs
# Host: build-32572216-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-05-06T20:20:10.002750
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM4
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567227316
FINISH 56985.0000000462840625
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1982.660469777562
CHI2R 1.0225170034953903
TRES 0.9930337181814973046
RAJ 4:59:59.99999863 1 0.00000186390448041385
DECJ 15:00:00.00003938 1 0.00015943900067729035
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.000000000000018124 1 3.5957472339125036196e-14
F1 -1.0000001721746152676e-15 1 8.1632490509007485924e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.000001816907212079 1 4.931235386074550482e-06
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.00040212877089185667 1 5.462679129400423e-06
TNDMAMP -13.035506633715427 0 0.0447740588825518
TNDMGAM 3.6941893394811083 0 0.2832598307932774
TNDMC 27
[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()
27 27
[18]:
<matplotlib.legend.Legend at 0x7ff047e3cbd0>
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-05-06T20:20:59.193325
# PINT_version: 0+untagged.322.g9210334
# User: docs
# Host: build-32572216-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-05-06T20:20:52.207581
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM5
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566827431
FINISH 56985.0000000457464583
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 2084.143465080816
CHI2R 1.078191135582419
TRES 1.0113875564766812932
RAJ 5:00:00.00000077 1 0.00000142072284024247
DECJ 14:59:59.99992485 1 0.00012459886351702087
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999997136 1 2.7839371607385480573e-14
F1 -1.0000009347858719237e-15 1 6.3517032863434461223e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.0
CM 1.15678794537224353 1 0.04986326327324507035
TNCHROMIDX 3.5
CMWXEPOCH 55000.0000000000000000
CMWXFREQ_0001 0.0002510040160586461
CMWXSIN_0001 1.3680173117105283 1 0.06502532440162262
CMWXCOS_0001 141.28098776875626 1 0.06963251609835268
CMWXFREQ_0002 0.0005020080321172922
CMWXSIN_0002 -35.87810802003495 1 0.05975806239636146
CMWXCOS_0002 72.85215736403693 1 0.05638723600694139
CMWXFREQ_0003 0.0007530120481759383
CMWXSIN_0003 -17.388044061702573 1 0.0594573987724311
CMWXCOS_0003 -9.200038784855405 1 0.05500131352906953
CMWXFREQ_0004 0.0010040160642345844
CMWXSIN_0004 18.0568866468496 1 0.05699483172554018
CMWXCOS_0004 7.786214766374167 1 0.057109622366869714
CMWXFREQ_0005 0.0012550200802932305
CMWXSIN_0005 14.636766246886609 1 0.056944569957652916
CMWXCOS_0005 -3.8555429569081197 1 0.05701802281109918
CMWXFREQ_0006 0.0015060240963518767
CMWXSIN_0006 12.57280307949233 1 0.058111349337888925
CMWXCOS_0006 5.194699657197391 1 0.05561739537011687
CMWXFREQ_0007 0.0017570281124105228
CMWXSIN_0007 -0.3964966917969004 1 0.05685584797980621
CMWXCOS_0007 0.6885062156776807 1 0.05671617714548486
CMWXFREQ_0008 0.0020080321284691688
CMWXSIN_0008 -4.134805013922262 1 0.05649813788365906
CMWXCOS_0008 -2.5343207718934484 1 0.056973965874742655
CMWXFREQ_0009 0.002259036144527815
CMWXSIN_0009 -0.6164560957445328 1 0.0578516525654616
CMWXCOS_0009 1.18501117096347 1 0.05579191512510528
CMWXFREQ_0010 0.002510040160586461
CMWXSIN_0010 0.18553444174326164 1 0.058030273249647504
CMWXCOS_0010 0.8099923985553042 1 0.055923383017513845
CMWXFREQ_0011 0.0027610441766451072
CMWXSIN_0011 -0.590713546873559 1 0.06947025080157661
CMWXCOS_0011 -2.1308250368979627 1 0.07157077477096273
CMWXFREQ_0012 0.0030120481927037534
CMWXSIN_0012 1.9349184836472526 1 0.057147548718208735
CMWXCOS_0012 0.6273631365489568 1 0.05652720443083328
CMWXFREQ_0013 0.0032630522087623995
CMWXSIN_0013 1.3446859115038918 1 0.05805753720958289
CMWXCOS_0013 -0.8422439674648834 1 0.05536579865134628
CMWXFREQ_0014 0.0035140562248210457
CMWXSIN_0014 -1.2789038790898013 1 0.05647479353911161
CMWXCOS_0014 -1.1748362854538554 1 0.056926989238805126
CMWXFREQ_0015 0.003765060240879692
CMWXSIN_0015 -0.03108189695994524 1 0.05620674355466351
CMWXCOS_0015 -0.7742887607095864 1 0.05715934757093877
CMWXFREQ_0016 0.0040160642569383375
CMWXSIN_0016 0.18136242409979286 1 0.056531272431861915
CMWXCOS_0016 -0.4819038964668708 1 0.05664592255140411
CMWXFREQ_0017 0.004267068272996984
CMWXSIN_0017 0.5422205960762 1 0.05809787094820577
CMWXCOS_0017 0.5898970478188008 1 0.0551200596057319
CMWXFREQ_0018 0.00451807228905563
CMWXSIN_0018 -0.22660367270167897 1 0.056560344015593914
CMWXCOS_0018 0.6582020584335326 1 0.05677500656744321
CMWXFREQ_0019 0.004769076305114276
CMWXSIN_0019 -0.6620242478438013 1 0.05617434812130137
CMWXCOS_0019 -1.1758931135613655 1 0.0569801520845545
CMWXFREQ_0020 0.005020080321172922
CMWXSIN_0020 0.36158895896547966 1 0.057076865915697474
CMWXCOS_0020 -0.5247084874531383 1 0.05609073521198633
CMWXFREQ_0021 0.005271084337231569
CMWXSIN_0021 -1.1596874475809786 1 0.056536127407327365
CMWXCOS_0021 0.2777819452889487 1 0.05661618109884885
CMWXFREQ_0022 0.0055220883532902144
CMWXSIN_0022 -0.9242395591405196 1 0.0566829584643604
CMWXCOS_0022 -0.47463611946893913 1 0.056557905034141036
CMWXFREQ_0023 0.005773092369348861
CMWXSIN_0023 -0.7349465517526166 1 0.0562875405952251
CMWXCOS_0023 -0.6110551269559346 1 0.05703194661686489
CMWXFREQ_0024 0.006024096385407507
CMWXSIN_0024 -0.4693873877486135 1 0.05733351960074186
CMWXCOS_0024 -0.07894104139328442 1 0.055669195351967495
CMWXFREQ_0025 0.006275100401466153
CMWXSIN_0025 0.5557457909974367 1 0.057764239932756374
CMWXCOS_0025 0.4315793654706563 1 0.05505493360129846
CMWXFREQ_0026 0.006526104417524799
CMWXSIN_0026 -0.3360892393278689 1 0.05546706894670295
CMWXCOS_0026 -0.49040959186776406 1 0.057364155103677726
CMWXFREQ_0027 0.006777108433583446
CMWXSIN_0027 -0.221345435236613 1 0.05673502382211537
CMWXCOS_0027 0.0395296318822882 1 0.05623999413034816
CMWXFREQ_0028 0.007028112449642091
CMWXSIN_0028 -0.43615259073303414 1 0.05670804755165782
CMWXCOS_0028 0.4047451106916018 1 0.055921636776296274
CMWXFREQ_0029 0.007279116465700738
CMWXSIN_0029 0.26740773946295987 1 0.05585086052496827
CMWXCOS_0029 -0.8110201231150055 1 0.0568753394972691
CMWXFREQ_0030 0.007530120481759384
CMWXSIN_0030 -0.200546846918701 1 0.05711692364293363
CMWXCOS_0030 0.2509676525715172 1 0.05544784054623583
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.0008356015602644258 1 4.253837354397047e-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-05-06T20:20:59.232875
# PINT_version: 0+untagged.322.g9210334
# User: docs
# Host: build-32572216-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-05-06T20:20:52.207581
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM5
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566827431
FINISH 56985.0000000457464583
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 2084.143465080816
CHI2R 1.078191135582419
TRES 1.0113875564766812932
RAJ 5:00:00.00000077 1 0.00000142072284024247
DECJ 14:59:59.99992485 1 0.00012459886351702087
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999997136 1 2.7839371607385480573e-14
F1 -1.0000009347858719237e-15 1 6.3517032863434461223e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.0
CM 1.15678794537224353 1 0.04986326327324507035
TNCHROMIDX 3.5
TNCHROMAMP -13.047630207103795 0 0.04018296312896732
TNCHROMGAM 3.579468312004918 0 0.20879414053522027
TNCHROMC 30
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.0008356015602644258 1 4.253837354397047e-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 0x7ff0540e3850>
[ ]: