Model Fitting¶

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

class pyspeckit.spectrum.fitters.Specfit(Spectrum, Registry=None)[source] [github] [bitbucket]

Bases: pyspeckit.spectrum.interactive.Interactive

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')[source] [github] [bitbucket]

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 Equivalent Width, or widths if components=True :
add_sliders(parlimitdict=None, **kwargs)[source] [github] [bitbucket]

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)[source] [github] [bitbucket]

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)[source] [github] [bitbucket]

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

clear(legend=True, components=True)[source] [github] [bitbucket]

Remove the fitted model from the plot

Also removes the legend by default

clear_all_connections(debug=False) [github] [bitbucket]

Prevent overlapping interactive sessions

clear_highlights() [github] [bitbucket]

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

copy(parent=None, registry=None)[source] [github] [bitbucket]

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)[source] [github] [bitbucket]

When spectrum.crop is called, this must be too

dof

degrees of freedom in fit

downsample(factor)[source] [github] [bitbucket]

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) [github] [bitbucket]

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

firstclick_guess() [github] [bitbucket]

Initialize self.guesses

fitter
fullsizemodel()[source] [github] [bitbucket]

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)[source] [github] [bitbucket]

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)[source] [github] [bitbucket]

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)[source] [github] [bitbucket]

compute the model over the full axis

get_model(xarr, pars=None, debug=False, add_baseline=None)[source] [github] [bitbucket]

Compute the model over a given axis

get_model_frompars(xarr, pars, debug=False, add_baseline=None)[source] [github] [bitbucket]

Compute the model over a given axis

get_model_xlimits(threshold='auto', peak_fraction=0.01, add_baseline=False, unit='pixels', units=None)[source] [github] [bitbucket]

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)[source] [github] [bitbucket]

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) [github] [bitbucket]

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) [github] [bitbucket]

Re-highlight the fitted region

kwargs are passed to matplotlib.plot

history_fitpars()[source] [github] [bitbucket]
integral(analytic=False, direct=False, threshold='auto', integration_limits=None, integration_limit_units='pixels', return_error=False, **kwargs)[source] [github] [bitbucket]

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) np.scalar or np.ndarray with the integral or integral & error :
mask

mask_sliced

Sliced (subset) Mask: True means “exclude”

measure_approximate_fwhm(threshold='error', emission=True, interpolate_factor=1, plot=False, grow_threshold=2, **kwargs)[source] [github] [bitbucket]

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 The approximated FWHM, if it can be computed : If there are <= 2 valid pixels, a fwhm cannot be computed :
model_mask(**kwargs)[source] [github] [bitbucket]

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? 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)[source] [github] [bitbucket]

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)[source] [github] [bitbucket]

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)[source] [github] [bitbucket]

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() 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)[source] [github] [bitbucket]

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)[source] [github] [bitbucket]

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)[source] [github] [bitbucket]

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)[source] [github] [bitbucket]

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)[source] [github] [bitbucket]

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)[source] [github] [bitbucket]

Print the best-fit parameters to the command line

refit(use_lmfit=False)[source] [github] [bitbucket]

Redo a fit using the current parinfo as input

register_fitter(*args, **kwargs)[source] [github] [bitbucket]

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()[source] [github] [bitbucket]

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) [github] [bitbucket]

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)[source] [github] [bitbucket]

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()[source] [github] [bitbucket]

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)[source] [github] [bitbucket]

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) [github] [bitbucket]

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.
xmax
xmin