astropy / specutils

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

Fit_continuum fitting method is not robust to noise #707

Open PatrickOgle opened 4 years ago

PatrickOgle commented 4 years ago

Using fitting.fit_continuum to fit a simple linear continuum to two continuum windows gives significantly different slopes for two different realizations of the noise. The chi^2 (by-eye) value sometimes appears to be unacceptable. Here are two examples for different realizations of the noise, in the following continuum windows: fitwindows1=[(5125.u.AA, 5175.u.AA),(5225.u.AA, 5275.u.AA)]

fit_continuum_2noise_realizations

PatrickOgle commented 4 years ago

fit_continuum calls fit_lines, which calls _fit_lines with fitter=fitting.LevMarLSQFitter (astropy.modeling.fitting.LevMarLSQFitter).

duytnguyendtn commented 4 years ago

The goal: get a "reasonable" fitting for a sharp emission line a good fraction of the time. (Better than 90%)

nmearl commented 4 years ago

Sounds like we should make a unit test to run a series of fit_continuums for different noise levels and ensure their result spread is acceptable.

alex180500 commented 3 years ago

I was trying to process a Black Body Radiation signal to delete the background noise but the Continuum method seems pretty off, am I doing something wrong?

Test

This is the code I implemented:

import matplotlib.pyplot as plt import numpy as np 
from astropy.modeling import models
from astropy import units as u
from specutils.spectra import Spectrum1D, SpectralRegion
from specutils.fitting import fit_generic_continuum

...

x = data[:, 0]
y = data[:, 1]
plt.plot(x, y, label = 'My Data', c='C0')

spectrum = Spectrum1D(flux=y*u.Jy, spectral_axis=x*u.um)
g1_fit = fit_generic_continuum(spectrum)
y_fit = g1_fit(x*u.um)

plt.plot(x, y_fit, label = ' Specutils Continuum Fit', c='C1')

plt.legend()
nmearl commented 3 years ago

Hi @alex180500. You'll need to exclude your line regions. You can do that by the following:

...
from specutils import SpectralRegion

spectrum = Spectrum1D(flux=y*u.Jy, spectral_axis=x*u.um)
g1_fit = fit_generic_continuum(spectrum, exclude_regions=[SpectralRegion(2 * u.um, 6 * u.um), SpectralRegion(58 * u.um, 65 * u.um)])
y_fit = g1_fit(x*u.um)

There are some examples in the documentation that you can look at as well.