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 = 15
[5]:
print(np.argmin(d_aics))
15
[6]:
# The Y axis is plotted in log scale only for better visibility.
plt.scatter(list(range(len(d_aics))), d_aics + 1)
plt.axvline(nharm_opt, color="red", label="Optimum number of harmonics")
plt.axvline(
int(m.TNREDC.value), color="black", ls="--", label="Injected number of harmonics"
)
plt.xlabel("Number of harmonics")
plt.ylabel("AIC - AIC$_\\min{} + 1$")
plt.legend()
plt.yscale("log")
# plt.savefig("sim3-aic.pdf")
[7]:
# Now create a new model with the optimum number of harmonics
m2 = deepcopy(m1)
Tspan = t.get_mjds().max() - t.get_mjds().min()
wavex_setup(m2, T_span=Tspan, n_freqs=nharm_opt, freeze_params=False)
ftr = WLSFitter(t, m2)
ftr.fit_toas(maxiter=10)
m2 = ftr.model
print(m2)
# Created: 2026-08-21T20:06:04.328203
# PINT_version: 1.1.6+68.g7d55691
# User: docs
# Host: build-34177976-project-85767-nanograv-pint
# OS: Linux-7.0.0-1004-aws-x86_64-with-glibc2.35
# Python: 3.11.15 (main, Jun 25 2026, 19:09:59) [GCC 11.4.0]
# Format: pint
# read_time: 2026-08-21T20:05:16.005458
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
# ell1h_shapiro: full
PSR SIM3
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.999999956712384286
FINISH 56985.000000046484340278
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1945.4218674718695
CHI2R 0.9910452712541363
TRES 0.9871112688146033
RAJ 4:59:59.99997587 1 0.00009480029673717161
DECJ 15:00:00.01191873 1 0.01088640122637804866
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.00000000000020244 1 4.2648431135885404138e-13
F1 -1.0001099951150653704e-15 1 1.764558762570160568e-19
PEPOCH 55000.000000000000000000
TZRMJD 55000.000000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.00010532999550494171 1 0.0008814004297344163
PLANET_SHAPIRO N
DM 14.999996318137779685 1 4.6410813833846497703e-06
WXEPOCH 55000.000000000000000000
WXFREQ_0001 0.0002510040160586012
WXSIN_0001 1.4981256979903787e-06 1 4.778784320506641e-07
WXCOS_0001 5.086986692052486e-06 1 1.0659140982960933e-05
WXFREQ_0002 0.0005020080321172024
WXSIN_0002 7.000786453308075e-07 1 2.428999180397423e-07
WXCOS_0002 -9.808476974406353e-07 1 2.70807770251443e-06
WXFREQ_0003 0.0007530120481758036
WXSIN_0003 -7.74327060748152e-07 1 1.687403825597217e-07
WXCOS_0003 2.463905791632959e-07 1 1.2414251607573212e-06
WXFREQ_0004 0.0010040160642344048
WXSIN_0004 -5.251724927430215e-07 1 1.3443500019355656e-07
WXCOS_0004 -3.866714285916211e-07 1 7.296629642673875e-07
WXFREQ_0005 0.0012550200802930059
WXSIN_0005 -4.016256589349985e-07 1 1.1783837962202977e-07
WXCOS_0005 2.67621431817593e-07 1 4.983441101599158e-07
WXFREQ_0006 0.0015060240963516072
WXSIN_0006 -6.589783429647492e-08 1 1.1000917212588875e-07
WXCOS_0006 -3.26653809071732e-07 1 3.7825884687031433e-07
WXFREQ_0007 0.0017570281124102084
WXSIN_0007 2.938179664113286e-07 1 1.0876111357731823e-07
WXCOS_0007 3.223648952995812e-07 1 3.1611237526575995e-07
WXFREQ_0008 0.0020080321284688097
WXSIN_0008 3.13693540298604e-08 1 1.1729310294923962e-07
WXCOS_0008 -2.3444646883790417e-07 1 2.937453884141349e-07
WXFREQ_0009 0.0022590361445274107
WXSIN_0009 -1.43137535998256e-07 1 1.4324665272738228e-07
WXCOS_0009 2.992931084203384e-07 1 3.1611908892025495e-07
WXFREQ_0010 0.0025100401605860118
WXSIN_0010 -2.0047593095448473e-07 1 2.454493015995973e-07
WXCOS_0010 -5.336225637495138e-07 1 4.7283919384542393e-07
WXFREQ_0011 0.0027610441766446133
WXSIN_0011 -1.41668188235869e-06 1 1.9947391054615076e-06
WXCOS_0011 -3.951893320058739e-06 1 3.3645431797409706e-06
WXFREQ_0012 0.0030120481927032143
WXSIN_0012 2.1317043516374468e-08 1 1.449067881138583e-07
WXCOS_0012 3.8576718408619216e-07 1 2.1147516398579864e-07
WXFREQ_0013 0.003263052208761816
WXSIN_0013 -2.6437981464973203e-08 1 7.024945574473811e-08
WXCOS_0013 -7.566080144052283e-08 1 8.626666022870887e-08
WXFREQ_0014 0.003514056224820417
WXSIN_0014 -2.090997528826084e-09 1 4.871601772640997e-08
WXCOS_0014 1.9342838905247063e-07 1 5.173618939881565e-08
WXFREQ_0015 0.003765060240879018
WXSIN_0015 -2.829339068706161e-08 1 3.9770101680644006e-08
WXCOS_0015 -9.054089271169725e-08 1 3.9643699762231166e-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)
15
[9]:
# We can create a `PLRedNoise` model from the `WaveX` model.
# This will estimate the spectral parameters from the `WaveX`
# amplitudes.
m3 = plrednoise_from_wavex(m2)
print(m3)
# Created: 2026-08-21T20:06:04.373937
# PINT_version: 1.1.6+68.g7d55691
# User: docs
# Host: build-34177976-project-85767-nanograv-pint
# OS: Linux-7.0.0-1004-aws-x86_64-with-glibc2.35
# Python: 3.11.15 (main, Jun 25 2026, 19:09:59) [GCC 11.4.0]
# Format: pint
# read_time: 2026-08-21T20:05:16.005458
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
# ell1h_shapiro: full
PSR SIM3
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.999999956712384286
FINISH 56985.000000046484340278
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1945.4218674718695
CHI2R 0.9910452712541363
TRES 0.9871112688146033
RAJ 4:59:59.99997587 1 0.00009480029673717161
DECJ 15:00:00.01191873 1 0.01088640122637804866
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.00000000000020244 1 4.2648431135885404138e-13
F1 -1.0001099951150653704e-15 1 1.764558762570160568e-19
PEPOCH 55000.000000000000000000
TZRMJD 55000.000000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.00010532999550494171 1 0.0008814004297344163
PLANET_SHAPIRO N
DM 14.999996318137779685 1 4.6410813833846497703e-06
TNREDAMP -12.912742958980543 0 0.1257112385758392
TNREDGAM 2.396298760449841 0 0.5759199414297901
TNREDC 15
[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()
15 15
[10]:
<matplotlib.legend.Legend at 0x70102c248fd0>
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 = 23
[14]:
# The Y axis is plotted in log scale only for better visibility.
plt.scatter(list(range(len(d_aics))), d_aics + 1)
plt.axvline(nharm_opt, color="red", label="Optimum number of harmonics")
plt.axvline(
int(m.TNDMC.value), color="black", ls="--", label="Injected number of harmonics"
)
plt.xlabel("Number of harmonics")
plt.ylabel("AIC - AIC$_\\min{} + 1$")
plt.legend()
plt.yscale("log")
# plt.savefig("sim3-aic.pdf")
[15]:
# Now create a new model with the optimum number of
# harmonics
m2 = deepcopy(m1)
Tspan = t.get_mjds().max() - t.get_mjds().min()
dmwavex_setup(m2, T_span=Tspan, n_freqs=nharm_opt, freeze_params=False)
ftr = WLSFitter(t, m2)
ftr.fit_toas(maxiter=10)
m2 = ftr.model
print(m2)
# Created: 2026-08-21T20:07:01.341447
# PINT_version: 1.1.6+68.g7d55691
# User: docs
# Host: build-34177976-project-85767-nanograv-pint
# OS: Linux-7.0.0-1004-aws-x86_64-with-glibc2.35
# Python: 3.11.15 (main, Jun 25 2026, 19:09:59) [GCC 11.4.0]
# Format: pint
# read_time: 2026-08-21T20:06:04.997655
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
# ell1h_shapiro: full
PSR SIM4
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.999999956709201387
FINISH 56985.000000046490081019
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 2054.66527870995
CHI2R 1.0552980373446073
TRES 1.0173437863245316
RAJ 4:59:59.99999796 1 0.00000191408616866844
DECJ 15:00:00.00026399 1 0.00016353477420529542
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999996306 1 3.6850956266004358225e-14
F1 -1.0000014494460157477e-15 1 8.560683320637460527e-22
PEPOCH 55000.000000000000000000
TZRMJD 55000.000000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.0002315100251818588 1 5.6200051801076834e-06
PLANET_SHAPIRO N
DM 15.0000025687513062105 1 4.9328022034904202875e-06
DMWXEPOCH 55000.000000000000000000
DMWXFREQ_0001 0.0002510040160586003
DMWXSIN_0001 -0.004490815852208149 1 5.936455752867055e-06
DMWXCOS_0001 -0.0011621965622247278 1 6.9113720034977375e-06
DMWXFREQ_0002 0.0005020080321172006
DMWXSIN_0002 0.0007311178838694184 1 4.658959499449195e-06
DMWXCOS_0002 -0.0005476088170965188 1 4.489294213372442e-06
DMWXFREQ_0003 0.0007530120481758009
DMWXSIN_0003 -0.00022621192213471668 1 4.529308954012159e-06
DMWXCOS_0003 0.0005574752901044574 1 4.22476817439032e-06
DMWXFREQ_0004 0.0010040160642344012
DMWXSIN_0004 0.0001981606644559563 1 4.368641276512823e-06
DMWXCOS_0004 -3.946872878135618e-05 1 4.259467362768826e-06
DMWXFREQ_0005 0.0012550200802930013
DMWXSIN_0005 0.00014326250933671377 1 4.3540798481599945e-06
DMWXCOS_0005 -0.00010032289699526552 1 4.2318733867909635e-06
DMWXFREQ_0006 0.0015060240963516017
DMWXSIN_0006 -0.00013375691823890474 1 4.270676444437832e-06
DMWXCOS_0006 0.0001891503931512572 1 4.252708409046796e-06
DMWXFREQ_0007 0.001757028112410202
DMWXSIN_0007 -3.9883347545526826e-05 1 4.349324464027162e-06
DMWXCOS_0007 3.5757155508060344e-05 1 4.155993814973536e-06
DMWXFREQ_0008 0.0020080321284688023
DMWXSIN_0008 -5.286816614964641e-05 1 4.227019592806262e-06
DMWXCOS_0008 -7.226817255311929e-05 1 4.289006944140098e-06
DMWXFREQ_0009 0.0022590361445274025
DMWXSIN_0009 2.772704431480763e-05 1 4.156438276830094e-06
DMWXCOS_0009 -5.9694668253846356e-06 1 4.3624271953783014e-06
DMWXFREQ_0010 0.0025100401605860027
DMWXSIN_0010 8.453272226965111e-05 1 4.228035738985233e-06
DMWXCOS_0010 7.375068257812361e-05 1 4.343038570913497e-06
DMWXFREQ_0011 0.0027610441766446033
DMWXSIN_0011 -1.1957097222060481e-05 1 6.9322567141798065e-06
DMWXCOS_0011 -3.989893542786022e-05 1 6.886045542262979e-06
DMWXFREQ_0012 0.0030120481927032035
DMWXSIN_0012 5.2767724053406646e-05 1 4.295710480442024e-06
DMWXCOS_0012 2.083517488636323e-05 1 4.214417628696752e-06
DMWXFREQ_0013 0.0032630522087618036
DMWXSIN_0013 1.2648607313548299e-05 1 4.256299975352864e-06
DMWXCOS_0013 -2.4901361533828016e-05 1 4.232576888256289e-06
DMWXFREQ_0014 0.003514056224820404
DMWXSIN_0014 7.993631787168164e-06 1 4.1387965079747556e-06
DMWXCOS_0014 -1.2788911635529563e-07 1 4.32453162249171e-06
DMWXFREQ_0015 0.0037650602408790044
DMWXSIN_0015 2.9994062665383558e-05 1 4.2620489330267836e-06
DMWXCOS_0015 1.7418916410070406e-05 1 4.206132636923229e-06
DMWXFREQ_0016 0.004016064256937605
DMWXSIN_0016 1.1003299691818269e-05 1 4.1378350001383284e-06
DMWXCOS_0016 -1.4580065979678517e-05 1 4.331120733671254e-06
DMWXFREQ_0017 0.004267068272996205
DMWXSIN_0017 -2.2684461216274025e-06 1 4.259383369222184e-06
DMWXCOS_0017 3.2377347426431546e-06 1 4.212824379560192e-06
DMWXFREQ_0018 0.004518072289054805
DMWXSIN_0018 -1.2810351658152087e-05 1 4.229982454888929e-06
DMWXCOS_0018 9.137887047601204e-06 1 4.237729604690205e-06
DMWXFREQ_0019 0.004769076305113406
DMWXSIN_0019 8.780201135706432e-06 1 4.308064347536147e-06
DMWXCOS_0019 -1.2681863423567164e-05 1 4.152694718037753e-06
DMWXFREQ_0020 0.005020080321172005
DMWXSIN_0020 -3.0798061635405427e-06 1 4.231388449078961e-06
DMWXCOS_0020 -4.800283968104186e-06 1 4.2253221077076595e-06
DMWXFREQ_0021 0.005271084337230606
DMWXSIN_0021 1.3977667091820053e-05 1 4.275931123615049e-06
DMWXCOS_0021 1.4559170112230357e-05 1 4.17913353669131e-06
DMWXFREQ_0022 0.005522088353289207
DMWXSIN_0022 -1.585238909024128e-05 1 4.269534422457546e-06
DMWXCOS_0022 8.323744008988871e-06 1 4.187913922629938e-06
DMWXFREQ_0023 0.005773092369347806
DMWXSIN_0023 -7.856910670395158e-06 1 4.206983763815748e-06
DMWXCOS_0023 -9.909292078097627e-06 1 4.239698127030762e-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)
23
[17]:
# We can create a `PLDMNoise` model from the `DMWaveX` model.
# This will estimate the spectral parameters from the `DMWaveX`
# amplitudes.
m3 = pldmnoise_from_dmwavex(m2)
print(m3)
# Created: 2026-08-21T20:07:01.396986
# PINT_version: 1.1.6+68.g7d55691
# User: docs
# Host: build-34177976-project-85767-nanograv-pint
# OS: Linux-7.0.0-1004-aws-x86_64-with-glibc2.35
# Python: 3.11.15 (main, Jun 25 2026, 19:09:59) [GCC 11.4.0]
# Format: pint
# read_time: 2026-08-21T20:06:04.997655
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
# ell1h_shapiro: full
PSR SIM4
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.999999956709201387
FINISH 56985.000000046490081019
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 2054.66527870995
CHI2R 1.0552980373446073
TRES 1.0173437863245316
RAJ 4:59:59.99999796 1 0.00000191408616866844
DECJ 15:00:00.00026399 1 0.00016353477420529542
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999996306 1 3.6850956266004358225e-14
F1 -1.0000014494460157477e-15 1 8.560683320637460527e-22
PEPOCH 55000.000000000000000000
TZRMJD 55000.000000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.0002315100251818588 1 5.6200051801076834e-06
PLANET_SHAPIRO N
DM 15.0000025687513062105 1 4.9328022034904202875e-06
TNDMAMP -13.08402371423738 0 0.04941856518747488
TNDMGAM 3.7886730982954306 0 0.2816837368057958
TNDMC 23
[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()
23 23
[18]:
<matplotlib.legend.Legend at 0x70102b5ed190>
Chromatic noise fitting
Let us now do a similar kind of analysis for chromatic noise.
[19]:
par_sim = """
PSR SIM5
RAJ 05:00:00 1
DECJ 15:00:00 1
PEPOCH 55000
F0 100 1
F1 -1e-15 1
PHOFF 0 1
DM 15
CM 1.2 1
TNCHROMIDX 3.5
TNCHROMAMP -13
TNCHROMGAM 3.5
TNCHROMC 30
TZRMJD 55000
TZRFRQ 1400
TZRSITE gbt
UNITS TDB
EPHEM DE440
CLOCK TT(BIPM2019)
"""
m = get_model(StringIO(par_sim))
[20]:
# Generate the simulated TOAs.
ntoas = 2000
toaerrs = np.random.uniform(0.5, 2.0, ntoas) * u.us
freqs = np.linspace(500, 1500, 8) * u.MHz
t = make_fake_toas_uniform(
startMJD=53001,
endMJD=57001,
ntoas=ntoas,
model=m,
freq=freqs,
obs="gbt",
error=toaerrs,
add_noise=True,
add_correlated_noise=True,
name="fake",
include_bipm=True,
multi_freqs_in_epoch=True,
)
[21]:
# Find the optimum number of harmonics by minimizing AIC.
m1 = deepcopy(m)
m1.remove_component("PLChromNoise")
m2 = deepcopy(m1)
nharm_opt = m.TNCHROMC.value
[22]:
# Now create a new model with the optimum number of
# harmonics
m2 = deepcopy(m1)
Tspan = t.get_mjds().max() - t.get_mjds().min()
cmwavex_setup(m2, T_span=Tspan, n_freqs=nharm_opt, freeze_params=False)
ftr = WLSFitter(t, m2)
ftr.fit_toas(maxiter=10)
m2 = ftr.model
print(m2)
# Created: 2026-08-21T20:07:11.340435
# PINT_version: 1.1.6+68.g7d55691
# User: docs
# Host: build-34177976-project-85767-nanograv-pint
# OS: Linux-7.0.0-1004-aws-x86_64-with-glibc2.35
# Python: 3.11.15 (main, Jun 25 2026, 19:09:59) [GCC 11.4.0]
# Format: pint
# read_time: 2026-08-21T20:07:02.059953
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
# ell1h_shapiro: full
PSR SIM5
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.999999956787824162
FINISH 56985.000000050841412037
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1924.7195246235137
CHI2R 0.9957162569185275
TRES 0.9877706004952963
RAJ 4:59:59.99999955 1 0.00000146238128818959
DECJ 15:00:00.00001146 1 0.00012322699426034855
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.000000000000006134 1 2.8040926342748202197e-14
F1 -1.0000000505979882195e-15 1 6.404970379774501803e-22
PEPOCH 55000.000000000000000000
TZRMJD 55000.000000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.000743291662428148 1 4.349327332289849e-06
PLANET_SHAPIRO N
DM 15.0
CM 1.1040354831558155475 1 0.050400879872884565835
TNCHROMIDX 3.5
CMWXEPOCH 55000.000000000000000000
CMWXFREQ_0001 0.0002510040160583312
CMWXSIN_0001 49.54975457024733 1 0.06591872777923594
CMWXCOS_0001 -226.57657261295776 1 0.0698246392126758
CMWXFREQ_0002 0.0005020080321166624
CMWXSIN_0002 -8.338435910723252 1 0.06059017509992533
CMWXCOS_0002 8.536701745408962 1 0.05653640908215501
CMWXFREQ_0003 0.0007530120481749936
CMWXSIN_0003 23.20721966913107 1 0.05888805535184548
CMWXCOS_0003 13.44928590439415 1 0.05648684933841949
CMWXFREQ_0004 0.0010040160642333248
CMWXSIN_0004 2.563114201412069 1 0.05781193611987641
CMWXCOS_0004 8.355521250004228 1 0.05669825638513653
CMWXFREQ_0005 0.001255020080291656
CMWXSIN_0005 -20.69415093544774 1 0.05718909815532079
CMWXCOS_0005 -8.688708053865552 1 0.05689264636726238
CMWXFREQ_0006 0.0015060240963499871
CMWXSIN_0006 1.4758708182365345 1 0.056164632034032776
CMWXCOS_0006 4.649165777068334 1 0.05754962887906138
CMWXFREQ_0007 0.0017570281124083184
CMWXSIN_0007 7.494319313833982 1 0.056876937452190385
CMWXCOS_0007 -0.4262247083328983 1 0.05682055941791669
CMWXFREQ_0008 0.0020080321284666495
CMWXSIN_0008 1.2583913973531684 1 0.0574099281917559
CMWXCOS_0008 5.335313968299185 1 0.056211495267241555
CMWXFREQ_0009 0.002259036144524981
CMWXSIN_0009 -0.8780216644039774 1 0.05655554108315636
CMWXCOS_0009 1.8622766790709941 1 0.05692496286274216
CMWXFREQ_0010 0.002510040160583312
CMWXSIN_0010 3.5937868107429423 1 0.057457481492545766
CMWXCOS_0010 3.473814008335264 1 0.056235209797714236
CMWXFREQ_0011 0.002761044176641643
CMWXSIN_0011 1.7776761806532637 1 0.07178070099601024
CMWXCOS_0011 0.5036729308491233 1 0.06944273619488539
CMWXFREQ_0012 0.0030120481926999743
CMWXSIN_0012 -1.8992076165545446 1 0.0577838279780376
CMWXCOS_0012 0.515642271759958 1 0.05562947215559465
CMWXFREQ_0013 0.0032630522087583056
CMWXSIN_0013 0.9023111994364673 1 0.05683112391548735
CMWXCOS_0013 0.17815994828800744 1 0.05649166204982716
CMWXFREQ_0014 0.003514056224816637
CMWXSIN_0014 -2.616935790281562 1 0.05619892354204903
CMWXCOS_0014 1.8917133542127316 1 0.05712659691394667
CMWXFREQ_0015 0.003765060240874968
CMWXSIN_0015 -0.4316384022194187 1 0.05683148796131015
CMWXCOS_0015 -0.904483426710036 1 0.056697937873983695
CMWXFREQ_0016 0.004016064256933299
CMWXSIN_0016 1.6593149039064974 1 0.057594812236980514
CMWXCOS_0016 -1.7525272206318072 1 0.055737649978773236
CMWXFREQ_0017 0.004267068272991631
CMWXSIN_0017 0.3554183313641823 1 0.05734374473433686
CMWXCOS_0017 1.216978659126739 1 0.05597132248874121
CMWXFREQ_0018 0.004518072289049962
CMWXSIN_0018 -1.030435342233712 1 0.05670819970638455
CMWXCOS_0018 1.6151044240784322 1 0.05677329701118686
CMWXFREQ_0019 0.0047690763051082925
CMWXSIN_0019 0.690878342154593 1 0.05650274457861476
CMWXCOS_0019 -0.24563258384188566 1 0.05684166683060986
CMWXFREQ_0020 0.005020080321166624
CMWXSIN_0020 -0.8227819896756196 1 0.05625936871827016
CMWXCOS_0020 -0.41486209903118004 1 0.05699996507770999
CMWXFREQ_0021 0.005271084337224955
CMWXSIN_0021 1.0474534912182711 1 0.05693584725924164
CMWXCOS_0021 0.8618797654046676 1 0.05635261418996345
CMWXFREQ_0022 0.005522088353283286
CMWXSIN_0022 -0.34511189858797325 1 0.055801300892170294
CMWXCOS_0022 0.920616971463264 1 0.05755264353434668
CMWXFREQ_0023 0.005773092369341618
CMWXSIN_0023 0.1454118259286491 1 0.05603029726127725
CMWXCOS_0023 0.45055791970256137 1 0.05739756908265092
CMWXFREQ_0024 0.0060240963853999486
CMWXSIN_0024 -0.23847012463655928 1 0.05505646449603128
CMWXCOS_0024 -0.43215665232423606 1 0.058159583976019814
CMWXFREQ_0025 0.00627510040145828
CMWXSIN_0025 -0.352394038392352 1 0.05647865044689391
CMWXCOS_0025 -0.29440801534666206 1 0.056954914678503256
CMWXFREQ_0026 0.006526104417516611
CMWXSIN_0026 -0.09726089542823896 1 0.05648599319460349
CMWXCOS_0026 0.5954256423489478 1 0.05710949890304955
CMWXFREQ_0027 0.006777108433574942
CMWXSIN_0027 -0.48854338113899476 1 0.057227349232574824
CMWXCOS_0027 0.8082914345429778 1 0.05627934748354936
CMWXFREQ_0028 0.007028112449633274
CMWXSIN_0028 0.5005038792418525 1 0.056492455014193914
CMWXCOS_0028 0.7106748205682492 1 0.0570445401821037
CMWXFREQ_0029 0.007279116465691605
CMWXSIN_0029 0.35349970820981613 1 0.05803632935484467
CMWXCOS_0029 -0.33326601250345195 1 0.055685298527841966
CMWXFREQ_0030 0.007530120481749936
CMWXSIN_0030 0.04161269029666645 1 0.055572132931179186
CMWXCOS_0030 0.4494221359830389 1 0.058234016423612124
Estimating the spectral parameters from the CMWaveX fit.
[23]:
# Get the Fourier amplitudes and powers and their uncertainties.
# Note that the `CMWaveX` amplitudes have the units of pc/cm^3/MHz^2.
# We multiply them by a constant factor to convert them to dimensions
# of time so that they are consistent with `PLChromNoise`.
scale = DMconst / 1400**m.TNCHROMIDX.value
idxs = np.array(m2.components["CMWaveX"].get_indices())
a = np.array(
[(scale * m2[f"CMWXSIN_{idx:04d}"].quantity).to_value("s") for idx in idxs]
)
da = np.array(
[(scale * m2[f"CMWXSIN_{idx:04d}"].uncertainty).to_value("s") for idx in idxs]
)
b = np.array(
[(scale * m2[f"CMWXCOS_{idx:04d}"].quantity).to_value("s") for idx in idxs]
)
db = np.array(
[(scale * m2[f"CMWXCOS_{idx:04d}"].uncertainty).to_value("s") for idx in idxs]
)
print(len(idxs))
P = (a**2 + b**2) / 2
dP = ((a * da) ** 2 + (b * db) ** 2) ** 0.5
f0 = (1 / Tspan).to_value(u.Hz)
fyr = (1 / u.year).to_value(u.Hz)
30
[24]:
# We can create a `PLChromNoise` model from the `CMWaveX` model.
# This will estimate the spectral parameters from the `CMWaveX`
# amplitudes.
m3 = plchromnoise_from_cmwavex(m2)
print(m3)
# Created: 2026-08-21T20:07:11.399269
# PINT_version: 1.1.6+68.g7d55691
# User: docs
# Host: build-34177976-project-85767-nanograv-pint
# OS: Linux-7.0.0-1004-aws-x86_64-with-glibc2.35
# Python: 3.11.15 (main, Jun 25 2026, 19:09:59) [GCC 11.4.0]
# Format: pint
# read_time: 2026-08-21T20:07:02.059953
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
# ell1h_shapiro: full
PSR SIM5
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.999999956787824162
FINISH 56985.000000050841412037
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1924.7195246235137
CHI2R 0.9957162569185275
TRES 0.9877706004952963
RAJ 4:59:59.99999955 1 0.00000146238128818959
DECJ 15:00:00.00001146 1 0.00012322699426034855
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.000000000000006134 1 2.8040926342748202197e-14
F1 -1.0000000505979882195e-15 1 6.404970379774501803e-22
PEPOCH 55000.000000000000000000
TZRMJD 55000.000000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.000743291662428148 1 4.349327332289849e-06
TNCHROMAMP -12.985999198590285 0 0.04008384287630784
TNCHROMGAM 3.4394138087142805 0 0.22014757457993323
TNCHROMC 30
PLANET_SHAPIRO N
DM 15.0
CM 1.1040354831558155475 1 0.050400879872884565835
TNCHROMIDX 3.5
[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 0x70103818fd50>
[ ]: