Model Fitting

General documentation for model fitting is below. There are several multi-component model fitting tools that have been developed based on pyspeckit:

• Mike Chen’s multi-vlsr

• Vlas Sokolov’s pyspecnest and multicube

• Jonny Henshaw’s SCOUSE

• Jared Keown’s astroclover

class pyspeckit.spectrum.fitters.Specfit(Spectrum, Registry=None)

Bases: Interactive

Declare interactive variables.

Must have a parent Spectrum class

Must declare button2action and button3action

EQW(plot=False, plotcolor='g', fitted=True, continuum=None, components=False, annotate=False, alpha=0.5, loc='lower left', xmin=None, xmax=None, xunits=None, continuum_as_baseline=False, xunit='pixel', midpt_location='plot-center')

Returns the equivalent width (integral of “baseline” or “continuum” minus the spectrum) over the selected range (the selected range defaults to self.xmin:self.xmax, so it may include multiple lines!)

Parameters

plot : bool

Plots a box indicating the EQW if plot==True (i.e., it will have a width equal to the equivalent width, and a height equal to the measured continuum)

fitted : bool

Use the fitted model? If false, uses the data

continuum : None or float

Can specify a fixed continuum with this keyword, otherwise will use the fitted baseline. WARNING: continuum=0 will still “work”, but will give numerically invalid results. Similarly, a negative continuum will work, but will yield results with questionable physical meaning.

continuum_as_baseline : bool

Replace the baseline with the specified continuum when computing the absorption depth of the line

components : bool

If your fit is multi-component, will attempt to acquire centroids for each component and print out individual EQWs

xmin : float

xmax : float

The range over which to compute the EQW

xunit : str

The units of xmin/xmax

midpt_location : ‘fitted’, ‘plot-center’

If ‘plot’ is set, this determines where the EQW will be drawn. It can be the fitted centroid or the plot-center, i.e. (xmin+xmax)/2

Returns

Equivalent Width, or widths if components=True :

add_sliders(parlimitdict=None, **kwargs)

Add a Sliders window in a new figure appropriately titled

Parameters

parlimitdict: dict :

Each parameter needs to have displayed limits; these are set in min-max pairs. If this is left empty, the widget will try to guess at reasonable limits, but the guessing is not very sophisticated yet.

.. todo:: Add a button in the navbar that makes this window pop up :

http://stackoverflow.com/questions/4740988/add-new-navigate-modes-in-matplotlib :

annotate(loc='upper right', labelspacing=0.25, markerscale=0.01, borderpad=0.1, handlelength=0.1, handletextpad=0.1, fontsize=10, frameon=False, chi2=None, optimal_chi2_kwargs={}, **kwargs)

Add a legend to the plot showing the fitted parameters

_clearlegend() will remove the legend

chi2 : {True or ‘reduced’ or ‘optimal’ or ‘allthree’}

kwargs passed to legend

button3action(event, debug=False, nwidths=1)

Disconnect the interactiveness Perform the fit (or die trying) Hide the guesses

clear(legend=True, components=True)

Remove the fitted model from the plot

Also removes the legend by default

clear_all_connections(debug=False)

Prevent overlapping interactive sessions

clear_highlights()

Hide and remove “highlight” colors from the plot indicating the selected region

copy(parent=None, registry=None)

Create a copy of the spectral fit - includes copies of the _full_model, the registry, the fitter, parinfo, modelpars, modelerrs, model, npeaks

Parameters

parent : pyspeckit.classes.Spectrum

A Spectrum instance that is the parent of the specfit instance. This needs to be specified at some point, but defaults to None to prevent overwriting a previous plot.

crop(x1pix, x2pix)

When spectrum.crop is called, this must be too

property dof

degrees of freedom in fit

downsample(factor)

Downsample the model spectrum (and the spectofit spectra) This should only be done when Spectrum.smooth is called

event_manager(event, force_over_toolbar=False, debug=False)

Decide what to do given input (click, keypress, etc.)

firstclick_guess()

Initialize self.guesses

property fitter

fullsizemodel()

If the model was fit to a sub-region of the spectrum, expand it (with zeros wherever the model was not defined) to fill the spectrum.

Examples

>>> noise = np.random.randn(100)
>>> xarr = np.linspace(-50,50,100)
>>> signal = np.exp(-(xarr-5)**2/(2*3.**2))
>>> sp = pyspeckit.Spectrum(data=noise + signal, xarr=xarr, xarrkwargs={'units':'km/s'})
>>> sp.specfit(xmin=-25,xmax=25)
>>> sp.specfit.model.shape
(48,)
>>> sp.specfit.fullsizemodel()
>>> sp.specfit.model.shape
(100,)

get_components(**kwargs)

If a model has been fitted, return the components of the model

Parameters

kwargs are passed to self.fitter.components :

get_emcee(nwalkers=None, **kwargs)

Get an emcee walker ensemble for the data & model using the current model type

Parameters

data : np.ndarray

error : np.ndarray

nwalkers : int

Number of walkers to use. Defaults to 2 * self.fitters.npars

Examples

>>> import pyspeckit
>>> x = pyspeckit.units.SpectroscopicAxis(np.linspace(-10,10,50), unit='km/s')
>>> e = np.random.randn(50)
>>> d = np.exp(-np.asarray(x)**2/2.)*5 + e
>>> sp = pyspeckit.Spectrum(data=d, xarr=x, error=np.ones(50)*e.std())
>>> sp.specfit(fittype='gaussian')
>>> emcee_ensemble = sp.specfit.get_emcee()
>>> p0 = emcee_ensemble.p0 * (np.random.randn(*emcee_ensemble.p0.shape) / 10. + 1.0)
>>> pos,logprob,state = emcee_ensemble.run_mcmc(p0,100)

get_full_model(debug=False, **kwargs)

compute the model over the full axis

get_model(xarr, pars=None, debug=False, add_baseline=None)

Compute the model over a given axis

get_model_frompars(xarr, pars, debug=False, add_baseline=None)

Compute the model over a given axis

get_model_xlimits(threshold='auto', peak_fraction=0.01, add_baseline=False, unit='pixels', units=None)

Return the x positions of the first and last points at which the model is above some threshold

Parameters

threshold : ‘auto’ or ‘error’ or float

If ‘auto’, the threshold will be set to peak_fraction * the peak model value. If ‘error’, uses the error spectrum as the threshold

peak_fraction : float

ignored unless threshold == ‘auto’

add_baseline : bool

Include the baseline when computing whether the model is above the threshold? default FALSE. Passed to get_full_model.

units : str

A valid unit type, e.g. ‘pixels’ or ‘angstroms’

get_pymc(**kwargs)

Create a pymc MCMC sampler from the current fitter. Defaults to ‘uninformative’ priors

kwargs are passed to the fitter’s get_pymc method, with parameters defined below.

Parameters

data : np.ndarray

error : np.ndarray

use_fitted_values : bool

Each parameter with a measured error will have a prior defined by the Normal distribution with sigma = par.error and mean = par.value

Examples

>>> x = pyspeckit.units.SpectroscopicAxis(np.linspace(-10,10,50), unit='km/s')
>>> e = np.random.randn(50)
>>> d = np.exp(-np.asarray(x)**2/2.)*5 + e
>>> sp = pyspeckit.Spectrum(data=d, xarr=x, error=np.ones(50)*e.std())
>>> sp.specfit(fittype='gaussian')
>>> MCuninformed = sp.specfit.get_pymc()
>>> MCwithpriors = sp.specfit.get_pymc(use_fitted_values=True)
>>> MCuninformed.sample(1000)
>>> MCuninformed.stats()['AMPLITUDE0']
>>> # WARNING: This will fail because width cannot be set <0, but it may randomly reach that...
>>> # How do you define a likelihood distribution with a lower limit?!
>>> MCwithpriors.sample(1000)
>>> MCwithpriors.stats()['AMPLITUDE0']

guesspeakwidth(event, debug=False, nwidths=1, **kwargs)

Interactively guess the peak height and width from user input

Width is assumed to be half-width-half-max

highlight_fitregion(drawstyle='steps-mid', color=(0, 0.8, 0, 0.5), linewidth=2, alpha=0.5, clear_highlights=True, **kwargs)

