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 = 13
[5]:
print(np.argmin(d_aics))
13
[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-01-17T17:58:30.430611
# PINT_version: 1.1.1+21.g5e4d4a3
# User: docs
# Host: build-26886578-project-85767-nanograv-pint
# OS: Linux-5.19.0-1028-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-01-17T17:57:42.347245
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM3
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567793287
FINISH 56985.0000000465485648
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 2052.14475644944
CHI2R 1.0432866072442502
TRES 1.0204795556440334719
RAJ 5:00:00.00010106 1 0.00006205289580283582
DECJ 15:00:00.00813585 1 0.00725206049468805122
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.00000000000096198 1 3.1544984608191554236e-13
F1 -1.0001287029299695832e-15 1 1.3909954744010472733e-19
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 14.999996262588832235 1 4.793012518440144986e-06
WXEPOCH 55000.0000000000000000
WXFREQ_0001 0.0002510040160586012
WXSIN_0001 7.06458888846353e-06 1 3.537235307858139e-07
WXCOS_0001 -5.416243466619228e-07 1 8.393786758866322e-06
WXFREQ_0002 0.0005020080321172024
WXSIN_0002 -5.963226644127301e-07 1 1.7949528082611834e-07
WXCOS_0002 -2.0725744889615108e-06 1 2.1268412788650935e-06
WXFREQ_0003 0.0007530120481758036
WXSIN_0003 5.109470121933902e-07 1 1.2535157595601008e-07
WXCOS_0003 7.924912740217189e-07 1 9.700797131422196e-07
WXFREQ_0004 0.0010040160642344048
WXSIN_0004 -2.801738270713388e-07 1 1.0055020320884033e-07
WXCOS_0004 -4.843879627093389e-07 1 5.662938090123987e-07
WXFREQ_0005 0.0012550200802930059
WXSIN_0005 9.438798001809711e-07 1 8.698749281740286e-08
WXCOS_0005 7.113030388274789e-07 1 3.836274347897225e-07
WXFREQ_0006 0.0015060240963516072
WXSIN_0006 -2.807091577639378e-07 1 8.095786611618818e-08
WXCOS_0006 1.3735854333079267e-07 1 2.8798021527049764e-07
WXFREQ_0007 0.0017570281124102084
WXSIN_0007 1.673517750350353e-07 1 7.900917341797437e-08
WXCOS_0007 3.6700373997919873e-07 1 2.3789425892902993e-07
WXFREQ_0008 0.0020080321284688097
WXSIN_0008 -2.414608377857738e-07 1 8.267119330741593e-08
WXCOS_0008 -2.594374501481749e-07 1 2.1622097963442515e-07
WXFREQ_0009 0.0022590361445274107
WXSIN_0009 3.6867236104821674e-07 1 9.911834851885475e-08
WXCOS_0009 3.5801498762999547e-07 1 2.2707632371233543e-07
WXFREQ_0010 0.0025100401605860118
WXSIN_0010 -4.686362495973888e-07 1 1.62799971949175e-07
WXCOS_0010 -4.281925238385445e-07 1 3.296940372574073e-07
WXFREQ_0011 0.0027610441766446133
WXSIN_0011 -4.902164990137589e-06 1 1.2758844615024466e-06
WXCOS_0011 -3.176497804212741e-06 1 2.230515569269665e-06
WXFREQ_0012 0.0030120481927032143
WXSIN_0012 3.464149326047503e-07 1 9.418223532049626e-08
WXCOS_0012 2.809314400341747e-07 1 1.3484745130657152e-07
WXFREQ_0013 0.003263052208761816
WXSIN_0013 -1.8480995805763308e-07 1 4.967295942032303e-08
WXCOS_0013 -1.596292992676261e-07 1 5.752271501152903e-08
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.0007298067796714427 1 0.0006935093127328683
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)
13
[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-01-17T17:58:30.465436
# PINT_version: 1.1.1+21.g5e4d4a3
# User: docs
# Host: build-26886578-project-85767-nanograv-pint
# OS: Linux-5.19.0-1028-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-01-17T17:57:42.347245
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM3
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567793287
FINISH 56985.0000000465485648
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 2052.14475644944
CHI2R 1.0432866072442502
TRES 1.0204795556440334719
RAJ 5:00:00.00010106 1 0.00006205289580283582
DECJ 15:00:00.00813585 1 0.00725206049468805122
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.00000000000096198 1 3.1544984608191554236e-13
F1 -1.0001287029299695832e-15 1 1.3909954744010472733e-19
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 14.999996262588832235 1 4.793012518440144986e-06
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.0007298067796714427 1 0.0006935093127328683
TNREDAMP -12.602938664716227 0 0.103597755998806
TNREDGAM 2.4298946429074677 0 0.3945047262902932
TNREDC 13.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()
13 13
[10]:
<matplotlib.legend.Legend at 0x7fc6cda76c50>
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-01-17T17:59:27.293259
# PINT_version: 1.1.1+21.g5e4d4a3
# User: docs
# Host: build-26886578-project-85767-nanograv-pint
# OS: Linux-5.19.0-1028-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-01-17T17:58:30.972851
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM4
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567048495
FINISH 56985.0000000463896412
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1779.0658122872765
CHI2R 0.9203651382758802
TRES 0.944048694098943704
RAJ 4:59:59.99999947 1 0.00000184330994858866
DECJ 14:59:59.99997992 1 0.00016072158671647763
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.9999999999999425 1 3.6193347663840826427e-14
F1 -1.0000002099195742136e-15 1 8.2297317055894663914e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 14.999999073198278685 1 4.855380110854511593e-06
DMWXEPOCH 55000.0000000000000000
DMWXFREQ_0001 0.0002510040160586067
DMWXSIN_0001 -0.001151553512040076 1 5.984553858163869e-06
DMWXCOS_0001 0.0010370406363315419 1 6.498251777857674e-06
DMWXFREQ_0002 0.0005020080321172134
DMWXSIN_0002 -0.0005368488522170889 1 4.627614475000041e-06
DMWXCOS_0002 0.000328203115843672 1 4.436996250541984e-06
DMWXFREQ_0003 0.0007530120481758201
DMWXSIN_0003 0.00029553927179974844 1 4.421367064911405e-06
DMWXCOS_0003 -0.00027638270495988445 1 4.227863861855463e-06
DMWXFREQ_0004 0.0010040160642344267
DMWXSIN_0004 -0.00013626797328903268 1 4.313912400251054e-06
DMWXCOS_0004 3.4953313694412306e-05 1 4.223134857095127e-06
DMWXFREQ_0005 0.0012550200802930334
DMWXSIN_0005 -0.0001670140906768246 1 4.296050015477791e-06
DMWXCOS_0005 -0.00015628542325731048 1 4.2071183429061005e-06
DMWXFREQ_0006 0.0015060240963516401
DMWXSIN_0006 -3.2157331910736075e-05 1 4.336183366533151e-06
DMWXCOS_0006 4.606175278956103e-05 1 4.132758370122652e-06
DMWXFREQ_0007 0.0017570281124102468
DMWXSIN_0007 -0.0001110051969181526 1 4.289874349588694e-06
DMWXCOS_0007 3.5675759387288884e-05 1 4.160780473490053e-06
DMWXFREQ_0008 0.0020080321284688535
DMWXSIN_0008 3.3684084310658505e-06 1 4.348390638538037e-06
DMWXCOS_0008 7.37817016280541e-06 1 4.0904557771319364e-06
DMWXFREQ_0009 0.00225903614452746
DMWXSIN_0009 -5.6980364226098324e-05 1 4.245955010227331e-06
DMWXCOS_0009 -3.9095416811597414e-05 1 4.208000209347333e-06
DMWXFREQ_0010 0.002510040160586067
DMWXSIN_0010 2.8427711744856107e-05 1 4.281437997952018e-06
DMWXCOS_0010 -7.028897566228149e-05 1 4.221911652899444e-06
DMWXFREQ_0011 0.0027610441766446735
DMWXSIN_0011 -7.090317853727789e-06 1 6.784743351299488e-06
DMWXCOS_0011 2.0683297377927098e-05 1 6.7734939840610395e-06
DMWXFREQ_0012 0.0030120481927032802
DMWXSIN_0012 2.9250672330573596e-05 1 4.245613268175153e-06
DMWXCOS_0012 -2.5313794527061933e-05 1 4.223317662088848e-06
DMWXFREQ_0013 0.003263052208761887
DMWXSIN_0013 2.6432539007699978e-05 1 4.298275662594632e-06
DMWXCOS_0013 5.6950571033387366e-05 1 4.139406191412351e-06
DMWXFREQ_0014 0.0035140562248204936
DMWXSIN_0014 -3.7286350888642905e-05 1 4.227919600068393e-06
DMWXCOS_0014 3.0173320911888647e-05 1 4.217200778018193e-06
DMWXFREQ_0015 0.0037650602408791007
DMWXSIN_0015 2.404131365623997e-05 1 4.206899359784394e-06
DMWXCOS_0015 1.0424869288420435e-05 1 4.232778473055851e-06
DMWXFREQ_0016 0.004016064256937707
DMWXSIN_0016 -1.7767169726084476e-05 1 4.196340383065675e-06
DMWXCOS_0016 2.723535136308069e-05 1 4.246411347847873e-06
DMWXFREQ_0017 0.004267068272996314
DMWXSIN_0017 4.212607977389691e-06 1 4.22139258109507e-06
DMWXCOS_0017 -2.3998051675143747e-05 1 4.21058549951389e-06
DMWXFREQ_0018 0.00451807228905492
DMWXSIN_0018 1.5395487445695763e-05 1 4.210913419860878e-06
DMWXCOS_0018 -6.955653454581643e-07 1 4.221614764945713e-06
DMWXFREQ_0019 0.004769076305113527
DMWXSIN_0019 -1.908410864285359e-05 1 4.231395563085037e-06
DMWXCOS_0019 2.7939812781680224e-05 1 4.19250570767075e-06
DMWXFREQ_0020 0.005020080321172134
DMWXSIN_0020 5.130527099010848e-05 1 4.219052660648784e-06
DMWXCOS_0020 -5.595172298459322e-06 1 4.196510335108274e-06
DMWXFREQ_0021 0.00527108433723074
DMWXSIN_0021 -1.3310694451436146e-05 1 4.224231388928734e-06
DMWXCOS_0021 -1.5112959803444295e-06 1 4.193153587950733e-06
DMWXFREQ_0022 0.005522088353289347
DMWXSIN_0022 5.305231750869442e-06 1 4.202744828034573e-06
DMWXCOS_0022 -8.237522305511471e-06 1 4.235303719607682e-06
DMWXFREQ_0023 0.005773092369347954
DMWXSIN_0023 -1.2208475233207471e-05 1 4.189498093983882e-06
DMWXCOS_0023 1.3093127903110844e-06 1 4.240476487443357e-06
DMWXFREQ_0024 0.0060240963854065604
DMWXSIN_0024 -3.853828587847458e-06 1 4.14626350508526e-06
DMWXCOS_0024 1.1721481583756743e-05 1 4.273653609683575e-06
DMWXFREQ_0025 0.006275100401465167
DMWXSIN_0025 -2.197115108962593e-05 1 4.202887842061431e-06
DMWXCOS_0025 -3.5975484842496197e-06 1 4.2268368114268095e-06
DMWXFREQ_0026 0.006526104417523774
DMWXSIN_0026 2.7110349403827702e-05 1 4.192930319940569e-06
DMWXCOS_0026 -8.757877166692213e-06 1 4.245146219205668e-06
DMWXFREQ_0027 0.0067771084335823805
DMWXSIN_0027 8.826950107312473e-06 1 4.218310648965201e-06
DMWXCOS_0027 -8.391257895128306e-06 1 4.214124783293381e-06
DMWXFREQ_0028 0.007028112449640987
DMWXSIN_0028 -1.1446148087056823e-05 1 4.2631708977214035e-06
DMWXCOS_0028 3.7088792343966757e-06 1 4.153467949010788e-06
DMWXFREQ_0029 0.007279116465699594
DMWXSIN_0029 7.541324023145176e-06 1 4.2127380543240914e-06
DMWXCOS_0029 1.6517039512730127e-06 1 4.2043941882568174e-06
DMWXFREQ_0030 0.007530120481758201
DMWXSIN_0030 6.765289557926009e-06 1 4.116719045493218e-06
DMWXCOS_0030 6.068828088844911e-06 1 4.304746992821773e-06
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.00022189744439221265 1 5.522203999675227e-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-01-17T17:59:27.344469
# PINT_version: 1.1.1+21.g5e4d4a3
# User: docs
# Host: build-26886578-project-85767-nanograv-pint
# OS: Linux-5.19.0-1028-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-01-17T17:58:30.972851
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM4
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567048495
FINISH 56985.0000000463896412
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1779.0658122872765
CHI2R 0.9203651382758802
TRES 0.944048694098943704
RAJ 4:59:59.99999947 1 0.00000184330994858866
DECJ 14:59:59.99997992 1 0.00016072158671647763
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.9999999999999425 1 3.6193347663840826427e-14
F1 -1.0000002099195742136e-15 1 8.2297317055894663914e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 14.999999073198278685 1 4.855380110854511593e-06
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.00022189744439221265 1 5.522203999675227e-06
TNDMAMP -13.025694895440711 0 0.04160906899283973
TNDMGAM 2.8107811106980343 0 0.22595626896409493
TNDMC 30.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()
30 30
[18]:
<matplotlib.legend.Legend at 0x7fc6bff3a610>
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-01-17T17:59:37.548385
# PINT_version: 1.1.1+21.g5e4d4a3
# User: docs
# Host: build-26886578-project-85767-nanograv-pint
# OS: Linux-5.19.0-1028-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-01-17T17:59:27.950585
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM5
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567579977
FINISH 56985.0000000495309028
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 2041.2224614404722
CHI2R 1.0559867881223344
TRES 1.0123359507226174067
RAJ 5:00:00.00000014 1 0.00000144898542143598
DECJ 15:00:00.00000950 1 0.00012535786514063503
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.999999999999990334 1 2.8589319012325210774e-14
F1 -1.0000003608688937947e-15 1 6.385972048014054523e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.0
CM 1.2480692242163977131 1 0.05287599072231057712
TNCHROMIDX 3.5
CMWXEPOCH 55000.0000000000000000
CMWXFREQ_0001 0.00025100401605841234
CMWXSIN_0001 44.62823350958511 1 0.06796280156782948
CMWXCOS_0001 -138.67643940675558 1 0.07487885237207727
CMWXFREQ_0002 0.0005020080321168247
CMWXSIN_0002 54.56572631609782 1 0.060438845672905375
CMWXCOS_0002 15.88147191047747 1 0.06234203854992498
CMWXFREQ_0003 0.000753012048175237
CMWXSIN_0003 -31.582683516599662 1 0.06136394741159831
CMWXCOS_0003 -15.475161116148405 1 0.05986858770909579
CMWXFREQ_0004 0.0010040160642336494
CMWXSIN_0004 2.4629383326217518 1 0.06213768783544967
CMWXCOS_0004 3.8126346640983946 1 0.057943716842553235
CMWXFREQ_0005 0.0012550200802920618
CMWXSIN_0005 -7.820209227843221 1 0.06026726696891069
CMWXCOS_0005 -16.7056670815133 1 0.059681407737503724
CMWXFREQ_0006 0.001506024096350474
CMWXSIN_0006 3.7446227510811783 1 0.05871578597293184
CMWXCOS_0006 -2.8816019087956577 1 0.06093827592750248
CMWXFREQ_0007 0.0017570281124088863
CMWXSIN_0007 -1.689724360324224 1 0.05868434915819912
CMWXCOS_0007 7.020226199649258 1 0.06070351114137215
CMWXFREQ_0008 0.0020080321284672987
CMWXSIN_0008 2.6170351923547726 1 0.05907422119207655
CMWXCOS_0008 -0.37766326837779723 1 0.060299909843576305
CMWXFREQ_0009 0.002259036144525711
CMWXSIN_0009 1.3618778426298324 1 0.06018788247274132
CMWXCOS_0009 -5.0420345512688005 1 0.05906014467667199
CMWXFREQ_0010 0.0025100401605841235
CMWXSIN_0010 0.3657683034480025 1 0.06034339449778846
CMWXCOS_0010 1.3167333874704636 1 0.05931737619860385
CMWXFREQ_0011 0.0027610441766425355
CMWXSIN_0011 -0.9170447658141379 1 0.0729913687110301
CMWXCOS_0011 0.40871228532576015 1 0.07508415637150753
CMWXFREQ_0012 0.003012048192700948
CMWXSIN_0012 2.9837316409240104 1 0.06128694139880227
CMWXCOS_0012 -2.657092760910861 1 0.0577068358499603
CMWXFREQ_0013 0.0032630522087593603
CMWXSIN_0013 -1.3495629370456663 1 0.058152091521208135
CMWXCOS_0013 2.216784536940846 1 0.061082666489208653
CMWXFREQ_0014 0.0035140562248177727
CMWXSIN_0014 -0.09176225271652373 1 0.05922752607367493
CMWXCOS_0014 -1.2255036296325938 1 0.06007232228130556
CMWXFREQ_0015 0.003765060240876185
CMWXSIN_0015 0.6957344478074341 1 0.06019694584427198
CMWXCOS_0015 -0.176028234285564 1 0.059117488085845216
CMWXFREQ_0016 0.0040160642569345975
CMWXSIN_0016 0.3626352115680632 1 0.05946697265342072
CMWXCOS_0016 2.02161135084759 1 0.05968054267822264
CMWXFREQ_0017 0.00426706827299301
CMWXSIN_0017 -0.19190197382372332 1 0.05898986660147084
CMWXCOS_0017 -2.279551950728274 1 0.060121616346394804
CMWXFREQ_0018 0.004518072289051422
CMWXSIN_0018 -0.994286072721405 1 0.05987721623698582
CMWXCOS_0018 -1.6775734124462043 1 0.059214122979946826
CMWXFREQ_0019 0.004769076305109835
CMWXSIN_0019 0.07070992955244562 1 0.058253117251057675
CMWXCOS_0019 0.23547908595610811 1 0.06054359265096031
CMWXFREQ_0020 0.005020080321168247
CMWXSIN_0020 -0.07790966979400349 1 0.05929935994222159
CMWXCOS_0020 -0.8555879503514198 1 0.05946724481606512
CMWXFREQ_0021 0.0052710843372266595
CMWXSIN_0021 -1.5745184924925681 1 0.05878899127869617
CMWXCOS_0021 0.8246553881918322 1 0.059938079053846405
CMWXFREQ_0022 0.005522088353285071
CMWXSIN_0022 0.6867938090650405 1 0.05913270135436648
CMWXCOS_0022 0.4396447891401464 1 0.0597534368461192
CMWXFREQ_0023 0.005773092369343483
CMWXSIN_0023 0.6463484485128523 1 0.05943599304795097
CMWXCOS_0023 -0.5341093317860084 1 0.05925133627104786
CMWXFREQ_0024 0.006024096385401896
CMWXSIN_0024 -0.6885409190876414 1 0.05909387632554548
CMWXCOS_0024 0.6634216067829631 1 0.05951476706735295
CMWXFREQ_0025 0.006275100401460308
CMWXSIN_0025 -0.6277713070724573 1 0.058705988044645996
CMWXCOS_0025 -0.11563383546023186 1 0.05989650075295663
CMWXFREQ_0026 0.006526104417518721
CMWXSIN_0026 -0.15134934688167545 1 0.058770267770688826
CMWXCOS_0026 -0.13026078676571934 1 0.059591028816331076
CMWXFREQ_0027 0.006777108433577133
CMWXSIN_0027 0.15357443804922916 1 0.05923416694460236
CMWXCOS_0027 -0.309540401688459 1 0.058879476786857635
CMWXFREQ_0028 0.007028112449635545
CMWXSIN_0028 -0.5051889796663018 1 0.05912067502996916
CMWXCOS_0028 0.005114545212036668 1 0.05901335644804258
CMWXFREQ_0029 0.007279116465693958
CMWXSIN_0029 -0.00465894190621673 1 0.05829307936795083
CMWXCOS_0029 0.786257569596469 1 0.0596385939887182
CMWXFREQ_0030 0.00753012048175237
CMWXSIN_0030 0.11482387952556963 1 0.05922948761395996
CMWXCOS_0030 -0.36699660093327247 1 0.05852100057154917
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.0006239601003764885 1 4.226972032699588e-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-01-17T17:59:37.596711
# PINT_version: 1.1.1+21.g5e4d4a3
# User: docs
# Host: build-26886578-project-85767-nanograv-pint
# OS: Linux-5.19.0-1028-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-01-17T17:59:27.950585
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM5
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567579977
FINISH 56985.0000000495309028
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 2041.2224614404722
CHI2R 1.0559867881223344
TRES 1.0123359507226174067
RAJ 5:00:00.00000014 1 0.00000144898542143598
DECJ 15:00:00.00000950 1 0.00012535786514063503
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.999999999999990334 1 2.8589319012325210774e-14
F1 -1.0000003608688937947e-15 1 6.385972048014054523e-22
PEPOCH 55000.0000000000000000
PLANET_SHAPIRO N
DM 15.0
CM 1.2480692242163977131 1 0.05287599072231057712
TNCHROMIDX 3.5
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF -0.0006239601003764885 1 4.226972032699588e-06
TNCHROMAMP -12.994003726423763 0 0.04011109051640436
TNCHROMGAM 3.3885647401766037 0 0.22880023730203486
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 0x7fc6ccbfa810>
[ ]: