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 = 18
[5]:
print(np.argmin(d_aics))
18
[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: 2024-12-18T17:11:27.452765
# PINT_version: 0+untagged.722.g91be315
# User: docs
# Host: build-26631508-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: 2024-12-18T17:10:25.441901
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM3
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999565982408
FINISH 56985.0000000463054861
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1971.5546966298277
CHI2R 1.0074372491721144
TRES 0.9792296514818439207
RAJ 4:59:59.99988424 1 0.00011666218796040997
DECJ 15:00:00.02819853 1 0.01157311722900301182
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.00000000000039768 1 5.64118466779309185e-13
F1 -1.0004049166875662421e-15 1 1.8035957488226566192e-19
PEPOCH 55000.0000000000000000
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.0008705617310160614 1 0.0009009690432405623
PLANET_SHAPIRO N
DM 15.000000286841134921 1 4.7736137901479509983e-06
WXEPOCH 55000.0000000000000000
WXFREQ_0001 0.00025100401605860533
WXSIN_0001 -5.10910643775862e-06 1 6.317046139305306e-07
WXCOS_0001 3.480261455961057e-05 1 1.089462711283908e-05
WXFREQ_0002 0.0005020080321172107
WXSIN_0002 2.0940955282251645e-06 1 3.1909185297362025e-07
WXCOS_0002 -6.273862707498201e-06 1 2.7682278339582713e-06
WXFREQ_0003 0.0007530120481758159
WXSIN_0003 -1.3317633326785057e-06 1 2.2283123681120997e-07
WXCOS_0003 2.7215260755531497e-06 1 1.2664328974886735e-06
WXFREQ_0004 0.0010040160642344213
WXSIN_0004 4.148871246247676e-07 1 1.7831662819372397e-07
WXCOS_0004 -2.456325567401831e-06 1 7.452199643158129e-07
WXFREQ_0005 0.0012550200802930267
WXSIN_0005 4.3958347705812804e-07 1 1.550716464646806e-07
WXCOS_0005 8.46362371968726e-07 1 5.09291229777495e-07
WXFREQ_0006 0.0015060240963516319
WXSIN_0006 -2.0556938565784705e-07 1 1.437277537076068e-07
WXCOS_0006 -7.32093704306219e-07 1 3.881572722625319e-07
WXFREQ_0007 0.0017570281124102373
WXSIN_0007 1.0421035066842556e-07 1 1.4314287329342078e-07
WXCOS_0007 5.126067374252659e-07 1 3.2580480054629137e-07
WXFREQ_0008 0.0020080321284688426
WXSIN_0008 -2.808455213782268e-07 1 1.558237229903988e-07
WXCOS_0008 -6.834376121747533e-07 1 3.008754551974722e-07
WXFREQ_0009 0.002259036144527448
WXSIN_0009 1.6214922997254517e-07 1 1.9400434829131743e-07
WXCOS_0009 8.11730425453154e-07 1 3.259404206589863e-07
WXFREQ_0010 0.0025100401605860534
WXSIN_0010 -1.1562394717105944e-07 1 3.3672269169703144e-07
WXCOS_0010 -1.3455595896739078e-06 1 4.951871357330058e-07
WXFREQ_0011 0.0027610441766446584
WXSIN_0011 -1.1817269701476921e-06 1 2.7757343206827652e-06
WXCOS_0011 -9.114818015519611e-06 1 3.5689725491322607e-06
WXFREQ_0012 0.0030120481927032637
WXSIN_0012 -3.1379951750813546e-08 1 2.026740593211293e-07
WXCOS_0012 4.804775391255473e-07 1 2.2806824238844415e-07
WXFREQ_0013 0.003263052208761869
WXSIN_0013 -5.219416459047305e-08 1 9.584696593463721e-08
WXCOS_0013 -3.174392494713679e-07 1 9.590091688520355e-08
WXFREQ_0014 0.0035140562248204745
WXSIN_0014 3.376400841883094e-08 1 6.050387962181686e-08
WXCOS_0014 2.3285622000327961e-07 1 5.75405541887561e-08
WXFREQ_0015 0.00376506024087908
WXSIN_0015 -1.6285610434333087e-08 1 4.6417502900735034e-08
WXCOS_0015 -4.6081144501191453e-08 1 4.3334240526488165e-08
WXFREQ_0016 0.004016064256937685
WXSIN_0016 4.962244136734745e-08 1 3.8914654450970734e-08
WXCOS_0016 2.0225854348266566e-09 1 3.8146327911366414e-08
WXFREQ_0017 0.00426706827299629
WXSIN_0017 -1.1387880382925195e-07 1 3.581328628461825e-08
WXCOS_0017 -8.196662853733965e-08 1 3.564533324234708e-08
WXFREQ_0018 0.004518072289054896
WXSIN_0018 5.64639383602681e-08 1 3.4087740829367463e-08
WXCOS_0018 -4.7146580863383764e-08 1 3.4366506173173226e-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)
18
[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: 2024-12-18T17:11:27.495721
# PINT_version: 0+untagged.722.g91be315
# User: docs
# Host: build-26631508-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: 2024-12-18T17:10:25.441901
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM3
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999565982408
FINISH 56985.0000000463054861
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 1971.5546966298277
CHI2R 1.0074372491721144
TRES 0.9792296514818439207
RAJ 4:59:59.99988424 1 0.00011666218796040997
DECJ 15:00:00.02819853 1 0.01157311722900301182
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.00000000000039768 1 5.64118466779309185e-13
F1 -1.0004049166875662421e-15 1 1.8035957488226566192e-19
PEPOCH 55000.0000000000000000
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.0008705617310160614 1 0.0009009690432405623
PLANET_SHAPIRO N
DM 15.000000286841134921 1 4.7736137901479509983e-06
TNREDAMP -12.710015250467169 0 0.09417242682653712
TNREDGAM 3.8392537550444636 0 0.42529882596538127
TNREDC 18.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()
18 18
[10]:
<matplotlib.legend.Legend at 0x7f694b0e01d0>
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 = 29
[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: 2024-12-18T17:12:38.871430
# PINT_version: 0+untagged.722.g91be315
# User: docs
# Host: build-26631508-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: 2024-12-18T17:11:28.156632
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM4
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567153472
FINISH 56985.0000000464802199
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 2051.370300978335
CHI2R 1.06013969042808
TRES 1.0054306033642953134
RAJ 4:59:59.99999893 1 0.00000186162579645673
DECJ 15:00:00.00018715 1 0.00015962181448812883
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.00000000000000528 1 3.6343774358014981518e-14
F1 -9.999999438059090705e-16 1 8.3090031639061892303e-22
PEPOCH 55000.0000000000000000
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.00025158241377740213 1 5.462397405967901e-06
PLANET_SHAPIRO N
DM 14.9999920268775959265 1 4.777980200906193114e-06
DMWXEPOCH 55000.0000000000000000
DMWXFREQ_0001 0.00025100401605860164
DMWXSIN_0001 0.001528911137361603 1 5.850758932184784e-06
DMWXCOS_0001 -0.000961363404530962 1 6.49842407630791e-06
DMWXFREQ_0002 0.0005020080321172033
DMWXSIN_0002 -0.0011388830597683656 1 4.512996081322161e-06
DMWXCOS_0002 0.0004808727370151203 1 4.375371376064127e-06
DMWXFREQ_0003 0.000753012048175805
DMWXSIN_0003 0.000749335495267624 1 4.423658793454535e-06
DMWXCOS_0003 0.0010346631323588098 1 4.030327115674731e-06
DMWXFREQ_0004 0.0010040160642344066
DMWXSIN_0004 9.341268845124963e-05 1 4.227862984616337e-06
DMWXCOS_0004 1.3294604018482131e-05 1 4.119564239045457e-06
DMWXFREQ_0005 0.0012550200802930083
DMWXSIN_0005 8.650881657523223e-05 1 4.202507365573531e-06
DMWXCOS_0005 0.0003927373471313164 1 4.138774783059289e-06
DMWXFREQ_0006 0.00150602409635161
DMWXSIN_0006 3.65269056591657e-07 1 4.133991066228027e-06
DMWXCOS_0006 0.00019483601473435508 1 4.1535607660836145e-06
DMWXFREQ_0007 0.0017570281124102117
DMWXSIN_0007 4.34047216368656e-05 1 4.147795848603989e-06
DMWXCOS_0007 1.1159461798053542e-05 1 4.104304118007916e-06
DMWXFREQ_0008 0.002008032128468813
DMWXSIN_0008 9.160282949963254e-06 1 4.156753135723396e-06
DMWXCOS_0008 -3.398693009945444e-05 1 4.12370197215394e-06
DMWXFREQ_0009 0.002259036144527415
DMWXSIN_0009 -6.56438424187592e-05 1 4.1533318640255745e-06
DMWXCOS_0009 -2.4911972675970918e-05 1 4.125716008302343e-06
DMWXFREQ_0010 0.0025100401605860165
DMWXSIN_0010 4.28659898641527e-06 1 4.152742400141663e-06
DMWXCOS_0010 -2.835433338888371e-05 1 4.175353448414525e-06
DMWXFREQ_0011 0.002761044176644618
DMWXSIN_0011 1.6991269617410833e-05 1 6.692632235885436e-06
DMWXCOS_0011 2.7305116607703203e-05 1 6.722435950061282e-06
DMWXFREQ_0012 0.00301204819270322
DMWXSIN_0012 9.164654799361362e-06 1 4.222368479930351e-06
DMWXCOS_0012 7.177741138489974e-05 1 4.0877214141254945e-06
DMWXFREQ_0013 0.0032630522087618214
DMWXSIN_0013 -8.010872518795092e-07 1 4.178202875577772e-06
DMWXCOS_0013 1.9406607570462305e-05 1 4.08248412013903e-06
DMWXFREQ_0014 0.0035140562248204233
DMWXSIN_0014 1.2021388612976657e-05 1 4.15020463327846e-06
DMWXCOS_0014 -6.328069152999817e-05 1 4.094274326070052e-06
DMWXFREQ_0015 0.003765060240879025
DMWXSIN_0015 3.697406059777383e-05 1 4.115348749504457e-06
DMWXCOS_0015 -7.180769433369052e-06 1 4.152076294671884e-06
DMWXFREQ_0016 0.004016064256937626
DMWXSIN_0016 7.994745762098073e-06 1 4.153470454096605e-06
DMWXCOS_0016 -2.8930892520455697e-05 1 4.12858788102982e-06
DMWXFREQ_0017 0.004267068272996228
DMWXSIN_0017 2.4722479880948792e-05 1 4.13957189110061e-06
DMWXCOS_0017 -6.820922145885122e-06 1 4.135238393831122e-06
DMWXFREQ_0018 0.00451807228905483
DMWXSIN_0018 5.036281405612548e-06 1 4.154823722372224e-06
DMWXCOS_0018 9.056386429914666e-06 1 4.1130949782391935e-06
DMWXFREQ_0019 0.004769076305113432
DMWXSIN_0019 7.931946668566553e-06 1 4.13259644976409e-06
DMWXCOS_0019 2.4531368495497435e-06 1 4.144849863579518e-06
DMWXFREQ_0020 0.005020080321172033
DMWXSIN_0020 -1.1986811613383546e-05 1 4.074342530344994e-06
DMWXCOS_0020 -2.1564958921618738e-05 1 4.202300518167373e-06
DMWXFREQ_0021 0.005271084337230635
DMWXSIN_0021 1.8003543840648036e-05 1 4.173014360213837e-06
DMWXCOS_0021 1.3954782325863442e-05 1 4.107851122249491e-06
DMWXFREQ_0022 0.005522088353289236
DMWXSIN_0022 9.802796691799707e-06 1 4.1505895552289615e-06
DMWXCOS_0022 2.622709426591976e-07 1 4.146537313931403e-06
DMWXFREQ_0023 0.005773092369347838
DMWXSIN_0023 3.329880430735563e-06 1 4.110433217175436e-06
DMWXCOS_0023 1.9328532842385752e-05 1 4.186277078219924e-06
DMWXFREQ_0024 0.00602409638540644
DMWXSIN_0024 -1.4536884067585144e-06 1 4.128826183610356e-06
DMWXCOS_0024 -3.855131588092268e-06 1 4.149060085415652e-06
DMWXFREQ_0025 0.006275100401465041
DMWXSIN_0025 -7.83539210383913e-06 1 4.165826979612949e-06
DMWXCOS_0025 -4.78004123328356e-06 1 4.115303671121349e-06
DMWXFREQ_0026 0.006526104417523643
DMWXSIN_0026 1.7972463943313383e-05 1 4.084017533210305e-06
DMWXCOS_0026 -4.302115850854455e-06 1 4.190570904239541e-06
DMWXFREQ_0027 0.006777108433582244
DMWXSIN_0027 1.950550344419311e-05 1 4.190753171771217e-06
DMWXCOS_0027 7.626228741384563e-07 1 4.0828374723332885e-06
DMWXFREQ_0028 0.007028112449640847
DMWXSIN_0028 -3.5979057652819476e-06 1 4.131865930930581e-06
DMWXCOS_0028 7.787377056277441e-06 1 4.157231459069434e-06
DMWXFREQ_0029 0.007279116465699448
DMWXSIN_0029 1.1168171041911957e-06 1 4.2221487720532615e-06
DMWXCOS_0029 1.6556235773208352e-05 1 4.071185007366203e-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)
29
[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: 2024-12-18T17:12:38.931013
# PINT_version: 0+untagged.722.g91be315
# User: docs
# Host: build-26631508-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: 2024-12-18T17:11:28.156632
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM4
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999567153472
FINISH 56985.0000000464802199
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 2051.370300978335
CHI2R 1.06013969042808
TRES 1.0054306033642953134
RAJ 4:59:59.99999893 1 0.00000186162579645673
DECJ 15:00:00.00018715 1 0.00015962181448812883
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 100.00000000000000528 1 3.6343774358014981518e-14
F1 -9.999999438059090705e-16 1 8.3090031639061892303e-22
PEPOCH 55000.0000000000000000
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.00025158241377740213 1 5.462397405967901e-06
PLANET_SHAPIRO N
DM 14.9999920268775959265 1 4.777980200906193114e-06
TNDMAMP -12.980545854867946 0 0.04334798376021426
TNDMGAM 3.6781680099523255 0 0.24252743414726424
TNDMC 29.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()
29 29
[18]:
<matplotlib.legend.Legend at 0x7f69484fd1d0>
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: 2024-12-18T17:12:51.802551
# PINT_version: 0+untagged.722.g91be315
# User: docs
# Host: build-26631508-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: 2024-12-18T17:12:39.441727
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM5
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566224190
FINISH 56985.0000000433934491
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 2030.3852452574922
CHI2R 1.0503803648512635
TRES 1.0090360839890963817
RAJ 4:59:59.99999853 1 0.00000145148823332992
DECJ 15:00:00.00020227 1 0.00012278158707374062
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999998131 1 2.7886952330336392031e-14
F1 -1.0000001951706878129e-15 1 6.308192678823079322e-22
PEPOCH 55000.0000000000000000
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.0005480946416878238 1 4.286351279443233e-06
PLANET_SHAPIRO N
DM 15.0
CM 1.1071966561231986351 1 0.051140809605846972163
TNCHROMIDX 3.5
CMWXEPOCH 55000.0000000000000000
CMWXFREQ_0001 0.0002510040160587901
CMWXSIN_0001 91.33060637006379 1 0.06676768203495426
CMWXCOS_0001 193.55806694966154 1 0.0710043362658943
CMWXFREQ_0002 0.0005020080321175802
CMWXSIN_0002 26.813793256433932 1 0.06014468071702132
CMWXCOS_0002 -0.7483703188456929 1 0.05926171406020733
CMWXFREQ_0003 0.0007530120481763702
CMWXSIN_0003 -7.819902805713786 1 0.05943946024629421
CMWXCOS_0003 -37.774781255686975 1 0.05783633564356317
CMWXFREQ_0004 0.0010040160642351603
CMWXSIN_0004 14.254686943790936 1 0.05872926490048164
CMWXCOS_0004 -18.411059321955502 1 0.057955348183790445
CMWXFREQ_0005 0.0012550200802939502
CMWXSIN_0005 -0.7526136759754258 1 0.05874396207282317
CMWXCOS_0005 -2.161422920872733 1 0.057743664904694424
CMWXFREQ_0006 0.0015060240963527404
CMWXSIN_0006 -0.37594648781296675 1 0.057517716291831396
CMWXCOS_0006 0.4373246824036515 1 0.05871162100102021
CMWXFREQ_0007 0.0017570281124115305
CMWXSIN_0007 -1.441152423702191 1 0.05665958659115779
CMWXCOS_0007 0.892845173867995 1 0.059579785860814305
CMWXFREQ_0008 0.0020080321284703206
CMWXSIN_0008 5.5662919353109155 1 0.058406388702938315
CMWXCOS_0008 -1.0920840667333254 1 0.057699377464156806
CMWXFREQ_0009 0.0022590361445291108
CMWXSIN_0009 1.685123619968438 1 0.05810162264079197
CMWXCOS_0009 5.647040806082003 1 0.05792992733521296
CMWXFREQ_0010 0.0025100401605879005
CMWXSIN_0010 3.3237709651464047 1 0.0566496918013604
CMWXCOS_0010 -2.8084425568884774 1 0.05964266671454535
CMWXFREQ_0011 0.0027610441766466906
CMWXSIN_0011 0.9967733081117067 1 0.07261351612212327
CMWXCOS_0011 0.07261786301483437 1 0.07129308867563579
CMWXFREQ_0012 0.0030120481927054807
CMWXSIN_0012 -1.8655978370630917 1 0.057884156651196866
CMWXCOS_0012 1.4547901567213195 1 0.05865007495752254
CMWXFREQ_0013 0.003263052208764271
CMWXSIN_0013 -0.39954132295971473 1 0.058731353973631335
CMWXCOS_0013 -1.0903478799271342 1 0.05785714351440871
CMWXFREQ_0014 0.003514056224823061
CMWXSIN_0014 0.9658059055056964 1 0.05876117308356738
CMWXCOS_0014 -0.6703729795320937 1 0.05782747580278207
CMWXFREQ_0015 0.003765060240881851
CMWXSIN_0015 0.6606342715753682 1 0.06032477026603267
CMWXCOS_0015 -0.2276431342286823 1 0.05633659881692205
CMWXFREQ_0016 0.004016064256940641
CMWXSIN_0016 -0.7741883028501193 1 0.05916353094070852
CMWXCOS_0016 0.33774177210233164 1 0.057133124592586945
CMWXFREQ_0017 0.004267068272999431
CMWXSIN_0017 -1.0512409594448395 1 0.05771155529214597
CMWXCOS_0017 -2.3009147674472126 1 0.05841674296352755
CMWXFREQ_0018 0.0045180722890582215
CMWXSIN_0018 0.15623635610078734 1 0.05816318182663128
CMWXCOS_0018 -0.1219316958364784 1 0.05814778470969154
CMWXFREQ_0019 0.004769076305117011
CMWXSIN_0019 0.2931771402797234 1 0.059132934428598065
CMWXCOS_0019 0.021588599843674002 1 0.056837541254797376
CMWXFREQ_0020 0.005020080321175801
CMWXSIN_0020 -0.22575185372936754 1 0.059189482035800225
CMWXCOS_0020 -0.044447206221996126 1 0.05693577459174301
CMWXFREQ_0021 0.0052710843372345915
CMWXSIN_0021 -0.49787281569213104 1 0.05817001397900512
CMWXCOS_0021 -0.3690320018773363 1 0.05792830239074001
CMWXFREQ_0022 0.005522088353293381
CMWXSIN_0022 0.5129874362262867 1 0.05852621425947568
CMWXCOS_0022 1.0420066009331985 1 0.0577788718404145
CMWXFREQ_0023 0.005773092369352172
CMWXSIN_0023 0.5281344887541117 1 0.0596434154710997
CMWXCOS_0023 -0.05757928966706748 1 0.05679774147108264
CMWXFREQ_0024 0.0060240963854109614
CMWXSIN_0024 0.23124519797233184 1 0.056348239353927895
CMWXCOS_0024 -1.117767931614879 1 0.05990008796654716
CMWXFREQ_0025 0.006275100401469751
CMWXSIN_0025 0.22054954812561098 1 0.05859039539976842
CMWXCOS_0025 0.02229956277140126 1 0.05777904647150758
CMWXFREQ_0026 0.006526104417528542
CMWXSIN_0026 -0.8563249448377719 1 0.05724048210240968
CMWXCOS_0026 -0.3721414745972516 1 0.05930214824249126
CMWXFREQ_0027 0.006777108433587331
CMWXSIN_0027 0.465010192060764 1 0.0589202438024504
CMWXCOS_0027 -0.4846832542990187 1 0.05756172240179164
CMWXFREQ_0028 0.007028112449646122
CMWXSIN_0028 0.07722222914797705 1 0.05902028686034095
CMWXCOS_0028 -0.1914614045443787 1 0.05729544570223176
CMWXFREQ_0029 0.007279116465704912
CMWXSIN_0029 -0.4419833848617347 1 0.0588702924366981
CMWXCOS_0029 0.10539791509070764 1 0.05749150044115044
CMWXFREQ_0030 0.007530120481763702
CMWXSIN_0030 0.21149402378473056 1 0.057330553457894284
CMWXCOS_0030 0.037284169700352865 1 0.059003545112163616
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: 2024-12-18T17:12:51.858972
# PINT_version: 0+untagged.722.g91be315
# User: docs
# Host: build-26631508-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: 2024-12-18T17:12:39.441727
# allow_tcb: False
# convert_tcb: False
# allow_T2: False
PSR SIM5
EPHEM DE440
CLOCK TT(BIPM2019)
UNITS TDB
START 53000.9999999566224190
FINISH 56985.0000000433934491
DILATEFREQ N
DMDATA N
NTOA 2000
CHI2 2030.3852452574922
CHI2R 1.0503803648512635
TRES 1.0090360839890963817
RAJ 4:59:59.99999853 1 0.00000145148823332992
DECJ 15:00:00.00020227 1 0.00012278158707374062
PMRA 0.0
PMDEC 0.0
PX 0.0
F0 99.99999999999998131 1 2.7886952330336392031e-14
F1 -1.0000001951706878129e-15 1 6.308192678823079322e-22
PEPOCH 55000.0000000000000000
TZRMJD 55000.0000000000000000
TZRSITE gbt
TZRFRQ 1400.0
PHOFF 0.0005480946416878238 1 4.286351279443233e-06
PLANET_SHAPIRO N
DM 15.0
CM 1.1071966561231986351 1 0.051140809605846972163
TNCHROMIDX 3.5
TNCHROMAMP -13.053338001368406 0 0.040178877918090265
TNCHROMGAM 3.57245027097845 0 0.21461188725779567
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 0x7f694b23bb50>
[ ]: