search

CHIMEFRB / fitburst

An open-source package of utilities for direct modeling of radio dynamic spectra.
https://chimefrb.github.io/fitburst/
MIT License
11 stars 2 forks source link

design and add structure for computing post-fit statistics #12

Closed emmanuelfonseca closed 2 years ago

emmanuelfonseca commented 3 years ago

the current fitburst only computes an F-test based on separate fits where scattering is either directly modeled or ignored entirely. however, there are a number of statistical tests to be performed, whether they be for model comparison and/or quality of fits (e.g., tests of Gaussianty in best-fit residuals). we should design and implement code to perform these tests, probably as a separate class that receives several inputs from the best-fit model and computes a variety of statistical measures.

emmanuelfonseca commented 3 years ago

@bwmeyers no rush, but what do you think is the best way to start with adding stats-related functionality? should we add methods to the LSFitter? or should we instead create a separate class that contains methods that derive stats from fitter output?

if we want to add things to LSFitter, we may want to create a FitterBaseClass as the tests will presumably be fitter-independent. (and given the Catalog refs' comments, it seems like there are opportunities to implement lots of generic tests.)

bwmeyers commented 3 years ago

Generally, I would advocate for a separate class. The reason being that if at some later stage it is decided that fitburst should move away from a least-squares approach and into another domain (e.g., GLM) then we won't have to re-write all the functionality into a new fitter.

So, this could be achieved as you suggest by implementing a FitterBaseClass that actually defines all of the tests to run on the output and then leaves the actual fitting procedure up to its child. Or you could just create an independent class (maybe FitburstResults or something) that take output from a fitter and computes a variety of tests, presents parameters, etc. I think I like the latter idea more because that way you can interrogate the actual results without having to interact with the fitter object.

emmanuelfonseca commented 3 years ago

i agree and i quite like that concept @bwmeyers : we can focus on creating a base class that defines methods and things needed for all solvers, whether it be least-squares, MCMC, etc. i suggest we make this the focus of next week, as it's very much important due to Catalog 1 comments. (i should be done making DM fittable and adding in dispersion-smearing removal over the next few days.)

bwmeyers commented 3 years ago

Did we end up finalizing the path forward for this? There seems to have been significant progress in the other related road blocks like DM fitting, etc. We should probably set up a project board specifically for this issue alone so that we can make a list of the different statistics and tests we need to perform and track their implementation.

bwmeyers commented 3 years ago

For easy access, I think these are the relevant statistical editor comments (some of which may be much longer timescale projects than others):

a) The fitting procedure is least squares, which is a maximum likelihood estimator if the errors (scatter) is Gaussian. This assumption should be checked with nonparametric tests for normality (e.g. Shapiro-Wilk). Casual examination of Figs 6-7 suggests this might be reasonable for the frequency-integrated projections but not for the time-integrated projections. If Gaussianity is violated, then the LS approach is biased by outliers. Two mitigations are: generalized linear modeling (GLM) with a fat-tailed distribution for errors, or a robust regression procedure that downweights outliers (e.g. M-estimation with Huber's psi function).

(Certainly something worth investigating, although it will also be very easy to go down some deep rabbit holes. Probably not feasible to implement on the time scale of paper re-submission?)

b) The statistical model for the temporal shape of the burst (eqn 2) was derived my McKinnon (2014) under very specific mathematical assumptions (see the first sentence of his sec 2). While it does seem to fit well, perhaps there is no problem to solve. But from a statistical perspective, there are well-established statistical probability distributions with decades or centuries of mathematical development that resemble the asymmetrical burst profiles in McKinnon's Fig 1. These include the gamma distribution, beta distribution and their many extensions (Wikipedia, see e.g. the inverse gamma distribution). These have analytical mean, median (usually), variance, skewness, kurtosis, maximum likelihood estimator and other properties that could be tabulated in the catalog and may be useful for scientific analysis.

(I think it was decided that actually the functional form fit is physically motivated, so I suspect an explanation of that is sufficient to push back on making any of these changes, which are generically true but not really appropriate for our purposes.)

c) After a best-fit model is found with fitburst, then a goodness-of-fit measure should be applied. The simplest is an Anderson-Darling 1-sample test of the cumulative data and the model. If the fit is bad, then that might be scientifically interesting.

(This is probably the most easily tractable issue to implement ~quickly.)

emmanuelfonseca commented 2 years ago

an initial version of this is now complete; it can be refined as needed.