Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Make whitened residuals available #1671

Merged
merged 7 commits into from
Nov 14, 2023
Merged
Show file tree
Hide file tree
Changes from 6 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CHANGELOG-unreleased.md
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@ the released changes.
- `ScaleToaError.register_toasigma_deriv_funcs` method to populate derivatives of scaled TOA uncertainties.
- `ScaleToaError.d_toasigma_d_EFAC` and `ScaleToaError.d_toasigma_d_EQUAD` methods.
- Separate `.fullname` for all observatories
- `Residuals.calc_whitened_resids()` method
### Fixed
- Wave model `validate()` can correctly use PEPOCH to assign WAVEEPOCH parameter
- Fixed RTD by specifying theme explicitly.
Expand Down
6 changes: 6 additions & 0 deletions src/pint/residuals.py
Original file line number Diff line number Diff line change
Expand Up @@ -500,6 +500,12 @@ def calc_time_resids(
)
return (phase_resids / self.get_PSR_freq(calctype=calctype)).to(u.s)

def calc_whitened_resids(self):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it better to have a separate method, or calc_time_resids() with a flag?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it should be a separate method since whitened residuals are dimensionless, whereas time resids have the dimentions of time.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Then probably add a "See Also" to the other residual methods

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Updated documentation.

r = self.calc_time_resids()
nr = sum(self.noise_resids.values())
sigma = self.get_data_error()
return ((r - nr) / sigma).to(u.dimensionless_unscaled)

def _calc_gls_chi2(self, lognorm=False):
"""Compute the chi2 when correlated noise is present in the timing model.
If the system is not singular, it uses Cholesky factorization to evaluate this.
Expand Down
45 changes: 44 additions & 1 deletion tests/test_residuals.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
from numpy.testing import assert_allclose
from pinttestdata import datadir

from pint.fitter import GLSFitter, WidebandTOAFitter, WLSFitter
from pint.fitter import Fitter, GLSFitter, WidebandTOAFitter, WLSFitter
from pint.models import get_model
from pint.residuals import CombinedResiduals, Residuals, WidebandTOAResiduals
import pint.residuals
Expand Down Expand Up @@ -438,3 +438,46 @@ def test_resid_mean_phase():
mean, _ = pint.residuals.weighted_mean(full, w)
assert mean == r.calc_phase_mean()
assert r.calc_phase_mean() == r2.calc_phase_mean()


@pytest.mark.parametrize(
dlakaplan marked this conversation as resolved.
Show resolved Hide resolved
"par",
[
"""
PSRJ J1234+5678
ELAT 0
ELONG 0
DM 10
F0 1
PEPOCH 58000
EFAC mjd 57000 58000 2
""",
"""
PSRJ J1234+5678
ELAT 0
ELONG 0
DM 10
F0 1
PEPOCH 58000
TNRedAmp -14.227505410948254
TNRedGam 4.91353
TNRedC 45
""",
],
)
def test_whitened_res(par):
m = get_model(StringIO(par))
t = make_fake_toas_uniform(
57000,
59000,
20,
model=m,
error=1 * u.us,
add_noise=True,
add_correlated_noise=m.has_correlated_errors,
)

ftr = Fitter.auto(t, m)
ftr.fit_toas()

assert np.isclose(ftr.resids.calc_whitened_resids().std(), 1, atol=0.75)
Loading