Re-highlight the fitted region

kwargs are passed to matplotlib.plot

history_fitpars()

integral(analytic=False, direct=False, threshold='auto', integration_limits=None, integration_limit_units='pixels', return_error=False, **kwargs)

Return the integral of the fitted spectrum

Parameters

analytic : bool

Return the analytic integral of the fitted function? .. WARNING:: This approach is only implemented for some models .. todo:: Implement error propagation for this approach

direct : bool

Return the integral of the spectrum (as opposed to the fit) over a range defined by the integration_limits if specified or threshold otherwise

threshold : ‘auto’ or ‘error’ or float

Determines what data to be included in the integral based off of where the model is greater than this number If ‘auto’, the threshold will be set to peak_fraction * the peak model value. If ‘error’, uses the error spectrum as the threshold See self.get_model_xlimits for details

integration_limits : None or 2-tuple

Manually specify the limits in integration_limit_units units

return_error : bool

Return the error on the integral if set. The error computed by sigma = sqrt(sum(sigma_i^2)) * dx

kwargs : :

passed to self.fitter.integral if not(direct)

Returns

np.scalar or np.ndarray with the integral or integral & error :

property mask

Mask: True means “exclude”

property mask_sliced

Sliced (subset) Mask: True means “exclude”

measure_approximate_fwhm(threshold='error', emission=True, interpolate_factor=1, plot=False, grow_threshold=2, **kwargs)

Measure the FWHM of a fitted line

This procedure is designed for multi-component blended lines; if the true FWHM is known (i.e., the line is well-represented by a single gauss/voigt/lorentz profile), use that instead. Do not use this for multiple independently peaked profiles.

This MUST be run AFTER a fit has been performed!

Parameters

threshold : ‘error’ | float

The threshold above which the spectrum will be interpreted as part of the line. This threshold is applied to the model. If it is ‘noise’, self.error will be used.

emission : bool

Is the line absorption or emission?

interpolate_factor : integer

Magnification factor for determining sub-pixel FWHM. If used, “zooms-in” by using linear interpolation within the line region

plot : bool

Overplot a line at the FWHM indicating the FWHM. kwargs are passed to matplotlib.plot

grow_threshold : int

Minimum number of valid points. If the total # of points above the threshold is <= to this number, it will be grown by 1 pixel on each side

Returns

The approximated FWHM, if it can be computed :

If there are <= 2 valid pixels, a fwhm cannot be computed :

model_mask(**kwargs)

Get a mask (boolean array) of the region where the fitted model is significant

Parameters

threshold : ‘auto’ or ‘error’ or float

The threshold to compare the model values to for selecting the mask region.

• auto: uses peak_fraction times the model peak

• error: use the spectrum error

• float: any floating point number as an absolute threshold

peak_fraction : float

Parameter used if threshold=='auto' to determine fraction of model peak to set threshold at

add_baseline : bool

Add the fitted baseline to the model before comparing to threshold?

Returns

mask : ndarray

A boolean mask array with the same size as the spectrum, set to True where the fitted model has values above a specified threshold

moments(fittype=None, **kwargs)

Return the moments

see the moments module

Parameters

fittype : None or str

The registered fit type to use for moment computation

multifit(fittype=None, renormalize='auto', annotate=None, show_components=None, verbose=True, color=None, guesses=None, parinfo=None, reset_fitspec=True, use_window_limits=None, use_lmfit=False, plot=True, **kwargs)

Fit multiple gaussians (or other profiles)

Parameters

fittype : str

What function will be fit? fittype must have been Registryed in the peakbgfitters dict. Uses default (‘gaussian’) if not specified

renormalize : ‘auto’ or bool

if ‘auto’ or True, will attempt to rescale small data (<1e-9) to be closer to 1 (scales by the median) so that the fit converges better

parinfo : parinfo structure

Guess structure; supercedes guesses

guesses : list or ‘moments’

Either a list of guesses matching the number of parameters * the number of peaks for the model, or ‘moments’ to fit a single spectrum with the moments as guesses

