Model Fitting

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

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')

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

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

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

mask

Mask: True means “exclude”

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?

Other Parameters:  

override: True | False

Whether to override any existing type if already present.

key: char

Key to select the fitter in interactive mode

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.

xmax

xmin

• Fitting a user-defined model to a spectrum

[edit this page on github]

[edit this page on bitbucket]