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:
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')[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
- Returns
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
- property 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
- property 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 orthreshold
otherwisethreshold : ‘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 detailsintegration_limits : None or 2-tuple
Manually specify the limits in
integration_limit_units
unitsreturn_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
ifnot(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)[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
- Returns
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 peakerror: 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 atadd_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)[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
structureGuess 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()
- 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)[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 plotterexclude: {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.
- property xmax¶
- property xmin¶