optimal_chi2(reduced=True, threshold='error', **kwargs)

Compute an “optimal” \(\chi^2\) statistic, i.e. one in which only pixels in which the model is statistically significant are included

Parameters

reduced : bool

Return the reduced \(\chi^2\)

threshold : ‘auto’ or ‘error’ or float

If ‘auto’, the threshold will be set to peak_fraction * the peak model value, where peak_fraction is a kwarg passed to get_model_xlimits reflecting the fraction of the model peak to consider significant If ‘error’, uses the error spectrum as the threshold

kwargs : dict

passed to get_model_xlimits()

Returns

chi2 : float

\(\chi^2\) statistic or reduced \(\chi^2\) statistic (\(\chi^2/n\))

\[\chi^2 = \sum( (d_i - m_i)^2 / e_i^2 )\]

peakbgfit(usemoments=True, annotate=None, vheight=True, height=0, negamp=None, fittype=None, renormalize='auto', color=None, use_lmfit=False, show_components=None, debug=False, use_window_limits=True, guesses=None, nsigcut_moments=None, plot=True, parinfo=None, **kwargs)

Fit a single peak (plus a background)

Parameters

usemoments : bool

The initial guess will be set by the fitter’s ‘moments’ function (this overrides ‘guesses’)

annotate : bool

Make a legend?

vheight : bool

Fit a (constant) background as well as a peak?

height : float

initial guess for background

negamp : bool

If True, assumes amplitude is negative. If False, assumes positive. If None, can be either.

fittype : bool

What function will be fit? fittype must have been Registryed in the peakbgfitters dict

renormalize : ‘auto’ or bool

if ‘auto’ or True, will attempt to rescale small data (<1e-9) to be closer to 1 (scales by the median) so that the fit converges better

nsigcut_moments : bool

pass to moment guesser; can do a sigma cut for moment guessing

plot_components(xarr=None, show_hyperfine_components=None, component_yoffset=0.0, component_lw=0.75, pars=None, component_fit_color='blue', component_kwargs={}, add_baseline=False, plotkwargs={}, **kwargs)

Overplot the individual components of a fit

Parameters

xarr : None

If none, will use the spectrum’s xarr. Otherwise, plot the specified xarr. This is useful if you want to plot a well-sampled model when the input spectrum is undersampled

show_hyperfine_components : None | bool

Keyword argument to pass to component codes; determines whether to return individual (e.g., hyperfine) components of a composite model

component_yoffset : float

Vertical (y-direction) offset to add to the components when plotting

component_lw : float

Line width of component lines

component_fitcolor : color

Color of component lines

component_kwargs : dict

Keyword arguments to pass to the fitter.components method

add_baseline : bool

Add the fit to the components before plotting. Makes sense to use if self.Spectrum.baseline.subtracted == False

pars : parinfo

A parinfo structure or list of model parameters. If none, uses best-fit

plot_fit(xarr=None, annotate=None, show_components=None, composite_fit_color='red', lw=0.5, composite_lw=0.75, pars=None, offset=None, use_window_limits=None, show_hyperfine_components=None, plotkwargs={}, **kwargs)

Plot the fit. Must have fitted something before calling this!

It will be automatically called whenever a spectrum is fit (assuming an axis for plotting exists)

kwargs are passed to the fitter’s components attribute

Parameters

xarr : None

If none, will use the spectrum’s xarr. Otherwise, plot the specified xarr. This is useful if you want to plot a well-sampled model when the input spectrum is undersampled

annotate : None or bool

Annotate the plot? If not specified, defaults to self.autoannotate

show_components : None or bool

show_hyperfine_components : None or bool

Show the individual gaussian components overlaid on the composite fit

use_window_limits : None or bool

If False, will reset the window to include the whole spectrum. If True, leaves the window as is. Defaults to self.use_window_limits if None.

pars : parinfo

A parinfo structure or list of model parameters. If none, uses best-fit

offset : None or float

Y-offset. If none, uses the default self.Spectrum.plotter offset, otherwise, uses the specified float.

plot_model(pars, offset=0.0, annotate=False, clear=False, **kwargs)

Plot a model from specified input parameters (see plot_fit for kwarg specification)

annotate is set to “false” because arbitrary annotations are not yet implemented

plotresiduals(fig=2, axis=None, clear=True, color='k', linewidth=0.5, drawstyle='steps-mid', yoffset=0.0, label=True, pars=None, zeroline=None, set_limits=True, **kwargs)

Plot residuals of the fit. Specify a figure or axis; defaults to figure(2).

Parameters

fig : int

Figure number. Overridden by axis

axis : axis

The axis to plot on

pars : None or parlist

If set, the residuals will be computed for the input parameters

zeroline : bool or None

Plot the “zero” line through the center of the residuals. If None, defaults to “True if yoffset!=0, False otherwise”

kwargs are passed to matplotlib plot :

print_fit(print_baseline=True, **kwargs)

Print the best-fit parameters to the command line

refit(use_lmfit=False)

Redo a fit using the current parinfo as input

register_fitter(*args, **kwargs)

Register a model fitter

Register a fitter function.

Parameters

name: string :

The fit function name.

function: function :

The fitter function. Single-fitters should take npars + 1 input parameters, where the +1 is for a 0th order baseline fit. They should accept an X-axis and data and standard fitting-function inputs (see, e.g., gaussfitter). Multi-fitters should take N * npars, but should also operate on X-axis and data arguments.

npars: int :

How many parameters does the function being fit accept?

savefit()

Save the fit parameters from a Gaussian fit to the FITS header

selectregion(xmin=None, xmax=None, xtype='wcs', highlight=False, fit_plotted_area=True, reset=False, verbose=False, debug=False, use_window_limits=None, exclude=None, **kwargs)

Pick a fitting region in either WCS units or pixel units

Parameters

*xmin / xmax* : [ float ]

The min/max X values to use in X-axis units (or pixel units if xtype is set). TAKES PRECEDENCE ALL OTHER BOOLEAN OPTIONS

*xtype* : [ string ]

A string specifying the xtype that xmin/xmax are specified in. It can be either ‘wcs’ or any valid xtype from pyspeckit.spectrum.units

*reset* : [ bool ]

Reset the selected region to the full spectrum? Only takes effect if xmin and xmax are not (both) specified. TAKES PRECEDENCE ALL SUBSEQUENT BOOLEAN OPTIONS

*fit_plotted_area* : [ bool ]

Use the plot limits as specified in :class:`pyspeckit.spectrum.plotters`? Note that this is not necessarily the same as the window plot limits!

*use_window_limits* : [ bool ]

Use the plot limits as displayed. Defaults to self.use_window_limits (pyspeckit.spectrum.interactive.use_window_limits). Overwrites xmin,xmax set by plotter

exclude: {list of length 2n,’interactive’, None} :

• interactive: start an interactive session to select the include/exclude regions

• list: parsed as a series of (startpoint, endpoint) in the spectrum’s X-axis units. Will exclude the regions between startpoint and endpoint

• None: No exclusion

seterrspec(usestd=None, useresiduals=True)

Simple wrapper function to set the error spectrum; will either use the input spectrum or determine the error using the RMS of the residuals, depending on whether the residuals exist.

setfitspec()

Set the spectrum that will be fit. This is primarily to remove NANs from consideration: if you simply remove the data from both the X-axis and the Y-axis, it will not be considered for the fit, and a linear X-axis is not needed for fitting.

However, it may be possible to do this using masked arrays instead of setting errors to be 1e10….

shift_pars(frame=None)

Shift the velocity / wavelength / frequency of the fitted parameters into a different frame

Right now this only takes care of redshift and only if redshift is defined. It should be extended to do other things later

start_interactive(debug=False, LoudDebug=False, reset_selection=False, print_message=True, clear_all_connections=True, **kwargs)

Initialize the interative session

Parameters

print_message : bool

Print the interactive help message?

clear_all_connections : bool

Clear all matplotlib event connections? (calls self.clear_all_connections())

reset_selection : bool

Reset the include mask to be empty, so that you’re setting up a fresh region.

property xmax

property xmin

• Fitting a user-defined model to a spectrum

[edit this page on github]

[edit this page on bitbucket]