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: 2025-08-22T08:59:17.969043
# PINT_version: 1.1.4+38.g086c46d
# User: docs
# Host: build-29289377-project-85767-nanograv-pint
# OS: Linux-6.8.0-1029-aws-x86_64-with-glibc2.35
# Python: 3.11.12 (main, May 6 2025, 10:45:53) [GCC 11.4.0]
# Format: pint
# read_time: 2025-08-22T08:58:41.635256
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM3
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566689468
FINISH 56985.0000000463735185
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1924.470530011436
CHI2R 0.9803721497765849
TRES 0.98776570295525825173
RAJ 5:00:00.00039125 1 0.00009334338002079630
DECJ 14:59:59.96742676 1 0.01079495150345503689
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.000000000001430925 1 4.315479363735996092e-13
F1 -9.994548189003365421e-16 1 1.7548178693745677244e-19
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 14.999996981863079358 1 4.7068327733927858687e-06
WXEPOCH 55000.0000000000000000
WXFREQ_0001 0.00025100401605860533
WXSIN_0001 4.089504590977486e-06 1 4.840800492721624e-07
WXCOS_0001 -3.023804568303233e-05 1 1.0600167866658755e-05
WXFREQ_0002 0.0005020080321172107
WXSIN_0002 -1.2248870622603246e-06 1 2.4525061182047507e-07
WXCOS_0002 7.035891901402007e-06 1 2.6928162902041242e-06
WXFREQ_0003 0.0007530120481758159
WXSIN_0003 1.118256502910999e-06 1 1.710918541148614e-07
WXCOS_0003 -2.9119353977896626e-06 1 1.2328312668655038e-06
WXFREQ_0004 0.0010040160642344213
WXSIN_0004 -1.6581559262911181e-07 1 1.36565697483751e-07
WXCOS_0004 1.8864033505908323e-06 1 7.258695638060184e-07
WXFREQ_0005 0.0012550200802930267
WXSIN_0005 6.233355512648275e-09 1 1.1834836906901909e-07
WXCOS_0005 -1.6689190476298632e-06 1 4.963202096846019e-07
WXFREQ_0006 0.0015060240963516319
WXSIN_0006 -2.0442724914102732e-07 1 1.107294465177277e-07
WXCOS_0006 1.6303039384825733e-06 1 3.775829004128417e-07
WXFREQ_0007 0.0017570281124102373
WXSIN_0007 6.514333390518089e-07 1 1.0938633911689482e-07
WXCOS_0007 -9.38008455096019e-07 1 3.1448483072908625e-07
WXFREQ_0008 0.0020080321284688426
WXSIN_0008 -3.314792368555245e-07 1 1.177303573953805e-07
WXCOS_0008 6.026077021266014e-07 1 2.9198626098760076e-07
WXFREQ_0009 0.002259036144527448
WXSIN_0009 4.788987043642917e-07 1 1.4489549366400598e-07
WXCOS_0009 -1.0864367293315992e-06 1 3.1394732454655935e-07
WXFREQ_0010 0.0025100401605860534
WXSIN_0010 -8.527288546792026e-07 1 2.466547027375384e-07
WXCOS_0010 1.5006244540114208e-06 1 4.7196345606906166e-07
WXFREQ_0011 0.0027610441766446584
WXSIN_0011 -6.990655384843641e-06 1 2.003801024014765e-06
WXCOS_0011 9.118850823052352e-06 1 3.346565749205455e-06
WXFREQ_0012 0.0030120481927032637
WXSIN_0012 5.356578366341532e-07 1 1.4676031060684668e-07
WXCOS_0012 -5.318516332327679e-07 1 2.0992755860220415e-07
WXFREQ_0013 0.003263052208761869
WXSIN_0013 -2.756745098026658e-07 1 7.134442378602165e-08
WXCOS_0013 2.201028469511223e-07 1 8.53836870762687e-08
WXFREQ_0014 0.0035140562248204745
WXSIN_0014 1.5300073502243007e-07 1 4.916288849727769e-08
WXCOS_0014 -2.712469593394749e-08 1 5.112608308097145e-08
WXFREQ_0015 0.00376506024087908
WXSIN_0015 4.245983787772075e-08 1 4.0001935763848024e-08
WXCOS_0015 8.850805989548093e-08 1 3.960593182629309e-08
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.00023240530042452836 1 0.0008767439654107115
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: 2025-08-22T08:59:17.998307
# PINT_version: 1.1.4+38.g086c46d
# User: docs
# Host: build-29289377-project-85767-nanograv-pint
# OS: Linux-6.8.0-1029-aws-x86_64-with-glibc2.35
# Python: 3.11.12 (main, May 6 2025, 10:45:53) [GCC 11.4.0]
# Format: pint
# read_time: 2025-08-22T08:58:41.635256
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM3
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566689468
FINISH 56985.0000000463735185
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1924.470530011436
CHI2R 0.9803721497765849
TRES 0.98776570295525825173
RAJ 5:00:00.00039125 1 0.00009334338002079630
DECJ 14:59:59.96742676 1 0.01079495150345503689
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.000000000001430925 1 4.315479363735996092e-13
F1 -9.994548189003365421e-16 1 1.7548178693745677244e-19
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 14.999996981863079358 1 4.7068327733927858687e-06
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.00023240530042452836 1 0.0008767439654107115
TNREDAMP -12.351721786140827 0 0.0975807974587293
TNREDGAM 2.8923389517723015 0 0.5519872565223239
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 0x70c4dcf8acd0>
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: 2025-08-22T09:00:00.771576
# PINT_version: 1.1.4+38.g086c46d
# User: docs
# Host: build-29289377-project-85767-nanograv-pint
# OS: Linux-6.8.0-1029-aws-x86_64-with-glibc2.35
# Python: 3.11.12 (main, May 6 2025, 10:45:53) [GCC 11.4.0]
# Format: pint
# read_time: 2025-08-22T08:59:18.296080
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM4
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566281019
FINISH 56985.0000000460385648
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1933.69812908528
CHI2R 1.0003611635205794
TRES 0.99031759072592760517
RAJ 5:00:00.00000174 1 0.00000189556170281422
DECJ 14:59:59.99982095 1 0.00016349563661705439
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.00000000000004888 1 3.674814399416990225e-14
F1 -9.999988654224247173e-16 1 8.387470731055820367e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 14.999997009504412865 1 5.059887230681033215e-06
DMWXEPOCH 55000.0000000000000000
DMWXFREQ_0001 0.00025100401605862414
DMWXSIN_0001 -0.002197629972818859 1 6.098022928458789e-06
DMWXCOS_0001 0.002170963595469702 1 7.037568006517105e-06
DMWXFREQ_0002 0.0005020080321172483
DMWXSIN_0002 -0.00010821404178614692 1 4.857739600941572e-06
DMWXCOS_0002 0.000194049317032426 1 4.601899980098925e-06
DMWXFREQ_0003 0.0007530120481758723
DMWXSIN_0003 -0.00036352474696750456 1 4.652203299116495e-06
DMWXCOS_0003 0.00012476183918775253 1 4.37329956412991e-06
DMWXFREQ_0004 0.0010040160642344966
DMWXSIN_0004 -6.78019617915454e-05 1 4.533120774753405e-06
DMWXCOS_0004 -0.00015447167008887633 1 4.3716915750022654e-06
DMWXFREQ_0005 0.0012550200802931206
DMWXSIN_0005 -0.00019304267994913116 1 4.42010610387566e-06
DMWXCOS_0005 -8.289722827667975e-05 1 4.394626798296678e-06
DMWXFREQ_0006 0.0015060240963517446
DMWXSIN_0006 0.00010224736056874352 1 4.340461969208331e-06
DMWXCOS_0006 0.00010677588830138642 1 4.457659700239173e-06
DMWXFREQ_0007 0.0017570281124103689
DMWXSIN_0007 0.0001076127724391922 1 4.428602069654563e-06
DMWXCOS_0007 2.209606463513009e-05 1 4.362504466569333e-06
DMWXFREQ_0008 0.002008032128468993
DMWXSIN_0008 -6.9702990044011186e-06 1 4.430597629687769e-06
DMWXCOS_0008 8.216136054399586e-06 1 4.339517005042538e-06
DMWXFREQ_0009 0.002259036144527617
DMWXSIN_0009 9.370781109030193e-05 1 4.4598995661442495e-06
DMWXCOS_0009 5.641274434444002e-05 1 4.2836498397716165e-06
DMWXFREQ_0010 0.002510040160586241
DMWXSIN_0010 -1.9137169968107504e-05 1 4.451047810352235e-06
DMWXCOS_0010 1.218553056609287e-05 1 4.340569222955986e-06
DMWXFREQ_0011 0.0027610441766448652
DMWXSIN_0011 -4.722723245885003e-05 1 7.030322034925955e-06
DMWXCOS_0011 -1.2754005244924282e-05 1 7.129613526484698e-06
DMWXFREQ_0012 0.0030120481927034893
DMWXSIN_0012 -3.2104632262240166e-05 1 4.418443567134403e-06
DMWXCOS_0012 -1.3290815790425055e-05 1 4.3681517162921685e-06
DMWXFREQ_0013 0.0032630522087621137
DMWXSIN_0013 1.0488154605334701e-05 1 4.387862762822781e-06
DMWXCOS_0013 -3.648803668387101e-06 1 4.366159395092879e-06
DMWXFREQ_0014 0.0035140562248207378
DMWXSIN_0014 -4.598016972056088e-05 1 4.341085175332989e-06
DMWXCOS_0014 -3.4067925609654e-05 1 4.414679266287372e-06
DMWXFREQ_0015 0.003765060240879362
DMWXSIN_0015 -4.0365873529119465e-05 1 4.4111469713916e-06
DMWXCOS_0015 2.2262390967406637e-05 1 4.3211646925538304e-06
DMWXFREQ_0016 0.004016064256937986
DMWXSIN_0016 2.0292272611604994e-05 1 4.359346841393871e-06
DMWXCOS_0016 2.1037573975441325e-05 1 4.405622530832685e-06
DMWXFREQ_0017 0.00426706827299661
DMWXSIN_0017 -1.2926190895637273e-06 1 4.44763658060513e-06
DMWXCOS_0017 -4.437799530072832e-05 1 4.314303748260531e-06
DMWXFREQ_0018 0.004518072289055234
DMWXSIN_0018 2.7275246038612768e-05 1 4.417784014910885e-06
DMWXCOS_0018 -2.3452170317579923e-05 1 4.352809787075193e-06
DMWXFREQ_0019 0.004769076305113858
DMWXSIN_0019 -1.854621013157318e-05 1 4.344096734285701e-06
DMWXCOS_0019 -1.3071605646776517e-05 1 4.410047214227243e-06
DMWXFREQ_0020 0.005020080321172482
DMWXSIN_0020 1.794847431294142e-05 1 4.351423406725283e-06
DMWXCOS_0020 1.1282608216879772e-05 1 4.39812936333917e-06
DMWXFREQ_0021 0.005271084337231106
DMWXSIN_0021 -1.2628201688511827e-06 1 4.499877910630195e-06
DMWXCOS_0021 8.23435578557348e-06 1 4.251839323549141e-06
DMWXFREQ_0022 0.0055220883532897305
DMWXSIN_0022 1.6399437999658405e-06 1 4.48118788535825e-06
DMWXCOS_0022 1.2746209811234782e-05 1 4.290623293826548e-06
DMWXFREQ_0023 0.0057730923693483545
DMWXSIN_0023 -1.0775800052453385e-05 1 4.406171961221124e-06
DMWXCOS_0023 1.0988709814495701e-05 1 4.368057610061865e-06
DMWXFREQ_0024 0.0060240963854069785
DMWXSIN_0024 2.128074401879036e-05 1 4.479599323546321e-06
DMWXCOS_0024 1.1548539253300997e-06 1 4.286057140969615e-06
DMWXFREQ_0025 0.006275100401465603
DMWXSIN_0025 -2.4195184469907023e-05 1 4.288866512228133e-06
DMWXCOS_0025 6.12485061755461e-07 1 4.477555113081536e-06
DMWXFREQ_0026 0.0065261044175242275
DMWXSIN_0026 6.6550465836663e-06 1 4.437368825068553e-06
DMWXCOS_0026 1.5342181087893907e-05 1 4.318368074597448e-06
DMWXFREQ_0027 0.0067771084335828515
DMWXSIN_0027 1.6041204873870308e-05 1 4.403005350888904e-06
DMWXCOS_0027 -1.791439295841358e-06 1 4.357639221483473e-06
DMWXFREQ_0028 0.0070281124496414755
DMWXSIN_0028 1.0994832527105089e-05 1 4.4038954025388726e-06
DMWXCOS_0028 2.5969963739524082e-06 1 4.3586423057983376e-06
DMWXFREQ_0029 0.0072791164657000995
DMWXSIN_0029 1.1673303008608406e-05 1 4.353322993867019e-06
DMWXCOS_0029 -7.68057250668106e-06 1 4.401707311616152e-06
DMWXFREQ_0030 0.007530120481758724
DMWXSIN_0030 5.991951563192312e-06 1 4.451044948214023e-06
DMWXCOS_0030 -1.0023098109008409e-05 1 4.310951758307237e-06
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.0005115072689343497 1 5.57673495052251e-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: 2025-08-22T09:00:00.811195
# PINT_version: 1.1.4+38.g086c46d
# User: docs
# Host: build-29289377-project-85767-nanograv-pint
# OS: Linux-6.8.0-1029-aws-x86_64-with-glibc2.35
# Python: 3.11.12 (main, May 6 2025, 10:45:53) [GCC 11.4.0]
# Format: pint
# read_time: 2025-08-22T08:59:18.296080
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM4
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566281019
FINISH 56985.0000000460385648
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1933.69812908528
CHI2R 1.0003611635205794
TRES 0.99031759072592760517
RAJ 5:00:00.00000174 1 0.00000189556170281422
DECJ 14:59:59.99982095 1 0.00016349563661705439
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.00000000000004888 1 3.674814399416990225e-14
F1 -9.999988654224247173e-16 1 8.387470731055820367e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 14.999997009504412865 1 5.059887230681033215e-06
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.0005115072689343497 1 5.57673495052251e-06
TNDMAMP -13.01805312819021 0 0.04213582744483032
TNDMGAM 3.047403361425293 0 0.2125275853525287
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 0x70c4dd677710>
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: 2025-08-22T09:00:08.361817
# PINT_version: 1.1.4+38.g086c46d
# User: docs
# Host: build-29289377-project-85767-nanograv-pint
# OS: Linux-6.8.0-1029-aws-x86_64-with-glibc2.35
# Python: 3.11.12 (main, May 6 2025, 10:45:53) [GCC 11.4.0]
# Format: pint
# read_time: 2025-08-22T09:00:01.221488
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM5
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567218403
FINISH 56985.0000000474333218
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1886.171478890896
CHI2R 0.9757741742839607
TRES 0.9646633082120204804
RAJ 4:59:59.99999995 1 0.00000140584319131311
DECJ 14:59:59.99994905 1 0.00012297273739955955
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.000000000000012414 1 2.8036533961911535114e-14
F1 -9.999992561770502888e-16 1 6.273047845942850402e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.0
CM 1.2278679202487329859 1 0.050308085868335149227
TNCHROMIDX 3.5
CMWXEPOCH 55000.0000000000000000
CMWXFREQ_0001 0.00025100401605854207
CMWXSIN_0001 -1.3496895265093807 1 0.06532580196917566
CMWXCOS_0001 -48.008288728547434 1 0.07108169995220924
CMWXFREQ_0002 0.0005020080321170841
CMWXSIN_0002 14.76785045662619 1 0.058852620042117436
CMWXCOS_0002 16.09277673623306 1 0.05932191972440692
CMWXFREQ_0003 0.0007530120481756262
CMWXSIN_0003 -3.8123133448939974 1 0.059903109740208495
CMWXCOS_0003 11.42916597211265 1 0.05670858578524261
CMWXFREQ_0004 0.0010040160642341683
CMWXSIN_0004 -0.9843471238681605 1 0.05917530989972974
CMWXCOS_0004 -4.450832437120306 1 0.056896621953557394
CMWXFREQ_0005 0.0012550200802927103
CMWXSIN_0005 2.5018940924958604 1 0.059306361974963986
CMWXCOS_0005 -5.980760421142669 1 0.05633947903230179
CMWXFREQ_0006 0.0015060240963512524
CMWXSIN_0006 5.175109403138978 1 0.05737224252288993
CMWXCOS_0006 5.955435674910116 1 0.05802675826330086
CMWXFREQ_0007 0.0017570281124097945
CMWXSIN_0007 -4.118120383659263 1 0.05720431676270974
CMWXCOS_0007 -1.3091609674269582 1 0.05802105194763716
CMWXFREQ_0008 0.0020080321284683365
CMWXSIN_0008 -1.5341453270709704 1 0.057132743873035956
CMWXCOS_0008 0.976159126776136 1 0.05774394605089907
CMWXFREQ_0009 0.0022590361445268786
CMWXSIN_0009 -0.7163882552295636 1 0.058044950030487534
CMWXCOS_0009 5.105562107223552 1 0.05690860859621168
CMWXFREQ_0010 0.0025100401605854207
CMWXSIN_0010 4.328340451890255 1 0.056038282937681906
CMWXCOS_0010 0.06742594178159615 1 0.05903873393647229
CMWXFREQ_0011 0.0027610441766439627
CMWXSIN_0011 1.7330734574314364 1 0.07142787434597471
CMWXCOS_0011 2.753241144454969 1 0.07024994492290826
CMWXFREQ_0012 0.003012048192702505
CMWXSIN_0012 -0.07367440846667538 1 0.058111639382278577
CMWXCOS_0012 -0.004203989285655284 1 0.0571187320195523
CMWXFREQ_0013 0.003263052208761047
CMWXSIN_0013 0.5739933461895095 1 0.05917568328053421
CMWXCOS_0013 -0.017949705051224216 1 0.055970781758711195
CMWXFREQ_0014 0.003514056224819589
CMWXSIN_0014 -0.8535256776045742 1 0.0576298619644878
CMWXCOS_0014 2.083722360024269 1 0.05719901680501403
CMWXFREQ_0015 0.003765060240878131
CMWXSIN_0015 -1.3407692531918796 1 0.05735735422696049
CMWXCOS_0015 0.25202332496585905 1 0.05743604776634145
CMWXFREQ_0016 0.004016064256936673
CMWXSIN_0016 -0.3266001522135103 1 0.056892572580679125
CMWXCOS_0016 1.140837522179714 1 0.05790774276426341
CMWXFREQ_0017 0.004267068272995215
CMWXSIN_0017 1.2503959092188994 1 0.057253602250343344
CMWXCOS_0017 -1.5726679816222613 1 0.05756668243490518
CMWXFREQ_0018 0.004518072289053757
CMWXSIN_0018 -0.36927371090510436 1 0.05725864342809286
CMWXCOS_0018 -0.61173825766449 1 0.05758641745945116
CMWXFREQ_0019 0.004769076305112299
CMWXSIN_0019 -0.8121106072992426 1 0.05803907519422637
CMWXCOS_0019 0.12386904911646283 1 0.05687465732842557
CMWXFREQ_0020 0.005020080321170841
CMWXSIN_0020 -0.3794532210215655 1 0.05753525545083635
CMWXCOS_0020 0.19606709209921766 1 0.05687545000788558
CMWXFREQ_0021 0.005271084337229383
CMWXSIN_0021 -0.39400204446603976 1 0.05683204162297632
CMWXCOS_0021 -0.588221052875795 1 0.057553034017263374
CMWXFREQ_0022 0.0055220883532879255
CMWXSIN_0022 0.894560980027318 1 0.05722798504030763
CMWXCOS_0022 0.8515742045217869 1 0.05712830333113069
CMWXFREQ_0023 0.005773092369346467
CMWXSIN_0023 0.14312991002890732 1 0.0566246927668895
CMWXCOS_0023 1.1095188504400288 1 0.058054190216531514
CMWXFREQ_0024 0.00602409638540501
CMWXSIN_0024 0.44498207479111995 1 0.0580807893095804
CMWXCOS_0024 -0.2776162791590537 1 0.056428978631417714
CMWXFREQ_0025 0.006275100401463551
CMWXSIN_0025 0.005287950795694433 1 0.056088877751159616
CMWXCOS_0025 0.3426403857226765 1 0.05858837458247478
CMWXFREQ_0026 0.006526104417522094
CMWXSIN_0026 0.2021905913580394 1 0.05850368742391591
CMWXCOS_0026 -0.5477982371535791 1 0.05618643179134904
CMWXFREQ_0027 0.006777108433580635
CMWXSIN_0027 -0.27042943918876644 1 0.05763166152944035
CMWXCOS_0027 0.23347856473162076 1 0.057064105036480635
CMWXFREQ_0028 0.007028112449639178
CMWXSIN_0028 0.28338401087609966 1 0.05749266352501264
CMWXCOS_0028 0.4288780126963341 1 0.05724655753313824
CMWXFREQ_0029 0.0072791164656977195
CMWXSIN_0029 0.3848480114484888 1 0.057612566836913154
CMWXCOS_0029 0.1500677771440271 1 0.05698078219090805
CMWXFREQ_0030 0.007530120481756262
CMWXSIN_0030 0.6358836684784619 1 0.05834341090165006
CMWXCOS_0030 0.15557422051645609 1 0.05632636230095492
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -5.880596721184931e-05 1 4.107496945141918e-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: 2025-08-22T09:00:08.400140
# PINT_version: 1.1.4+38.g086c46d
# User: docs
# Host: build-29289377-project-85767-nanograv-pint
# OS: Linux-6.8.0-1029-aws-x86_64-with-glibc2.35
# Python: 3.11.12 (main, May 6 2025, 10:45:53) [GCC 11.4.0]
# Format: pint
# read_time: 2025-08-22T09:00:01.221488
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM5
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567218403
FINISH 56985.0000000474333218
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1886.171478890896
CHI2R 0.9757741742839607
TRES 0.9646633082120204804
RAJ 4:59:59.99999995 1 0.00000140584319131311
DECJ 14:59:59.99994905 1 0.00012297273739955955
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.000000000000012414 1 2.8036533961911535114e-14
F1 -9.999992561770502888e-16 1 6.273047845942850402e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.0
CM 1.2278679202487329859 1 0.050308085868335149227
TNCHROMIDX 3.5
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -5.880596721184931e-05 1 4.107496945141918e-06
TNCHROMAMP -13.120694254601394 0 0.04018401065208207
TNCHROMGAM 2.807318348209132 0 0.2593840949151346
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 0x70c4ced00950>
[ ]: