How specutils and muler propagate uncertainty

In this notebook we explore how specutils does uncertainty propagation, and how to combine two spectra with different—but overlapping—extents. That is, if a spectrum has close to, but not exactly the same wavelengths, how does specutils combine them? These two issues are important for combining spectra of the same objects taken from different exposures and different nights, when the target’s rest-frame wavelength solution can change greater than one pixel.

from specutils import Spectrum1D
import numpy as np
from astropy.nddata import StdDevUncertainty
import astropy.units as u
import matplotlib.pyplot as plt
%config InlineBackend.figure_format='retina'

Spectrum 1: \(S/N=20\)

First we’ll make a spectrum with signal-to-noise ratio equal to 20, and mean of 1.0.

N_points = 300
fake_wavelength = np.linspace(500, 600, num=N_points)*u.nm
mean_val, sigma = 1.0, 0.05
snr = mean_val / sigma
known_uncertainties = np.repeat(sigma, N_points) * u.Watt / u.cm**2
fake_flux = np.random.normal(loc=mean_val, scale=known_uncertainties) * u.Watt / u.cm**2

spec1 = Spectrum1D(spectral_axis=fake_wavelength,
                 flux=fake_flux,
                 uncertainty=StdDevUncertainty(known_uncertainties))

known_uncertainties.value[0:7]

array([0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05])

plt.axhline(1.0, linestyle='dashed', color='k', zorder=10)
plt.errorbar(spec1.wavelength.value, spec1.flux.value, yerr=spec1.uncertainty.array,
             linestyle='none', marker='o', ecolor='k', alpha=0.5)
plt.ylim(0, 1.5);

Spectrum 2: \(S/N = 50\) and conspicuously offset in wavelength

Now we’ll make a spectrum with signal-to-noise ratio equal to 50, and mean of 0.5. The wavelength axes are offset by 10 nanometers.

N_points2 = N_points
fake_wavelength2 = np.linspace(510, 610, num=N_points2)*u.nm
mean_val2, sigma2 = 0.5, 0.01
snr = mean_val2 / sigma2
known_uncertainties2 = np.repeat(sigma2, N_points2) * u.Watt / u.cm**2
fake_flux2 = np.random.normal(loc=mean_val2, scale=known_uncertainties2) * u.Watt / u.cm**2

spec2 = Spectrum1D(spectral_axis=fake_wavelength2,
                 flux=fake_flux2,
                 uncertainty=StdDevUncertainty(known_uncertainties2))

Add Spectrum 1 and Spectrum 2

What happens?

We expect the uncertainties to add in quadrature:

$ \sigma{net} = :nbsphinx-math:`sqrt{sigma_1^2 + sigma2^2}`$ $ :nbsphinx-math:`sigma`{net} = \sqrt{0.05^2 + 0.01^2} $

np.hypot(0.05, 0.01)

0.050990195135927854

spec_net = spec1 + spec2

spec_net.uncertainty[0:7]

StdDevUncertainty([0.0509902, 0.0509902, 0.0509902, 0.0509902, 0.0509902,
                   0.0509902, 0.0509902])

Woohoo! Specutils automatically propagates the error correctly! You can turn this error propagation off (I’m not sure why you would want to) by calling the method with a kwarg:

spec_net_no_error_propagation = spec1.add(spec2, propagate_uncertainties=False)

spec_net_no_error_propagation.uncertainty[0:7]

StdDevUncertainty([0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05])

Wait, but what about the offset? How did it deal with the non-overlapping edges?

plt.errorbar(spec1.wavelength.value, spec1.flux.value, yerr=spec1.uncertainty.array,
             linestyle='none', marker='o', label='Spec1: $S/N=20$')

plt.errorbar(spec2.wavelength.value, spec2.flux.value, yerr=spec2.uncertainty.array,
             linestyle='none', marker='o', markersize=1, label='Spec2: $S/N=50$')

plt.errorbar(spec_net.wavelength.value, spec_net.flux.value, yerr=spec_net.uncertainty.array,
             linestyle='none', marker='o',
             label='Spec net = Spec + Spec2 : $\sigma_{net}=\sqrt{\sigma^2 + \sigma_2^2}$')
plt.legend(loc='best')
plt.ylim(0, 2.5)

(0.0, 2.5)

Whoa! Specutils pretends like the signals are aligned to Spectrum 1. That is probably not the desired behavior for such an extreme offset as this one, but may be “good enough” for spectra that are either exactly aligned or within a pixel. It depends on your science application. PRV applications should not just round-to-the-nearest pixel, since they are trying to infer changes at sub-pixel levels.

What if you add the two spectra the other way around? Math says addition should be communitative…

spec_alt = spec2 + spec1

plt.errorbar(spec1.wavelength.value, spec1.flux.value, yerr=spec1.uncertainty.array,
             linestyle='none', marker='o', label='Spec1: $S/N=20$')

plt.errorbar(spec2.wavelength.value, spec2.flux.value, yerr=spec2.uncertainty.array,
             linestyle='none', marker='o', markersize=1, label='Spec2: $S/N=50$')

plt.errorbar(spec_alt.wavelength.value, spec_alt.flux.value, yerr=spec_alt.uncertainty.array,
             linestyle='none', marker='o', label='Spec alt = Spec2 + Spec : $\sigma_{net}=\sqrt{\sigma^2 + \sigma_2^2}$')

plt.legend(loc='best')
plt.ylim(0, 2.5);

Weird, so the wavelengths of the result are taken from the bounds of the first argument. This means “addition is not commutative” in specutils. Let’s see why:

spec_net.spectral_axis

$[500,~500.33445,~500.6689,~\dots,~599.3311,~599.66555,~600] \; \mathrm{nm}$

spec_alt.spectral_axis

$[510,~510.33445,~510.6689,~\dots,~609.3311,~609.66555,~610] \; \mathrm{nm}$

spec1.add(spec2, compare_wcs=None).spectral_axis

$[0,~1,~2,~\dots,~297,~298,~299] \; \mathrm{pix}$

Pixels instead of nanometers!

spec2.add(spec1, compare_wcs='first_found').spectral_axis

$[510,~510.33445,~510.6689,~\dots,~609.3311,~609.66555,~610] \; \mathrm{nm}$

spec1.add(spec2, compare_wcs='first_found').spectral_axis

$[500,~500.33445,~500.6689,~\dots,~599.3311,~599.66555,~600] \; \mathrm{nm}$

A ha! The compare_wcs kwarg controls what-to-do with the mis-matched spectral axes. Basically the sum of the two spectra just takes the wavelength labels from the first spectrum and uses those, when compare_wcs='first_found'—the default—is provided. It doesn’t actually interpolate or anything fancy…

And then adding them together.

Can we resample mis-aligned spectra and still get reasonable error propagation? How would that work?

from specutils.manipulation import FluxConservingResampler, LinearInterpolatedResampler

resampler = FluxConservingResampler(extrapolation_treatment='nan_fill')

%%capture
#This method throws a warning for an unknown reason...
resampled_spec2 = resampler(spec2, spec1.spectral_axis)

np.all(resampled_spec2.wavelength == spec1.wavelength)

True

resampled_spec2.uncertainty[0:50]

InverseVariance([           nan,            nan,            nan,
                            nan,            nan,            nan,
                            nan,            nan,            nan,
                            nan,            nan,            nan,
                            nan,            nan,            nan,
                            nan,            nan,            nan,
                            nan,            nan,            nan,
                            nan,            nan,            nan,
                            nan,            nan,            nan,
                            nan,            nan,            nan,
                 11043.15260749, 11043.15260748, 11043.15260749,
                 11043.15260749, 11043.15260749, 11043.15260749,
                 11043.15260749, 11043.15260748, 11043.15260749,
                 11043.15260749, 11043.15260748, 11043.15260749,
                 11043.15260748, 11043.15260749, 11043.15260749,
                 11043.15260748, 11043.15260749, 11043.15260748,
                 11043.15260748, 11043.15260749])

Hmmm… this process causes a new type of uncertainty “inverse variance” instead of std deviation… I’m not sure why! Hmm… We’ll manually convert? Variance is just standard-deviation squared. Inverse is just one-over-that…

new_sigma = np.sqrt(1/resampled_spec2.uncertainty.array)

resampled_spec2.uncertainty = StdDevUncertainty(new_sigma)

spec_final = spec1.add(resampled_spec2, propagate_uncertainties=True)

spec_final.uncertainty[0:50]

StdDevUncertainty([       nan,        nan,        nan,        nan,
                          nan,        nan,        nan,        nan,
                          nan,        nan,        nan,        nan,
                          nan,        nan,        nan,        nan,
                          nan,        nan,        nan,        nan,
                          nan,        nan,        nan,        nan,
                          nan,        nan,        nan,        nan,
                          nan,        nan, 0.05089748, 0.05089748,
                   0.05089748, 0.05089748, 0.05089748, 0.05089748,
                   0.05089748, 0.05089748, 0.05089748, 0.05089748,
                   0.05089748, 0.05089748, 0.05089748, 0.05089748,
                   0.05089748, 0.05089748, 0.05089748, 0.05089748,
                   0.05089748, 0.05089748])

Voila! It worked! I am not sure why specutils does not go the whole way and do this step manually.

plt.errorbar(spec1.wavelength.value, spec1.flux.value, yerr=spec1.uncertainty.array,
             linestyle='none', marker='o', label='Spec1: $S/N=20$')

plt.errorbar(spec2.wavelength.value, spec2.flux.value, yerr=spec2.uncertainty.array,
             linestyle='none', marker='o', markersize=1, label='Spec2: $S/N=50$')
plt.errorbar(resampled_spec2.wavelength.value, resampled_spec2.flux.value, yerr=new_sigma,
             linestyle='none', marker='o', markersize=1, label='Resampled Spec2: $S/N=50$')

plt.errorbar(spec_final.wavelength.value, spec_final.flux.value, yerr=spec_final.uncertainty.array,
             linestyle='none', marker='o', label='Spec final = ResampledSpec2 + Spec1 : $\sigma_{net}=?$')

plt.legend(loc='best')
plt.ylim(0, 2.5)

(0.0, 2.5)

Okay! Now we understand how specutils combines spectra, both for error propagation and for comparing the spectral axes of the two spectra.