centroid: implement propagated uncertainties #938
Conversation
* follows same pattern as in line_flux of setting an uncertainty attribute on the returned quantity instead of returning an extra value * implements an private uncertainty._convert_uncertainty function to convert the Spectrum1D.uncertainty to any assumed type * hardcodes dispersion uncertainty to zero for now, but keeps in the equations in case we want to support passing/storing dispersion uncertainties in the future. Alternatively, the equations would simplify somewhat significantly if we drop support for uncertainties on dispersion entirely.
This comment was marked as resolved.
This comment was marked as resolved.
|
I'll also just add the same comment I did in #939 - I think at some point attaching the uncertainty to the quantity object should be revisited, as it also is easily dropped when converting units, but that's well outside the scope of these two PRs. |
There was a problem hiding this comment.
Mostly some minor suggestions, overall this looks generally good.
There is one comment about the "must be 1d or 2d" that is concerning, but that might be a pre-existing limitation that shouldn't be solved in this PR (I don't know what the current behavior is in >2D!).
On the broader topic of how to encode the uncertainties: I think the ship has sailed somewhat in that line_flux is following the .uncertainty -> stddev-style uncertainties. I think that's a fine starting point, as it is infinitely better than no uncertainties and a "better" way does not yet exist.
Going forward I think we may want to encode these instead either as
- StdDevUncertainty objects, with a deprecation period on the "old" way of just passing stddev style uncertainties without a wrapper object, or
- use
astropy.uncertainties. Unfrotunately that does not yet have a clean way to handle stddev/var/ivar-style uncertainties, though, so that's probably a project for the future.
Regardless, though, I think it would be good to add an uncertainty_type attribute which is set to std. That's how uncertainty objects encode it, and is sort of a soft duck-typing that might give us a fighting chance of following option 1 without a ton of pain. If you do that, though, we should also make that change in #939 and in line_flux for consistency.
| @@ -186,3 +187,26 @@ def _snr_derived(spectrum, region=None): | |||
| return signal / noise | |||
| else: | |||
| return 0.0 | |||
|
|
|||
|
|
|||
| def _convert_uncertainty(uncertainty, to_class): | |||
There was a problem hiding this comment.
This should really be in astropy.nddata. And hopefully will be in 5.1 with stropy/astropy#12057 ! We should change this to use that if the version is > 5.0.x, and fall back on this implementation (or a simple copy-paste of 12057 into here?) if not. That way we won't be forced to require Astropy 5.1 which isn't even out yet.
There was a problem hiding this comment.
@eteq I'm revisiting this now - are there good reasons to keep _convert_uncertainty rather than just bumping the astropy requirement to 5.1 in the next specutils release and using the represent_as machinery now available on the uncertainty objects?
| dispersion = np.tile(dispersion, [flux.shape[0], 1]) | ||
| elif len(flux.shape) > 2: | ||
| raise ValueError("spectrum must be 1D or 2D") |
There was a problem hiding this comment.
I'm a bit confused on this. Shouldn't it be OK for it to be any dimension, and the centroiding is just applied only along the spectral dimension?
There was a problem hiding this comment.
probably, but then we need an additional argument to control which index is the spectral dimension. That is probably out of scope here... but you could also argue that the new ValueError is out of scope, so I'm happy to remove that and/or open a new ticket to support any number of dimensions.
Co-authored-by: Ricky O'Steen <39831871+rosteen@users.noreply.github.com>
Co-authored-by: Erik Tollerud <erik.tollerud@gmail.com>
There was a problem hiding this comment.
As @kecnry asked, I independently derived this just to check the propagation math - if you want to check my work here it is:
I made the slightly different (simpler) assumption of no uncertainty on the dispersion, but it comes to the same result, so huzzah! All is good here. I did leave one inline suggestion for a minor performance optimization, though, which you can take or leave as you wish.
Two other large-scale thoughts:
- this derivation is only correct if the pixel uncertainties are uncorrelated. That's fine, because the correlated case is way more complicated anyway, and I don't think we need to say we'll support that right now. But we should just note in the docs somewhere that the uncertainty propogation assumes all pixel uncertainties are un-correlated.
- I agree we can leave the
uncertainty_typechange to a follow-on PR.
So I'm approving with the proviso that @kecnry looks at my suggestion and accepts it or not as you like, and at least consider putting something in the docstrings about the uncorrelated assuption.
Co-authored-by: Erik Tollerud <erik.tollerud@gmail.com>
|
Dev failure is also present on main (need to take a look at that...) so I'm merging. |
|
Leaving a note here since we've had some realizations that are getting implemented in #1057 : this derivation above is subtly incorrect: it hinged on an assumption about the order of summation that was incorrect. Attached is a whiteboard picture of a correct derivation (cross-checked with a few different other online sources).
|
line_fluxof setting an uncertainty attribute on the returned quantity if the inputSpectrum1Dhas assigned uncertainties.uncertainty._convert_uncertaintyfunction to convert theSpectrum1D.uncertaintyto any assumed type. This is now used inline_fluxas well.