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

PINT Noise Fitting Examples

[1]:
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 Fitter

import numpy as np
from io import StringIO
from astropy import units as u
from matplotlib import pyplot as plt
[2]:
setup_log(level="WARNING")
[2]:
1

Fitting for EFAC and EQUAD

[3]:
# Let us begin by simulating a dataset with an EFAC and an EQUAD.
# Note that the EFAC and the EQUAD are set as fit parameters ("1").
par = """
    PSR             TEST1
    RAJ             05:00:00    1
    DECJ            15:00:00    1
    PEPOCH          55000
    F0              100         1
    F1              -1e-15      1
    EFAC tel gbt    1.3         1
    EQUAD tel gbt   1.1         1
    TZRMJD          55000
    TZRFRQ          1400
    TZRSITE         gbt
    EPHEM           DE440
    CLOCK           TT(BIPM2019)
    UNITS           TDB
"""

m = get_model(StringIO(par))

ntoas = 200

# EFAC and EQUAD cannot be measured separately if all TOA uncertainties
# are the same. So we must set a different toa uncertainty for each TOA.
# This is how it is in real datasets anyway.
toaerrs = np.random.uniform(0.5, 2, ntoas) * u.us

t = make_fake_toas_uniform(
    startMJD=54000,
    endMJD=56000,
    ntoas=ntoas,
    model=m,
    obs="gbt",
    error=toaerrs,
    add_noise=True,
    include_bipm=True,
)
[4]:
# Now create the fitter. The `Fitter.auto()` function creates a
# Downhill fitter. Noise parameter fitting is only available in
# Downhill fitters.
ftr = Fitter.auto(t, m)
[5]:
# Now do the fitting.
ftr.fit_toas()
[6]:
# Print the post-fit model. We can see that the EFAC and EQUAD have been
# and the uncertainties are listed.
print(ftr.model)
# Created: 2024-06-05T07:29:12.476891
# 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                                 TEST1
EPHEM                               DE440
CLOCK                        TT(BIPM2019)
UNITS                                 TDB
START              53999.9999999852943056
FINISH             56000.0000000046727199
DILATEFREQ                              N
DMDATA                                  N
NTOA                                  200
CHI2                   199.99866689971557
CHI2R                  1.0362625227964537
TRES                2.1062198359451234314
RAJ                      4:59:59.99998855 1 0.00000761331760776661
DECJ                    15:00:00.00067522 1 0.00066414853789284706
PMRA                                  0.0
PMDEC                                 0.0
PX                                    0.0
F0                 100.000000000000044076 1 2.9970737003565711744e-13
F1              -1.0000142093174891752e-15 1 1.3576408974814889524e-20
PEPOCH             55000.0000000000000000
EFAC            tel gbt         1.447335388931792 1 0.15005801009405703
EQUAD           tel gbt        0.8770440884943274 1 0.21931025268215176
PLANET_SHAPIRO                          N
TZRMJD             55000.0000000000000000
TZRSITE                               gbt
TZRFRQ                             1400.0

[7]:
# Let us plot the injected and measured noise parameters together to
# compare them.
plt.scatter(m.EFAC1.value, m.EQUAD1.value, label="Injected", marker="o", color="blue")
plt.errorbar(
    ftr.model.EFAC1.value,
    ftr.model.EQUAD1.value,
    xerr=ftr.model.EFAC1.uncertainty_value,
    yerr=ftr.model.EQUAD1.uncertainty_value,
    marker="+",
    label="Measured",
    color="red",
)
plt.xlabel("EFAC_tel_gbt")
plt.ylabel("EQUAD_tel_gbt (us)")
plt.legend()
plt.show()
../_images/examples_noise-fitting-example_8_0.png

Fitting for ECORRs

[8]:
# Note the explicit offset (PHOFF) in the par file below.
# Implicit offset subtraction is typically not accurate enough when
# ECORR (or any other type of correlated noise) is present.
# i.e., PHOFF should be a free parameter when ECORRs are being fit.
par = """
    PSR             TEST2
    RAJ             05:00:00    1
    DECJ            15:00:00    1
    PEPOCH          55000
    F0              100         1
    F1              -1e-15      1
    PHOFF           0           1
    EFAC tel gbt    1.3         1
    ECORR tel gbt   1.1         1
    TZRMJD          55000
    TZRFRQ          1400
    TZRSITE         gbt
    EPHEM           DE440
    CLOCK           TT(BIPM2019)
    UNITS           TDB
"""

m = get_model(StringIO(par))

# ECORRs only apply when there are multiple TOAs per epoch.
# This can be simulated by providing multiple frequencies and
# setting the `multi_freqs_in_epoch` option. The `add_correlated_noise`
# option should also be set because correlated noise components
# are not simulated by default.

ntoas = 500
toaerrs = np.random.uniform(0.5, 2, ntoas) * u.us
freqs = np.linspace(1300, 1500, 4) * u.MHz

t = make_fake_toas_uniform(
    startMJD=54000,
    endMJD=56000,
    ntoas=ntoas,
    model=m,
    obs="gbt",
    error=toaerrs,
    freq=freqs,
    add_noise=True,
    add_correlated_noise=True,
    include_bipm=True,
    multi_freqs_in_epoch=True,
)
[9]:
ftr = Fitter.auto(t, m)
[10]:
ftr.fit_toas()
[10]:
True
[11]:
print(ftr.model)
# Created: 2024-06-05T07:29:31.829712
# 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                                 TEST2
EPHEM                               DE440
CLOCK                        TT(BIPM2019)
UNITS                                 TDB
START              53999.9999999862439353
FINISH             55984.0000000565202315
DILATEFREQ                              N
DMDATA                                  N
NTOA                                  500
CHI2                   500.02161856141134
CHI2R                  1.0163041027670963
TRES                1.7191506628907308634
RAJ                      4:59:59.99999766 1 0.00000582737861881591
DECJ                    15:00:00.00057374 1 0.00051081667296093815
PMRA                                  0.0
PMDEC                                 0.0
PX                                    0.0
F0                  100.00000000000024986 1 2.2680829620039313538e-13
F1              -9.999987965925370547e-16 1 1.0318240650939894726e-20
PEPOCH             55000.0000000000000000
EFAC            tel gbt        1.3465133040546062 1 0.048988698728936166
ECORR           tel gbt        1.0533720869724128 1 0.09793319690591715
TZRMJD             55000.0000000000000000
TZRSITE                               gbt
TZRFRQ                             1400.0
PHOFF               2.897066700660491e-05 1 1.7244508781765097e-05
PLANET_SHAPIRO                          N

[12]:
# Let us plot the injected and measured noise parameters together to
# compare them.
plt.scatter(m.EFAC1.value, m.ECORR1.value, label="Injected", marker="o", color="blue")
plt.errorbar(
    ftr.model.EFAC1.value,
    ftr.model.ECORR1.value,
    xerr=ftr.model.EFAC1.uncertainty_value,
    yerr=ftr.model.ECORR1.uncertainty_value,
    marker="+",
    label="Measured",
    color="red",
)
plt.xlabel("EFAC_tel_gbt")
plt.ylabel("ECORR_tel_gbt (us)")
plt.legend()
plt.show()
../_images/examples_noise-fitting-example_14_0.png