Fit white noise parameters

[abhisrkckl]

WLSNoiseFitter can do maximum-likelihood estimates of white noise noise parameters (EFACs and EQUADs) along with their uncertainties.

This is done by maximizing the normalized log-likelihood function over the timing noise parameters and the noise parameters alternately. The timing model maximization is done using DownhillWLSFitter. The noise model maximization is done using the Newton-CG algorithm available in scipy. This requires derivatives wrt noise parameters, which are also implemented here. The noise parameter uncertainties can be computed using the Hessian (obtained using numdifftools).

The drawback of this approach is that we can't estimate the covariances between timing model and noise parameters.

Edit: I am also including Residuals.lnlikelihood() and utils.akaike_information_criterion in this PR since the underlying code for those functions are already present here.

[codecov]

Codecov Report

Attention: 11 lines in your changes are missing coverage. Please review.

Comparison is base (41868ad) 68.40% compared to head (9b4739e) 68.50%.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1639      +/-   ##
==========================================
+ Coverage   68.40%   68.50%   +0.10%     
==========================================
  Files         104      104              
  Lines       24151    24264     +113     
  Branches     4308     4329      +21     
==========================================
+ Hits        16520    16622     +102     
- Misses       6548     6555       +7     
- Partials     1083     1087       +4     

Files	Coverage Δ
src/pint/fitter.py	84.08% <100.00%> (+0.45%)	⬆️
src/pint/models/parameter.py	82.44% <ø> (ø)
src/pint/simulation.py	85.62% <100.00%> (+0.17%)	⬆️
src/pint/residuals.py	76.20% <96.42%> (+1.39%)	⬆️
src/pint/utils.py	57.08% <71.42%> (+0.11%)	⬆️
src/pint/models/noise_model.py	93.49% <89.18%> (-0.55%)	⬇️
src/pint/models/timing_model.py	83.71% <71.42%> (-0.11%)	⬇️

... and 1 file with indirect coverage changes

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[abhisrkckl]

Example:

from pint.models import get_model
from pint.simulation import make_fake_toas_uniform
from pint.noise_fitter import WLSNoiseFitter
from io import StringIO

par = """
    ELAT    1.3     1
    ELONG   2.5     1
    F0      100     1
    F1      1e-13   1
    PEPOCH  55000
    EFAC mjd 50000 53000 2      1
    EQUAD mjd 53000 55000 1    1
"""
m = get_model(StringIO(par))
t = make_fake_toas_uniform(50000, 55000, 1000, m, add_noise=True)

ftr = WLSNoiseFitter(t, m, uncertainty_method="hess")
ftr.fit_toas()

Pre-fit par file:

DILATEFREQ                              N
DMDATA                                  N
NTOA                                    0
CHI2                                  0.0
ELONG                   2.500000000000000 1 0.00000000000000000000
ELAT                    1.300000000000000 1 0.00000000000000000000
PMELONG                               0.0
PMELAT                                0.0
PX                                    0.0
ECL                              IERS2010
F0                                  100.0 1 0.0
F1                                  1e-13 1 0.0
PEPOCH             55000.0000000000000000
PLANET_SHAPIRO                          N
EFAC            mjd 50000.0 53000.0                       2.0 1 0.0
EQUAD           mjd 53000.0 55000.0                       1.0 1 0.0
TZRMJD             55000.0000000000000000
TZRSITE                               ssb
TZRFRQ                                inf

Post-fit par file:

EPHEM                               DE421
CLOCK                             TT(TAI)
START              49999.9999999475177662
FINISH             54999.9999999909061806
DILATEFREQ                              N
DMDATA                                  N
NTOA                                 1000
CHI2                                  0.0
ELONG                   2.500000006598401 1 0.00000000864327838779
ELAT                    1.299999534194669 1 0.00000038051831122624
PMELONG                               0.0
PMELAT                                0.0
PX                                    0.0
ECL                              IERS2010
F0                  99.999999999999908185 1 1.66510871960953868e-13
F1                9.99999995394515813e-14 1 7.7723152376714238183e-22
PEPOCH             55000.0000000000000000
PLANET_SHAPIRO                          N
EFAC            mjd 50000.0 53000.0        2.0050491397917654 1 0.05798063163434304
EQUAD           mjd 53000.0 55000.0        0.9787394583552456 1 0.07081435918133107
TZRMJD             55000.0000000000000000
TZRSITE                               ssb
TZRFRQ                                inf

[dlakaplan]

Prospects for Wideband data?

[dlakaplan]

and how does this work with correlated errors?

[dlakaplan]

Rather than defining yet another fitter type, can you make the DownhillWLS fitter just use the appropriate code when free noise parameters are included?

[abhisrkckl]

Prospects for Wideband data?

This method should in principle work for wideband also. But I have not gotten that far yet. I think I'll do the wideband version as a separate PR.

[abhisrkckl]

and how does this work with correlated errors?

Should work. It should be a separate PR I think.

[abhisrkckl]

Rather than defining yet another fitter type, can you make the DownhillWLS fitter just use the appropriate code when free noise parameters are included?

I am testing things out at the moment, and don't want to mess with existing code. I'll move the functionality to DownhillWLSFitter once I am confident everything is in order.

[abhisrkckl]

Here is a comparison between frequentist estimate using PINT (blue curve is the Gaussian distribution corresponding to the measured value and uncertainty) and Bayesian estimate using ENTERPRISE (black histogram).

[dlakaplan]

Looks good. There is a slight bias in MJD2_EFAC: is that just random or is it always like that?

[dlakaplan]

It's also worth discussing whether adding numdifftools as a requirement is needed. Or if it can just be recommended (i.e., import within the function, error if it can't)

[abhisrkckl]

Looks good. There is a slight bias in MJD2_EFAC: is that just random or is it always like that?

Here is the result of another run with the same setup:

It looks like the posterior is more like a chi2 distribution with a slightly longer tail to the right. If that's the case, the Gaussian will be slightly biased to lower values because its mean is at the maximum likelihood value.

I have run this multiple times with different noise realizations and it looks like this slight bias is consistent.

[abhisrkckl]

Here is a comparison of injected value vs measured value for different realizations of the same par file. It seems that the values are correctly recovered within 2-sigma level.

[dlakaplan]

Has anybody given feedback on slack about the new requirement? If not, we can probably just include it and I can go ahead and merge this.

[abhisrkckl]

Has anybody given feedback on slack about the new requirement?

No.

[abhisrkckl]

Please merge this if you think this is OK.

[dlakaplan]

Is there a notebook or code snipped in the docstring describing how to use this? How to retrieve the full likelihood (not just chi^2) and compute AIC?

[abhisrkckl]

Is there a notebook or code snipped in the docstring describing how to use this? How to retrieve the full likelihood (not just chi^2) and compute AIC?

The AIC function is simply used as pint.utils.akaike_information_criterion(model, toas). It has a docstring.

Similarly, the likelihood function is residuals.lnlikelihood() similar to residuals.calc_chi2(). I think this is pretty self-explanatory.

[dlakaplan]

But if somebody is coming into this and wants to understand the full cycle of fit/evaluate, having a short worked example helps.

[abhisrkckl]

OK... I'll add an example.