astropy / specutils

An Astropy coordinated package for astronomical spectroscopy. Maintainers: @rosteen @keflavich @eteq
http://specutils.readthedocs.io/en/latest/
166 stars 126 forks source link

Have a reasonable uncertainty treatment in `line_flux` #627

Closed eteq closed 3 months ago

eteq commented 4 years ago

The line_flux function in specutils.analysis should yield a measurement that includes uncertainties when a spectrum has one of the uncertainty types that is consistent with "gaussian" uncertainties (i.e. stddev, variance, or inverse variance).  The implementation is pretty straightforward: the pixels in the flux measurement should be weighted by their inverse-variance to yield the flux estimate.

Getting an uncertainty of the measurement is a bit more complicated, but is still well-defined following standard error-propogation rules for gaussian uncertainties. 

A question I don't have a clear answer to is: how best to pass on the uncertainty in the derived value. I'll post a follow-on comment with my current idea.

eteq commented 4 years ago

To address the "how to out the uncertainty" question: We don't want to change the API to have a second return value because that's awkward in the case of no uncertainties. So it would be better if the return value carried its own uncertainty.

However we currently don't have a way of doing that consistently in Astropy. As a workaround, perhaps the following:

>>> flux = line_flux(spectrum, ...)
>>> fluxunc = flux.uncertainty
>>> assert isinstance(fluxunc, Quantity)

that is actually pretty straightofrward to "hack" inside line_flux because Quantity objects can get arbitrary attributes assigned. E.g.:

flux = [1,2,3]*u.Jy
flux.uncertainty = StdDevUncertainty([.1,.15, .18]*u.Jy)

would work right now. I'm inclined to suggest just doing that and later on adapting the API when there's something in astropy itself to support this.

An alternative is to use the astropy.uncertainties machinery (which conveniently I am a maintainer of), but it currently has some performance penalties we may not be comfortable with. I think we can get it to behave a lot like my above suggestion if we want, though, so there's a potential upgrade path.

hcferguson commented 4 years ago

This seems like a good solution to me. Although I'm confused how flux.uncertainty can simultaneously be an instance of Quantity and a StdDevUncertainty object. Can it?

dhomeier commented 4 years ago

This seems like a good solution to me. Although I'm confused how flux.uncertainty can simultaneously be an instance of Quantity and a StdDevUncertainty object. Can it?

Does not seem to preserve the original property in current Astropy at least:

>>> flux = [1, 2, 3] * u.Jy
>>> isinstance(flux, u.Quantity)
True
>>> flux.uncertainty = StdDevUncertainty([.1, .15, .18] * u.Jy)
>>> isinstance(flux.uncertainty, u.Quantity)
False
>>> isinstance(flux.uncertainty, StdDevUncertainty)
True
nmearl commented 4 years ago

Does not seem to preserve the original property in current Astropy at least:

StdDevUncertainty (and NDUncertanty) in general are not subclasses of Quantity.

rosteen commented 3 months ago

Uncertainty in line flux was implemented in https://github.com/astropy/specutils/pull/669, I imagine this should have been closed.