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-04-25T19:13:42.992733
# PINT_version: 1.1.3
# User: docs
# Host: build-27973916-project-85767-nanograv-pint
# OS: Linux-6.8.0-1021-aws-x86_64-with-glibc2.35
# Python: 3.11.10 (main, Nov 5 2024, 16:13:19) [GCC 11.4.0]
# Format: pint
# read_time: 2025-04-25T19:13:06.797993
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM3
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999568564121
FINISH 56985.0000000466048843
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1976.2487353715876
CHI2R 1.0067492284114048
TRES 0.98076659973533944877
RAJ 4:59:59.99990117 1 0.00009439613531363219
DECJ 15:00:00.02266803 1 0.01073789928559386320
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.00000000000026429 1 4.4819105808922956456e-13
F1 -1.0004490525657857629e-15 1 1.7551099430937645479e-19
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.000002815362470379 1 4.6212567333903515116e-06
WXEPOCH 55000.0000000000000000
WXFREQ_0001 0.00025100401605860257
WXSIN_0001 3.054360650828639e-06 1 5.017540429333399e-07
WXCOS_0001 1.3493751948023366e-05 1 1.059981667362874e-05
WXFREQ_0002 0.0005020080321172051
WXSIN_0002 -1.4190733692127331e-06 1 2.5429396352700877e-07
WXCOS_0002 -6.92846836918644e-06 1 2.6914495576043164e-06
WXFREQ_0003 0.0007530120481758077
WXSIN_0003 -7.462905501423177e-07 1 1.7664107262971507e-07
WXCOS_0003 3.443404168626045e-06 1 1.2317094758295847e-06
WXFREQ_0004 0.0010040160642344103
WXSIN_0004 4.167819821197001e-07 1 1.4042181158382986e-07
WXCOS_0004 -2.050898709407139e-06 1 7.237727991740395e-07
WXFREQ_0005 0.0012550200802930128
WXSIN_0005 2.853975712494987e-07 1 1.2240056699590148e-07
WXCOS_0005 1.5056183419026351e-06 1 4.935774043444172e-07
WXFREQ_0006 0.0015060240963516154
WXSIN_0006 -6.387285580091275e-08 1 1.13360336137856e-07
WXCOS_0006 -3.89064544947718e-07 1 3.736451831043143e-07
WXFREQ_0007 0.001757028112410218
WXSIN_0007 -2.5320802886009586e-07 1 1.1101745231122115e-07
WXCOS_0007 8.170095337962644e-07 1 3.1264749651414255e-07
WXFREQ_0008 0.0020080321284688205
WXSIN_0008 -2.3902703871678225e-08 1 1.193373171241267e-07
WXCOS_0008 -6.441700485263734e-07 1 2.9016270330149514e-07
WXFREQ_0009 0.0022590361445274233
WXSIN_0009 9.56588303270652e-08 1 1.4641125671892112e-07
WXCOS_0009 7.79666356050766e-07 1 3.1136310400788766e-07
WXFREQ_0010 0.0025100401605860257
WXSIN_0010 -1.0914543947815197e-07 1 2.499135592837588e-07
WXCOS_0010 -1.2690455527125327e-06 1 4.667854272728541e-07
WXFREQ_0011 0.0027610441766446284
WXSIN_0011 -6.530015747887636e-07 1 2.015375799165559e-06
WXCOS_0011 -7.3643829791897535e-06 1 3.319781272371602e-06
WXFREQ_0012 0.003012048192703231
WXSIN_0012 1.458991431767975e-07 1 1.463012511798286e-07
WXCOS_0012 5.818081144971718e-07 1 2.0892458360646624e-07
WXFREQ_0013 0.0032630522087618336
WXSIN_0013 -6.284650049850665e-08 1 7.078372293872416e-08
WXCOS_0013 -1.2933801538385428e-07 1 8.605439831054999e-08
WXFREQ_0014 0.003514056224820436
WXSIN_0014 -8.55663825210462e-09 1 4.811005968952105e-08
WXCOS_0014 -4.9537140886517964e-08 1 5.052625462500755e-08
WXFREQ_0015 0.0037650602408790387
WXSIN_0015 -1.1665186284655846e-07 1 3.919844582919343e-08
WXCOS_0015 3.218571711559628e-08 1 3.903480570429775e-08
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.001287105284327685 1 0.0008765017463313435
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-04-25T19:13:43.021112
# PINT_version: 1.1.3
# User: docs
# Host: build-27973916-project-85767-nanograv-pint
# OS: Linux-6.8.0-1021-aws-x86_64-with-glibc2.35
# Python: 3.11.10 (main, Nov 5 2024, 16:13:19) [GCC 11.4.0]
# Format: pint
# read_time: 2025-04-25T19:13:06.797993
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM3
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999568564121
FINISH 56985.0000000466048843
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1976.2487353715876
CHI2R 1.0067492284114048
TRES 0.98076659973533944877
RAJ 4:59:59.99990117 1 0.00009439613531363219
DECJ 15:00:00.02266803 1 0.01073789928559386320
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.00000000000026429 1 4.4819105808922956456e-13
F1 -1.0004490525657857629e-15 1 1.7551099430937645479e-19
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.000002815362470379 1 4.6212567333903515116e-06
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.001287105284327685 1 0.0008765017463313435
TNREDAMP -12.754846796174277 0 0.15985219666642375
TNREDGAM 3.5453335104692547 0 0.713657440956131
TNREDC 15.0
[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 0x7847c5450850>
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: 2025-04-25T19:14:24.909100
# PINT_version: 1.1.3
# User: docs
# Host: build-27973916-project-85767-nanograv-pint
# OS: Linux-6.8.0-1021-aws-x86_64-with-glibc2.35
# Python: 3.11.10 (main, Nov 5 2024, 16:13:19) [GCC 11.4.0]
# Format: pint
# read_time: 2025-04-25T19:13:43.424876
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM4
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567362847
FINISH 56985.0000000470324769
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1988.0673756948056
CHI2R 1.0253055057734943
TRES 1.0104911616518964136
RAJ 5:00:00.00000269 1 0.00000188995862122353
DECJ 15:00:00.00002903 1 0.00016632249092022995
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999998996 1 3.713480229936589964e-14
F1 -9.99999221191338389e-16 1 8.522448293184007202e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.0000038461467818185 1 4.993013686443469371e-06
DMWXEPOCH 55000.0000000000000000
DMWXFREQ_0001 0.0002510040160585682
DMWXSIN_0001 -0.0013989478080285213 1 6.079678815169478e-06
DMWXCOS_0001 -0.0022767092064315004 1 6.945307321478476e-06
DMWXFREQ_0002 0.0005020080321171364
DMWXSIN_0002 0.0005649108917920255 1 4.8126321631943316e-06
DMWXCOS_0002 0.00022444207885372576 1 4.514621423666806e-06
DMWXFREQ_0003 0.0007530120481757046
DMWXSIN_0003 0.0005241054090921158 1 4.5429848962227915e-06
DMWXCOS_0003 -0.0003381117027895699 1 4.352599445546296e-06
DMWXFREQ_0004 0.0010040160642342728
DMWXSIN_0004 0.00023224543996598634 1 4.440991115865019e-06
DMWXCOS_0004 0.00010840471872800818 1 4.320025931937416e-06
DMWXFREQ_0005 0.0012550200802928409
DMWXSIN_0005 -0.00013567140321968453 1 4.509694968486752e-06
DMWXCOS_0005 -0.000275844573014689 1 4.229226957304002e-06
DMWXFREQ_0006 0.0015060240963514092
DMWXSIN_0006 0.00010130295358760314 1 4.423497892063346e-06
DMWXCOS_0006 -1.4230720801793442e-05 1 4.290316741987764e-06
DMWXFREQ_0007 0.0017570281124099773
DMWXSIN_0007 6.474806540169774e-05 1 4.3871249633577485e-06
DMWXCOS_0007 0.00010538146164091029 1 4.299692735009792e-06
DMWXFREQ_0008 0.0020080321284685456
DMWXSIN_0008 -7.971946361676778e-06 1 4.274553609538613e-06
DMWXCOS_0008 -2.783652142786104e-06 1 4.39748342535812e-06
DMWXFREQ_0009 0.0022590361445271137
DMWXSIN_0009 -3.9827911276379846e-05 1 4.356948466034298e-06
DMWXCOS_0009 -1.6818424259214715e-05 1 4.3254165605035866e-06
DMWXFREQ_0010 0.0025100401605856817
DMWXSIN_0010 -3.506631677816747e-05 1 4.354675829767269e-06
DMWXCOS_0010 -3.8081632086391444e-05 1 4.388586769873062e-06
DMWXFREQ_0011 0.0027610441766442503
DMWXSIN_0011 5.508618394676314e-05 1 7.029024601309757e-06
DMWXCOS_0011 -2.3939204303059854e-05 1 7.018802614414904e-06
DMWXFREQ_0012 0.0030120481927028184
DMWXSIN_0012 -2.420354753388627e-05 1 4.283759938020829e-06
DMWXCOS_0012 -8.511124293828144e-05 1 4.416953835429475e-06
DMWXFREQ_0013 0.0032630522087613864
DMWXSIN_0013 3.421901147512207e-06 1 4.427570997524464e-06
DMWXCOS_0013 1.4478436817401263e-05 1 4.235157646224127e-06
DMWXFREQ_0014 0.0035140562248199545
DMWXSIN_0014 -1.8957590787785936e-05 1 4.448515391659194e-06
DMWXCOS_0014 1.529928001329849e-05 1 4.218750413418105e-06
DMWXFREQ_0015 0.003765060240878523
DMWXSIN_0015 9.42736576853696e-06 1 4.331151020212839e-06
DMWXCOS_0015 -8.987953430616455e-06 1 4.351089718537813e-06
DMWXFREQ_0016 0.004016064256937091
DMWXSIN_0016 -5.5543977219161385e-05 1 4.262469746301342e-06
DMWXCOS_0016 3.6822695613169156e-06 1 4.4142191991440844e-06
DMWXFREQ_0017 0.00426706827299566
DMWXSIN_0017 6.537998577290224e-05 1 4.333810364324099e-06
DMWXCOS_0017 1.5177204458247716e-05 1 4.348366182862225e-06
DMWXFREQ_0018 0.004518072289054227
DMWXSIN_0018 -8.597416765250975e-06 1 4.285024019505573e-06
DMWXCOS_0018 -1.096727724225024e-05 1 4.3875379029727474e-06
DMWXFREQ_0019 0.004769076305112796
DMWXSIN_0019 -2.7953414177332615e-06 1 4.3530566484747266e-06
DMWXCOS_0019 -7.151968710079922e-06 1 4.320225700585614e-06
DMWXFREQ_0020 0.0050200803211713635
DMWXSIN_0020 5.164507116897369e-06 1 4.382374257099054e-06
DMWXCOS_0020 2.927334961909498e-06 1 4.288641055217795e-06
DMWXFREQ_0021 0.005271084337229932
DMWXSIN_0021 7.3091977910654735e-06 1 4.287700971001209e-06
DMWXCOS_0021 -6.0085397935665685e-06 1 4.38408411136867e-06
DMWXFREQ_0022 0.0055220883532885005
DMWXSIN_0022 -1.9422889475773494e-05 1 4.335461666156039e-06
DMWXCOS_0022 -1.2061649071132232e-06 1 4.33585907136291e-06
DMWXFREQ_0023 0.005773092369347068
DMWXSIN_0023 5.938441547404043e-06 1 4.3279942050491516e-06
DMWXCOS_0023 1.9584130965625197e-05 1 4.335435385599081e-06
DMWXFREQ_0024 0.006024096385405637
DMWXSIN_0024 1.5827332155931208e-05 1 4.323041182883758e-06
DMWXCOS_0024 2.3586681511780647e-06 1 4.333118809462631e-06
DMWXFREQ_0025 0.006275100401464205
DMWXSIN_0025 -1.1968531716558146e-05 1 4.292393536244339e-06
DMWXCOS_0025 2.0216226482832047e-07 1 4.3745046942174305e-06
DMWXFREQ_0026 0.006526104417522773
DMWXSIN_0026 -4.4916156274207886e-06 1 4.387840342311987e-06
DMWXCOS_0026 -1.3506724957877616e-05 1 4.2740581547692735e-06
DMWXFREQ_0027 0.006777108433581341
DMWXSIN_0027 -8.863509600362029e-06 1 4.278936714099419e-06
DMWXCOS_0027 1.9753389098061763e-05 1 4.3839822233807365e-06
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.0005436342472844697 1 5.6125915473997405e-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: 2025-04-25T19:14:24.947219
# PINT_version: 1.1.3
# User: docs
# Host: build-27973916-project-85767-nanograv-pint
# OS: Linux-6.8.0-1021-aws-x86_64-with-glibc2.35
# Python: 3.11.10 (main, Nov 5 2024, 16:13:19) [GCC 11.4.0]
# Format: pint
# read_time: 2025-04-25T19:13:43.424876
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM4
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567362847
FINISH 56985.0000000470324769
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1988.0673756948056
CHI2R 1.0253055057734943
TRES 1.0104911616518964136
RAJ 5:00:00.00000269 1 0.00000188995862122353
DECJ 15:00:00.00002903 1 0.00016632249092022995
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999998996 1 3.713480229936589964e-14
F1 -9.99999221191338389e-16 1 8.522448293184007202e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.0000038461467818185 1 4.993013686443469371e-06
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.0005436342472844697 1 5.6125915473997405e-06
TNDMAMP -12.987210336411174 0 0.04367894922830129
TNDMGAM 3.169867322937473 0 0.24290002480366138
TNDMC 27.0
[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 0x7847b8082a90>
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-04-25T19:14:32.875967
# PINT_version: 1.1.3
# User: docs
# Host: build-27973916-project-85767-nanograv-pint
# OS: Linux-6.8.0-1021-aws-x86_64-with-glibc2.35
# Python: 3.11.10 (main, Nov 5 2024, 16:13:19) [GCC 11.4.0]
# Format: pint
# read_time: 2025-04-25T19:14:25.477311
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM5
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567653473
FINISH 56985.0000000495860532
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 2017.0695480597647
CHI2R 1.0434917475735979
TRES 1.0102796249802010473
RAJ 4:59:59.99999937 1 0.00000145047232593044
DECJ 15:00:00.00009919 1 0.00012360638394895532
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.00000000000004686 1 2.8114953931020545323e-14
F1 -1.0000002371994813805e-15 1 6.382797738281504891e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.0
CM 1.1285158512690496573 1 0.049822983015583417743
TNCHROMIDX 3.5
CMWXEPOCH 55000.0000000000000000
CMWXFREQ_0001 0.00025100401605840914
CMWXSIN_0001 83.37798450642805 1 0.0657067593934502
CMWXCOS_0001 -164.06904011784508 1 0.06836945323202252
CMWXFREQ_0002 0.0005020080321168183
CMWXSIN_0002 2.5668561591654173 1 0.059231814449480176
CMWXCOS_0002 31.020329020243356 1 0.0567995426655353
CMWXFREQ_0003 0.0007530120481752274
CMWXSIN_0003 -40.28378558599709 1 0.058254118093319145
CMWXCOS_0003 13.753346201242964 1 0.055744250243425233
CMWXFREQ_0004 0.0010040160642336366
CMWXSIN_0004 0.8255012524669064 1 0.05611995886752616
CMWXCOS_0004 -12.710257550922279 1 0.05756558224188158
CMWXFREQ_0005 0.0012550200802920457
CMWXSIN_0005 -0.16268203508220155 1 0.055477721132935
CMWXCOS_0005 -13.56838635712535 1 0.05783943671518137
CMWXFREQ_0006 0.0015060240963504549
CMWXSIN_0006 8.392306596631409 1 0.0581512138064508
CMWXCOS_0006 5.932754155842204 1 0.05503338628032997
CMWXFREQ_0007 0.0017570281124088638
CMWXSIN_0007 1.3870163340664057 1 0.05527183031225447
CMWXCOS_0007 5.589591070067888 1 0.057887519752226044
CMWXFREQ_0008 0.002008032128467273
CMWXSIN_0008 -2.150174631329295 1 0.05697716387695915
CMWXCOS_0008 -1.524906870224698 1 0.056220336385828414
CMWXFREQ_0009 0.002259036144525682
CMWXSIN_0009 3.3387319667245143 1 0.057251062620511585
CMWXCOS_0009 5.9000450496295 1 0.05598746892658823
CMWXFREQ_0010 0.0025100401605840914
CMWXSIN_0010 -3.134291058102196 1 0.05458619876169419
CMWXCOS_0010 -0.4212742110010921 1 0.058839808557841
CMWXFREQ_0011 0.0027610441766425004
CMWXSIN_0011 1.4744156618205675 1 0.0707753652722478
CMWXCOS_0011 -2.61630905779263 1 0.0693593535286471
CMWXFREQ_0012 0.0030120481927009097
CMWXSIN_0012 -1.0463574858265463 1 0.05504872937021639
CMWXCOS_0012 1.660327208652486 1 0.0579540771778797
CMWXFREQ_0013 0.0032630522087593187
CMWXSIN_0013 0.30677414292251765 1 0.057889043562484546
CMWXCOS_0013 0.4418585349959677 1 0.05495026702076665
CMWXFREQ_0014 0.0035140562248177276
CMWXSIN_0014 -2.626625399204619 1 0.05764291368890012
CMWXCOS_0014 -1.2152154488269662 1 0.05544910032869357
CMWXFREQ_0015 0.003765060240876137
CMWXSIN_0015 0.9588676767470993 1 0.05767539121054689
CMWXCOS_0015 -1.3286966604308792 1 0.0553232338142267
CMWXFREQ_0016 0.004016064256934546
CMWXSIN_0016 -0.25450589943131446 1 0.05494530633850789
CMWXCOS_0016 -1.323232777245967 1 0.05802606295513827
CMWXFREQ_0017 0.004267068272992955
CMWXSIN_0017 0.1605430238696972 1 0.05768032013855451
CMWXCOS_0017 -0.5668176208976322 1 0.05520099682868873
CMWXFREQ_0018 0.004518072289051364
CMWXSIN_0018 -0.7702832028396797 1 0.057181346570356284
CMWXCOS_0018 -0.03417687856740607 1 0.055642747493785644
CMWXFREQ_0019 0.004769076305109773
CMWXSIN_0019 0.7277055506144118 1 0.056689735432128574
CMWXCOS_0019 0.04304504942796155 1 0.056157066683740375
CMWXFREQ_0020 0.005020080321168183
CMWXSIN_0020 -1.3206706263560286 1 0.054772964329660084
CMWXCOS_0020 0.8358995829113862 1 0.05798566827884903
CMWXFREQ_0021 0.005271084337226592
CMWXSIN_0021 -0.6429839439117915 1 0.05580628917259059
CMWXCOS_0021 -0.23173665275156025 1 0.05702893283365541
CMWXFREQ_0022 0.005522088353285001
CMWXSIN_0022 -0.9222126917109118 1 0.05711309664237067
CMWXCOS_0022 1.478538970881841 1 0.055612364303876595
CMWXFREQ_0023 0.00577309236934341
CMWXSIN_0023 -0.5186515612351832 1 0.057160940452283704
CMWXCOS_0023 -1.1338222467076449 1 0.05554570951245871
CMWXFREQ_0024 0.0060240963854018195
CMWXSIN_0024 1.3871201227021281 1 0.05516514823719787
CMWXCOS_0024 -0.25156054534254124 1 0.05766392093906659
CMWXFREQ_0025 0.006275100401460228
CMWXSIN_0025 1.6744178836270205 1 0.057037677583923765
CMWXCOS_0025 0.66972976639629 1 0.05590216076173842
CMWXFREQ_0026 0.006526104417518637
CMWXSIN_0026 -0.12105220389506248 1 0.05640950583359847
CMWXCOS_0026 0.626237086615664 1 0.05662146741699743
CMWXFREQ_0027 0.006777108433577046
CMWXSIN_0027 0.04399724449603705 1 0.05712600076577524
CMWXCOS_0027 -0.297964579456211 1 0.05589178338062079
CMWXFREQ_0028 0.007028112449635455
CMWXSIN_0028 0.5954096651122338 1 0.05658687240648319
CMWXCOS_0028 0.367240431431695 1 0.05637539865280734
CMWXFREQ_0029 0.007279116465693865
CMWXSIN_0029 -0.47968385825476034 1 0.0573681243616053
CMWXCOS_0029 0.48300521858572315 1 0.055768912471678624
CMWXFREQ_0030 0.007530120481752274
CMWXSIN_0030 0.315262591207279 1 0.05550689913636328
CMWXCOS_0030 0.3062030482664468 1 0.057560113761878584
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.0005294421856029251 1 4.2796726743750345e-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-04-25T19:14:32.914360
# PINT_version: 1.1.3
# User: docs
# Host: build-27973916-project-85767-nanograv-pint
# OS: Linux-6.8.0-1021-aws-x86_64-with-glibc2.35
# Python: 3.11.10 (main, Nov 5 2024, 16:13:19) [GCC 11.4.0]
# Format: pint
# read_time: 2025-04-25T19:14:25.477311
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM5
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567653473
FINISH 56985.0000000495860532
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 2017.0695480597647
CHI2R 1.0434917475735979
TRES 1.0102796249802010473
RAJ 4:59:59.99999937 1 0.00000145047232593044
DECJ 15:00:00.00009919 1 0.00012360638394895532
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.00000000000004686 1 2.8114953931020545323e-14
F1 -1.0000002371994813805e-15 1 6.382797738281504891e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.0
CM 1.1285158512690496573 1 0.049822983015583417743
TNCHROMIDX 3.5
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.0005294421856029251 1 4.2796726743750345e-06
TNCHROMAMP -12.950768163728512 0 0.04001848951602807
TNCHROMGAM 3.2280205046946087 0 0.20594058936415496
TNCHROMC 30.0
[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 0x7847c43e66d0>
[ ]: