Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Inconsistent treatment of uncertainties in FluxConservingResampler? #1151

Open
kepplerm opened this issue Jul 6, 2024 · 0 comments
Open

Comments

@kepplerm
Copy link

kepplerm commented Jul 6, 2024

Hello,

I would like to raise a question regarding the implementation of the flux conserving resampling of spectra within specutils.

According to the documentation, the implementation follows the method described in this paper.

It seems to me that there is possibly some confusion of the uncertainty representation (standard deviation vs variance), and thus also of the error propagation in the specutils implementation.

My reasoning is the following: One input to the function is the uncertainty, called errs. According to the docstring, this is to be provided in representation of variance. In general, the variance var is defined as the square of the standard deviation std, hence var == std**2.

In the current implementation, error propagation is performed in the following way (see lines 212-217 of specutils/manipulation/resample.py):

            [...]
            if errs is not None:
                final_err = np.sum((errs[first_bin:final_bin+1] * p_ij) ** 2) / (sum_pij * sum_pij)
                output_errs[i] = np.sqrt(final_err)

    if errs is not None:
        output_errs = InverseVariance(np.reciprocal(output_errs))

Comparing this to the definition of the error propagation in the paper, this seems only correct if errs and output_errs are assumed to be represented as standard deviation, not variance (as stated in the docstring). See Equ. 4 of the paper linked above, or the corresponding screenshot below:

Screenshot 2024-07-05 at 5 03 49 PM

In my view, the corrected implementation should read in the following way for a consistent treatment of errs as variance, not standard deviation:

            [...]
            if errs is not None:
                final_err = np.sum(errs[first_bin:final_bin+1] * (p_ij** 2)) / (sum_pij * sum_pij)
                output_errs[i] = final_err

        if errs is not None:
            output_errs = InverseVariance(np.reciprocal(output_errs))

Apologies in case I got anything wrong! Does anyone have any feedback?

Thank you in advance!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant