Fitting a continuum model as a modelΒΆ

This example shows the initialization of a pyspeckit object from numpy arrays, as in Creating a Spectrum from scratch, but it adds the twist of including a steep continuum.

We fit the continuum using the polynomial continuum model, which gives access to the error on the polynomial fit parameters. No such parameters are accessible via the pyspeckit.Spectrum.baseline tools because they use numpy.poly1d to fit the data.

import numpy as np
import pyspeckit

xaxis = np.linspace(-50,150,100)
sigma = 10.
center = 50.

baseline = np.poly1d([0.1, 0.25])(xaxis)

synth_data = np.exp(-(xaxis-center)**2/(sigma**2 * 2.)) + baseline

# Add noise
stddev = 0.1
noise = np.random.randn(xaxis.size)*stddev
error = stddev*np.ones_like(synth_data)
data = noise+synth_data

# this will give a "blank header" warning, which is fine
sp = pyspeckit.Spectrum(data=data, error=error, xarr=xaxis,
                        xarrkwargs={'unit':'km/s'},
                        unit='erg/s/cm^2/AA')

sp.plotter()

sp.specfit.Registry.add_fitter('polycontinuum',
                               pyspeckit.models.polynomial_continuum.poly_fitter(),
                               2)

sp.specfit(fittype='polycontinuum', guesses=(0,0), exclude=[30, 70])

# subtract the model fit to create a new spectrum
sp_contsub = sp.copy()
sp_contsub.data -= sp.specfit.get_full_model()
sp_contsub.plotter()

# Fit with automatic guesses
sp_contsub.specfit(fittype='gaussian')

# Fit with input guesses
# The guesses initialize the fitter
# This approach uses the 0th, 1st, and 2nd moments
data = sp_contsub.data
amplitude_guess = data.max()
center_guess = (data*xaxis).sum()/data.sum()
width_guess = (data.sum() / amplitude_guess / (2*np.pi))**0.5
guesses = [amplitude_guess, center_guess, width_guess]
sp_contsub.specfit(fittype='gaussian', guesses=guesses)

sp_contsub.plotter(errstyle='fill')
sp_contsub.specfit.plot_fit()