Models¶
See Parameters for information on how to restrict/modify model parameters.
The generic SpectralModel class is a wrapper for model functions. A model should take in an X-axis and some number of parameters. In order to declare a SpectralModel, you give SpectralModel the function name and the number of parameters it requires. The rest of the options are optional, though parnames & shortvarnames are strongly recommended. If you do not specify fitunits, your fitting code must deal with units internally.
Here are some examples of how to make your own fitters:
hill5_fitter = model.SpectralModel(hill5_model, 5,
parnames=['tau', 'v_lsr', 'v_infall', 'sigma', 'tpeak'],
parlimited=[(True,False),(False,False),(True,False),(True,False), (True,False)],
parlimits=[(0,0), (0,0), (0,0), (0,0), (0,0)],
# specify the parameter names (TeX is OK)
shortvarnames=("\\tau","v_{lsr}","v_{infall}","\\sigma","T_{peak}"),
fitunits='Hz' )
gaussfitter = model.SpectralModel(gaussian, 3,
parnames=['amplitude','shift','width'],
parlimited=[(False,False),(False,False),(True,False)],
parlimits=[(0,0), (0,0), (0,0)],
shortvarnames=('A',r'\Delta x',r'\sigma'))
Then you can register these fitters.
Fitting¶
Once you have a model defined, you can fit it using the
pyspeckit.Spectrum.specfit
module. Documents on fitting have not been
prepared yet, but you can learn most of the tricks by looking at the various
fitting examples and the Parameters documentation.
See also Model Fitting.
Generic Models and Tools¶
Specific Models¶
API Documentation for Models¶
We include the API documentation for the generic model and fitter wrappers here.
Module API¶
- class pyspeckit.spectrum.models.model.AstropyModel(model, shortvarnames=None, **kwargs)[source] [github] [bitbucket]¶
Override the SpectralModel initialization
- fitter(xax, data, err=None, quiet=True, veryverbose=False, debug=False, parinfo=None, params=None, npeaks=None, **kwargs)[source] [github] [bitbucket]¶
Run the fitter using mpfit.
kwargs will be passed to _make_parinfo and mpfit.
- Parameters
xax : SpectroscopicAxis
The X-axis of the spectrum
data : ndarray
The data to fit
err : ndarray (optional)
The error on the data. If unspecified, will be uniform unity
parinfo : ParinfoList
The guesses, parameter limits, etc. See
pyspeckit.spectrum.parinfo
for detailsquiet : bool
pass to mpfit. If False, will print out the parameter values for each iteration of the fitter
veryverbose : bool
print out a variety of mpfit output parameters
debug : bool
raise an exception (rather than a warning) if chi^2 is nan
- n_modelfunc(pars=None, debug=False, **kwargs)[source] [github] [bitbucket]¶
Only deals with single-peak functions
- class pyspeckit.spectrum.models.model.SpectralModel(modelfunc, npars, shortvarnames=('A', '\\Delta x', '\\sigma'), fitunit=None, centroid_par=None, fwhm_func=None, fwhm_pars=None, integral_func=None, use_lmfit=False, guess_types=('amplitude', 'center', 'width'), **kwargs)[source] [github] [bitbucket]¶
A wrapper class for a spectra model. Includes internal functions to generate multi-component models, annotations, integrals, and individual components. The declaration can be complex, since you should name individual variables, set limits on them, set the units the fit will be performed in, and set the annotations to be used. Check out some of the hyperfine codes (hcn, n2hp) for examples.
Spectral Model Initialization
Create a Spectral Model class for data fitting
- Parameters
modelfunc : function
the model function to be fitted. Should take an X-axis (spectroscopic axis) as an input followed by input parameters. Returns an array with the same shape as the input X-axis
npars : int
number of parameters required by the model
use_lmfit: bool :
Use lmfit instead of mpfit to do the fitting
parnames : list (optional)
a list or tuple of the parameter names
parvalues : list (optional)
the initial guesses for the input parameters (defaults to ZEROS)
parlimits : list (optional)
the upper/lower limits for each variable (defaults to ZEROS)
parfixed : list (optional)
Can declare any variables to be fixed (defaults to ZEROS)
parerror : list (optional)
technically an output parameter. Specifying it here will have no effect. (defaults to ZEROS)
partied : list (optional)
not the past tense of party. Can declare, via text, that some parameters are tied to each other. Defaults to zeros like the others, but it’s not clear if that’s a sensible default
fitunit : str (optional)
convert X-axis to these units before passing to model
parsteps : list (optional)
minimum step size for each paremeter (defaults to ZEROS)
npeaks : list (optional)
default number of peaks to assume when fitting (can be overridden)
shortvarnames : list (optional)
TeX names of the variables to use when annotating
amplitude_types : tuple
A tuple listing the types of the different parameters when guessing. The valid values are ‘amplitude’, ‘width’, and ‘center’. These are handled by parse_3par_guesses, which translate these into input guess lists for the fitter. For a “standard” 3-parameter Gaussian fitter, nothing changes, but for other models that have more than 3 parameters, some translation is needed.
- Returns
A tuple containing (model best-fit parameters, the model, parameter :
errors, chi^2 value) :
- analytic_centroids(centroidpar=None)[source] [github] [bitbucket]¶
Return the analytic centroids of the model components
- Parameters
centroidpar : None or string
The name of the parameter in the fit that represents the centroid some models have default centroid parameters - these will be used if centroidpar is unspecified
- Returns
List of the centroid values (even if there’s only 1) :
- analytic_fwhm(parinfo=None)[source] [github] [bitbucket]¶
Return the FWHMa of the model components if a fwhm_func has been defined
Done with incomprehensible list comprehensions instead of nested for loops… readability sacrificed for speed and simplicity. This is unpythonic.
- analytic_integral(modelpars=None, npeaks=None, npars=None)[source] [github] [bitbucket]¶
Placeholder for analyic integrals; these must be defined for individual models
- annotations(shortvarnames=None, debug=False)[source] [github] [bitbucket]¶
Return a list of TeX-formatted labels
The values and errors are formatted so that only the significant digits are displayed. Rounding is performed using the decimal package.
- Parameters
shortvarnames : list
A list of variable names (tex is allowed) to include in the annotations. Defaults to self.shortvarnames
Examples
>>> # Annotate a Gaussian >>> sp.specfit.annotate(shortvarnames=['A','\Delta x','\sigma'])
- component_integrals(xarr, dx=None)[source] [github] [bitbucket]¶
Compute the integrals of each component
- components(xarr, pars, **kwargs)[source] [github] [bitbucket]¶
Return a numpy ndarray of shape [npeaks x modelshape] of the independent components of the fits
- computed_centroid(xarr=None)[source] [github] [bitbucket]¶
Return the computed centroid of the model
- Parameters
xarr : None or np.ndarray
The X coordinates of the model over which the centroid should be computed. If unspecified, the centroid will be in pixel units
- fitter(xax, data, err=None, quiet=True, veryverbose=False, debug=False, parinfo=None, **kwargs)[source] [github] [bitbucket]¶
Run the fitter using mpfit.
kwargs will be passed to _make_parinfo and mpfit.
- Parameters
xax : SpectroscopicAxis
The X-axis of the spectrum
data : ndarray
The data to fit
err : ndarray (optional)
The error on the data. If unspecified, will be uniform unity
parinfo : ParinfoList
The guesses, parameter limits, etc. See
pyspeckit.spectrum.parinfo
for detailsquiet : bool
pass to mpfit. If False, will print out the parameter values for each iteration of the fitter
veryverbose : bool
print out a variety of mpfit output parameters
debug : bool
raise an exception (rather than a warning) if chi^2 is nan
- get_emcee_ensemblesampler(xarr, data, error, nwalkers, **kwargs)[source] [github] [bitbucket]¶
Get an emcee walker ensemble for the data & model
- Parameters
data : np.ndarray
error : np.ndarray
nwalkers : int
Number of walkers to use
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') >>> nwalkers = sp.specfit.fitter.npars * 2 >>> emcee_ensemble = sp.specfit.fitter.get_emcee_ensemblesampler(sp.xarr, sp.data, sp.error, nwalkers) >>> p0 = np.array([sp.specfit.parinfo.values] * nwalkers) >>> p0 *= np.random.randn(*p0.shape) / 10. + 1.0 >>> pos,logprob,state = emcee_ensemble.run_mcmc(p0,100)
- get_emcee_sampler(xarr, data, error, **kwargs)[source] [github] [bitbucket]¶
Get an emcee walker for the data & model
- Parameters
xarr : pyspeckit.units.SpectroscopicAxis
data : np.ndarray
error : np.ndarray
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_sampler = sp.specfit.fitter.get_emcee_sampler(sp.xarr, sp.data, sp.error) >>> p0 = sp.specfit.parinfo >>> emcee_sampler.run_mcmc(p0,100)
- get_pymc(xarr, data, error, use_fitted_values=False, inf=inf, use_adaptive=False, return_dict=False, **kwargs)[source] [github] [bitbucket]¶
Create a pymc MCMC sampler. Defaults to ‘uninformative’ priors
- 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
use_adaptive : bool
Use the Adaptive Metropolis-Hastings sampler?
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.fitter.get_pymc(sp.xarr, sp.data, sp.error) >>> MCwithpriors = sp.specfit.fitter.get_pymc(sp.xarr, sp.data, sp.error, 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']
- integral(modelpars, dx=None, **kwargs)[source] [github] [bitbucket]¶
Extremely simple integrator: IGNORES modelpars; just sums self.model
- lmfitfun(x, y, err=None, debug=False)[source] [github] [bitbucket]¶
Wrapper function to compute the fit residuals in an lmfit-friendly format
- lmfitter(xax, data, err=None, parinfo=None, quiet=True, debug=False, **kwargs)[source] [github] [bitbucket]¶
Use lmfit instead of mpfit to do the fitting
- Parameters
xax : SpectroscopicAxis
The X-axis of the spectrum
data : ndarray
The data to fit
err : ndarray (optional)
The error on the data. If unspecified, will be uniform unity
parinfo : ParinfoList
The guesses, parameter limits, etc. See
pyspeckit.spectrum.parinfo
for detailsquiet : bool
If false, print out some messages about the fitting
- logp(xarr, data, error, pars=None)[source] [github] [bitbucket]¶
Return the log probability of the model. If the parameter is out of range, return -inf
- mpfitfun(x, y, err=None)[source] [github] [bitbucket]¶
Wrapper function to compute the fit residuals in an mpfit-friendly format
- n_modelfunc(pars=None, debug=False, **kwargs)[source] [github] [bitbucket]¶
Simple wrapper to deal with N independent peaks for a given spectral model
- parse_3par_guesses(guesses)[source] [github] [bitbucket]¶
Try to convert a set of interactive guesses (peak, center, width) into guesses appropriate to the model.
- slope(xinp)[source] [github] [bitbucket]¶
Find the local slope of the model at location x (x must be in xax’s units)
- pyspeckit.spectrum.models.model.parse_offset_guess(gname, gval)[source] [github] [bitbucket]¶
Utility function for handling guesses. Allows guess types to be specified as ‘amplitude*2’ or ‘width+3’.
Adds a variable height (background) component to any model