This Jupyter notebook can be downloaded from rednoise-fit-example.ipynb, or viewed as a python script at rednoise-fit-example.py.

Red noise and DM 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 (
    dmwavex_setup,
    find_optimal_nharms,
    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 =  19
[5]:
print(np.argmin(d_aics))
19
[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")
../_images/examples_rednoise-fit-example_9_0.png
[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-06-05T07:30:41.191566
# PINT_version: 1.0+259.g224e5f1
# User: docs
# Host: build-24596653-project-85767-nanograv-pint
# OS: Linux-5.19.0-1028-aws-x86_64-with-glibc2.35
# Python: 3.11.6 (main, Feb  1 2024, 16:47:41) [GCC 11.4.0]
# Format: pint
PSR                                  SIM3
EPHEM                               DE440
CLOCK                        TT(BIPM2019)
UNITS                                 TDB
START              53000.9999999566359260
FINISH             56985.0000000463809144
DILATEFREQ                              N
DMDATA                                  N
NTOA                                 2000
CHI2                    1968.957720156762
CHI2R                  1.0071394988014128
TRES               0.99889189160602530386
RAJ                      5:00:00.00021733 1 0.00012190233068143402
DECJ                    15:00:00.00022464 1 0.01210114975198083552
PMRA                                  0.0
PMDEC                                 0.0
PX                                    0.0
F0                  100.00000000000134907 1 5.8476521001962098555e-13
F1              -9.999217218613239942e-16 1 1.8869604538516016345e-19
PEPOCH             55000.0000000000000000
PLANET_SHAPIRO                          N
DM                  14.999999984505273675 1 4.8551778324263834287e-06
WXEPOCH            55000.0000000000000000
WXFREQ_0001        0.00025100401605860257
WXSIN_0001         -2.874626780795822e-06 1 6.55028055751508e-07
WXCOS_0001         3.5804889889744396e-06 1 1.1396339300046183e-05
WXFREQ_0002         0.0005020080321172051
WXSIN_0002         1.0834571808736603e-06 1 3.305805190859114e-07
WXCOS_0002         1.1319679196187132e-06 1 2.893949191509805e-06
WXFREQ_0003         0.0007530120481758077
WXSIN_0003          7.535081674145011e-07 1 2.3022610033327065e-07
WXCOS_0003        -1.8003701768468572e-06 1 1.324801533491352e-06
WXFREQ_0004         0.0010040160642344103
WXSIN_0004        -2.8836173263730066e-07 1 1.839533495294969e-07
WXCOS_0004          5.650688245607681e-07 1 7.786498959872268e-07
WXFREQ_0005         0.0012550200802930128
WXSIN_0005        -3.2911784923097456e-07 1 1.59013993282829e-07
WXCOS_0005         -3.176474944691965e-07 1 5.326677636012073e-07
WXFREQ_0006         0.0015060240963516154
WXSIN_0006         -4.627028554571487e-07 1 1.4843701893007818e-07
WXCOS_0006         3.5657393989573524e-07 1 4.0553986877500255e-07
WXFREQ_0007          0.001757028112410218
WXSIN_0007          3.393172884232886e-07 1 1.4700334921391352e-07
WXCOS_0007        -1.4789241106682563e-07 1 3.3901139639592024e-07
WXFREQ_0008         0.0020080321284688205
WXSIN_0008         -5.055418163318304e-07 1 1.6095276174054566e-07
WXCOS_0008         1.4286618707469562e-08 1 3.1387516795370065e-07
WXFREQ_0009         0.0022590361445274233
WXSIN_0009          3.389026772160248e-07 1 2.0113382872699153e-07
WXCOS_0009        -1.8772896371178376e-08 1 3.3964325698596996e-07
WXFREQ_0010         0.0025100401605860257
WXSIN_0010         -6.563158445329482e-07 1 3.4873766411651243e-07
WXCOS_0010        -2.9450919928129823e-07 1 5.1700048418491e-07
WXFREQ_0011         0.0027610441766446284
WXSIN_0011         -7.352245968506125e-06 1 2.879703411844396e-06
WXCOS_0011        -1.4444366418212426e-06 1 3.7257594738506983e-06
WXFREQ_0012          0.003012048192703231
WXSIN_0012         4.4413042270818973e-07 1 2.113031415078226e-07
WXCOS_0012         1.3170236753397323e-07 1 2.387380735476677e-07
WXFREQ_0013         0.0032630522087618336
WXSIN_0013        -1.9821887183659818e-07 1 9.869162361214918e-08
WXCOS_0013         -4.983911441396678e-08 1 1.0080271450341067e-07
WXFREQ_0014          0.003514056224820436
WXSIN_0014           8.81337153877407e-08 1 6.381312696041468e-08
WXCOS_0014        -4.8706553309407484e-08 1 6.013778553839037e-08
WXFREQ_0015         0.0037650602408790387
WXSIN_0015         -9.156385726448888e-08 1 4.794522021439862e-08
WXCOS_0015        -5.0459163603222935e-08 1 4.4868218438665034e-08
WXFREQ_0016          0.004016064256937641
WXSIN_0016         1.1302749001396256e-07 1 4.1021738687671915e-08
WXCOS_0016         1.0101416841496889e-07 1 3.9248719107148235e-08
WXFREQ_0017          0.004267068272996243
WXSIN_0017         -6.016530032899183e-08 1 3.648259407628985e-08
WXCOS_0017           -6.9771265113416e-08 1 3.682754958251583e-08
WXFREQ_0018          0.004518072289054847
WXSIN_0018          6.089034594124779e-08 1 3.4894911851984295e-08
WXCOS_0018          6.534128557102885e-08 1 3.520564149044771e-08
WXFREQ_0019          0.004769076305113449
WXSIN_0019         -6.400675177485786e-10 1 3.3861806919037895e-08
WXCOS_0019          8.363813613462873e-08 1 3.468546001672525e-08
TZRMJD             55000.0000000000000000
TZRSITE                               gbt
TZRFRQ                             1400.0
PHOFF               0.0006974283187515753 1 0.000942421910155596

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)
19
[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-06-05T07:30:41.240345
# PINT_version: 1.0+259.g224e5f1
# User: docs
# Host: build-24596653-project-85767-nanograv-pint
# OS: Linux-5.19.0-1028-aws-x86_64-with-glibc2.35
# Python: 3.11.6 (main, Feb  1 2024, 16:47:41) [GCC 11.4.0]
# Format: pint
PSR                                  SIM3
EPHEM                               DE440
CLOCK                        TT(BIPM2019)
UNITS                                 TDB
START              53000.9999999566359260
FINISH             56985.0000000463809144
DILATEFREQ                              N
DMDATA                                  N
NTOA                                 2000
CHI2                    1968.957720156762
CHI2R                  1.0071394988014128
TRES               0.99889189160602530386
RAJ                      5:00:00.00021733 1 0.00012190233068143402
DECJ                    15:00:00.00022464 1 0.01210114975198083552
PMRA                                  0.0
PMDEC                                 0.0
PX                                    0.0
F0                  100.00000000000134907 1 5.8476521001962098555e-13
F1              -9.999217218613239942e-16 1 1.8869604538516016345e-19
PEPOCH             55000.0000000000000000
PLANET_SHAPIRO                          N
DM                  14.999999984505273675 1 4.8551778324263834287e-06
TZRMJD             55000.0000000000000000
TZRSITE                               gbt
TZRFRQ                             1400.0
PHOFF               0.0006974283187515753 1 0.000942421910155596
TNREDAMP               -12.84613974181327 0 0.09093304797501635
TNREDGAM                2.635906702321785 0 0.4467454952508081
TNREDC                               19.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()
19 19
[10]:
<matplotlib.legend.Legend at 0x7f781de39390>
../_images/examples_rednoise-fit-example_14_2.png

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 =  28
[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")
../_images/examples_rednoise-fit-example_20_0.png
[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-06-05T07:31:49.515920
# PINT_version: 1.0+259.g224e5f1
# User: docs
# Host: build-24596653-project-85767-nanograv-pint
# OS: Linux-5.19.0-1028-aws-x86_64-with-glibc2.35
# Python: 3.11.6 (main, Feb  1 2024, 16:47:41) [GCC 11.4.0]
# Format: pint
PSR                                  SIM4
EPHEM                               DE440
CLOCK                        TT(BIPM2019)
UNITS                                 TDB
START              53000.9999999566798032
FINISH             56985.0000000462515278
DILATEFREQ                              N
DMDATA                                  N
NTOA                                 2000
CHI2                   1940.9066974345815
CHI2R                  1.0020168804515135
TRES                 0.974971419443211205
RAJ                      4:59:59.99999934 1 0.00000184446431703325
DECJ                    14:59:59.99995008 1 0.00015889481056295898
PMRA                                  0.0
PMDEC                                 0.0
PX                                    0.0
F0                  99.999999999999991285 1 3.5340049636668685853e-14
F1              -9.999993527486125767e-16 1 8.187051942619115337e-22
PEPOCH             55000.0000000000000000
PLANET_SHAPIRO                          N
DM                  14.999993912193997757 1 4.874979983667827708e-06
DMWXEPOCH          55000.0000000000000000
DMWXFREQ_0001      0.00025100401605861357
DMWXSIN_0001         0.004849691131118741 1 5.8527737059155375e-06
DMWXCOS_0001        0.0013163045026811692 1 6.748662930486138e-06
DMWXFREQ_0002       0.0005020080321172271
DMWXSIN_0002       -0.0004067226040496782 1 4.710027071471525e-06
DMWXCOS_0002       -0.0004574283398274072 1 4.405535276264252e-06
DMWXFREQ_0003       0.0007530120481758408
DMWXSIN_0003       -0.0001737135873643711 1 4.3988699914506135e-06
DMWXCOS_0003      -0.00021746952883309445 1 4.344843284215329e-06
DMWXFREQ_0004       0.0010040160642344543
DMWXSIN_0004       -4.261579825964447e-06 1 4.359923597394595e-06
DMWXCOS_0004       0.00022378968671608732 1 4.288580333989539e-06
DMWXFREQ_0005        0.001255020080293068
DMWXSIN_0005      -0.00017829384510868818 1 4.305539903047698e-06
DMWXCOS_0005      -0.00023102385822171569 1 4.287235210945394e-06
DMWXFREQ_0006       0.0015060240963516815
DMWXSIN_0006       -5.479719100072904e-05 1 4.4509946906769015e-06
DMWXCOS_0006       -7.248845164412066e-05 1 4.130083231742793e-06
DMWXFREQ_0007        0.001757028112410295
DMWXSIN_0007       0.00022269328977934762 1 4.3148909438821635e-06
DMWXCOS_0007       -7.205073903378553e-05 1 4.276133475291716e-06
DMWXFREQ_0008       0.0020080321284689086
DMWXSIN_0008       -9.877208676864651e-06 1 4.3862447698879325e-06
DMWXCOS_0008       -4.868331351975415e-06 1 4.190251557942024e-06
DMWXFREQ_0009        0.002259036144527522
DMWXSIN_0009       -4.244657967126982e-05 1 4.247735904636898e-06
DMWXCOS_0009       1.5464844648060617e-05 1 4.330500739538359e-06
DMWXFREQ_0010        0.002510040160586136
DMWXSIN_0010        3.986177417192025e-05 1 4.278870504442912e-06
DMWXCOS_0010       2.2164698598897792e-06 1 4.370464810977569e-06
DMWXFREQ_0011       0.0027610441766447494
DMWXSIN_0011       2.4681798504108808e-06 1 6.8264109305655224e-06
DMWXCOS_0011       4.2082903570656184e-05 1 6.891512432717892e-06
DMWXFREQ_0012        0.003012048192703363
DMWXSIN_0012       2.5761016308428335e-05 1 4.381698257125023e-06
DMWXCOS_0012      -2.0838082690045576e-06 1 4.188488111148695e-06
DMWXFREQ_0013       0.0032630522087619767
DMWXSIN_0013       -8.056409655878467e-05 1 4.2554318321821665e-06
DMWXCOS_0013       -5.589741841825081e-05 1 4.287001345201776e-06
DMWXFREQ_0014         0.00351405622482059
DMWXSIN_0014      -2.6315860678508233e-06 1 4.22118734684745e-06
DMWXCOS_0014        8.015174593412657e-05 1 4.320216316061192e-06
DMWXFREQ_0015       0.0037650602408792035
DMWXSIN_0015       -4.178751395275333e-06 1 4.26368821026133e-06
DMWXCOS_0015        7.904882836920868e-06 1 4.277008049122148e-06
DMWXFREQ_0016        0.004016064256937817
DMWXSIN_0016       -2.891499682393343e-05 1 4.3550512811838695e-06
DMWXCOS_0016        9.801863157597333e-06 1 4.194290266080651e-06
DMWXFREQ_0017        0.004267068272996431
DMWXSIN_0017       1.2447695915581254e-05 1 4.2637441379013285e-06
DMWXCOS_0017        5.756160437585231e-06 1 4.274170599587617e-06
DMWXFREQ_0018        0.004518072289055044
DMWXSIN_0018       -5.287747254811066e-07 1 4.295055990214306e-06
DMWXCOS_0018       1.3630218905795319e-05 1 4.242354606332036e-06
DMWXFREQ_0019        0.004769076305113658
DMWXSIN_0019      -1.5384482036713192e-05 1 4.315961373827117e-06
DMWXCOS_0019      -1.1973908652031285e-05 1 4.227098472030087e-06
DMWXFREQ_0020        0.005020080321172272
DMWXSIN_0020        8.099035834374283e-06 1 4.2047228263026205e-06
DMWXCOS_0020        3.874574075930171e-06 1 4.332937005035777e-06
DMWXFREQ_0021        0.005271084337230885
DMWXSIN_0021       -1.599997680945171e-05 1 4.272621124418996e-06
DMWXCOS_0021      -2.3338857837056942e-05 1 4.27764072984849e-06
DMWXFREQ_0022        0.005522088353289499
DMWXSIN_0022       -3.007060839424949e-08 1 4.3280345953838585e-06
DMWXCOS_0022       -5.945844931490797e-06 1 4.2304225917092695e-06
DMWXFREQ_0023       0.0057730923693481125
DMWXSIN_0023      -4.8764646454461545e-06 1 4.308784824294432e-06
DMWXCOS_0023       -6.813361911844917e-06 1 4.2345743453827374e-06
DMWXFREQ_0024        0.006024096385406726
DMWXSIN_0024        6.839539049098091e-06 1 4.178676246899459e-06
DMWXCOS_0024      -1.6846093901909762e-06 1 4.360430020729325e-06
DMWXFREQ_0025         0.00627510040146534
DMWXSIN_0025      -3.1887443289411425e-06 1 4.236606162590421e-06
DMWXCOS_0025       1.0944873282459872e-05 1 4.295710041400588e-06
DMWXFREQ_0026        0.006526104417523953
DMWXSIN_0026       1.0512787087578663e-05 1 4.347149006387475e-06
DMWXCOS_0026          5.8125641264771e-07 1 4.190399613730798e-06
DMWXFREQ_0027        0.006777108433582567
DMWXSIN_0027         4.64734226768418e-06 1 4.3603455280402124e-06
DMWXCOS_0027        8.328038959391896e-06 1 4.182535739666894e-06
DMWXFREQ_0028         0.00702811244964118
DMWXSIN_0028         4.89870450572724e-06 1 4.35930886653655e-06
DMWXCOS_0028      -1.0751317373526608e-05 1 4.179602191773488e-06
TZRMJD             55000.0000000000000000
TZRSITE                               gbt
TZRFRQ                             1400.0
PHOFF              0.00012391400801415394 1 5.480560925276534e-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)
28
[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-06-05T07:31:49.574674
# PINT_version: 1.0+259.g224e5f1
# User: docs
# Host: build-24596653-project-85767-nanograv-pint
# OS: Linux-5.19.0-1028-aws-x86_64-with-glibc2.35
# Python: 3.11.6 (main, Feb  1 2024, 16:47:41) [GCC 11.4.0]
# Format: pint
PSR                                  SIM4
EPHEM                               DE440
CLOCK                        TT(BIPM2019)
UNITS                                 TDB
START              53000.9999999566798032
FINISH             56985.0000000462515278
DILATEFREQ                              N
DMDATA                                  N
NTOA                                 2000
CHI2                   1940.9066974345815
CHI2R                  1.0020168804515135
TRES                 0.974971419443211205
RAJ                      4:59:59.99999934 1 0.00000184446431703325
DECJ                    14:59:59.99995008 1 0.00015889481056295898
PMRA                                  0.0
PMDEC                                 0.0
PX                                    0.0
F0                  99.999999999999991285 1 3.5340049636668685853e-14
F1              -9.999993527486125767e-16 1 8.187051942619115337e-22
PEPOCH             55000.0000000000000000
PLANET_SHAPIRO                          N
DM                  14.999993912193997757 1 4.874979983667827708e-06
TZRMJD             55000.0000000000000000
TZRSITE                               gbt
TZRFRQ                             1400.0
PHOFF              0.00012391400801415394 1 5.480560925276534e-06
TNDMAMP               -13.009910798440412 0 0.043169779230571734
TNDMGAM                 3.611582591488531 0 0.26058268834491294
TNDMC                                28.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()
28 28
[18]:
<matplotlib.legend.Legend at 0x7f781c9b6c90>
../_images/examples_rednoise-fit-example_25_2.png