keflavich
commented
6 years ago 1.a) Yes, that is Andae's main point, but in practice they really show that the approximation does work in many cases, even though it is not rigorously justified.
1.b) Yes, lmfit does this under the hood. This: np.sqrt(np.diag(result.covar)) is not the default return. See https://github.com/lmfit/lmfit-py/blob/549c333c72c735698e20521d8714466a7ab74da8/lmfit/minimizer.py#L355, https://github.com/lmfit/lmfit-py/blob/549c333c72c735698e20521d8714466a7ab74da8/lmfit/minimizer.py#L1369, and astropy doing the same thing here: https://github.com/astropy/astropy/blob/7459062f222a15f9accb5e03a15cc2abb3580583/astropy/modeling/fitting.py#L770
I'm not inclined to include a deeper discussion of this topic in the pyspeckit paper, since we're just using tools that, under the hood, use an approximation to derive errors. I do want to alert users that the errors aren't rigorously defined, but this is a question far beyond this particular tool.
2) I like this idea, but it would be best as an appendix.
3) OK, good point.
4) Yep, need that.
5) We can hint at that, but not said the way you've said it. There are plenty of other assumptions that go into the MCMC analyses.
Some immediate fixes of typos and formatting are already in the PR. What follows are comments on the specific parts of the paper.
I am confused in this paragraph on the error estimation: "This covariance matrix is not directly the covariance of the parameters, and must be rescaled to deliver an approximate error. The standard rescaling is to multiply the covariance by the sum of the squared errors divided by the degree of freedom of the fit, usually referred to as χ2 /N." Where does the covariance matrix scaling come from? First of all, isn't Andrae et al. (2010) main point is that the degree of freedom of a nonlinear model is undefined in the first place, even if it is Taylor expansion-friendly? And then, looking at the way lmfit computes the errors, I don't see anything about the covariance matrix rescaling. If an example from https://lmfit.github.io/lmfit-py/model.html is run, then the uncertainties on the parameters turn out to be exactly the
np.sqrt(np.diag(result.covar)), with no apparent dependence on the reduced chi^2. Does it do it under the hood somewhere?A lot of new users coming from the radio side would be wondering how the molecular line emission profiles are computed. Should we add a section which starts off with a solution of the radiative transfer equation and goes on to explain how a hyperfine structure models are being generated (ideally with an example synthetic spectrum fit)? This echoes the "if this is the beefiest part of pyspeckit" comment closing Section 5.
I think this saying that
Cubes.fiteach"This tool is not the best tested component of pyspeckit" might scare people off! It can be read like there is a separate, not as reliable, way of fitting spectral cubes. We should instead say that we're using the same fitting procedure (minus the interactive UI bits) for the individual pixels on the cube as we do for normal spectra.The main development work on pyspeckit is on GitHub, but there's no GitHub link in the paper! Let's either add it as a footnote in the abstract or in a separate Contributing section?
After "These requirements are frequently not satisfied; see Andrae et al. (2010) and Andrae (2010) for details." - should we close the section by saying that for the error estimation with the least assumptions made one could use mcmc methods in pyspeckit?