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.

Related Documents

• Parameters

Generic Models and Tools

• Gaussian model
• Voigt Profile model
• LTE Molecule Model
• Hyperfine Line model
• Model Grid
• Spectral Template Model

Specific Models

• Ammonia Temperature and Hyperfine model
• Formaldehyde model
• HCN model
• Hill5 Infall Model
• N2H+ Model
• Ionized Hydrogen 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)

Override the SpectralModel initialization

fitter(xax, data, err=None, quiet=True, veryverbose=False, debug=False, parinfo=None, params=None, npeaks=None, **kwargs)

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 details

quiet : 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)

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)

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)

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)

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)

Placeholder for analyic integrals; these must be defined for individual models

annotations(shortvarnames=None, debug=False)

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)

Compute the integrals of each component

components(xarr, pars, **kwargs)

Return a numpy ndarray of shape [npeaks x modelshape] of the independent components of the fits

computed_centroid(xarr=None)

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)

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 details

quiet : 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)

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)

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)

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)

Extremely simple integrator: IGNORES modelpars; just sums self.model

lmfitfun(x, y, err=None, debug=False)

Wrapper function to compute the fit residuals in an lmfit-friendly format

lmfitter(xax, data, err=None, parinfo=None, quiet=True, debug=False, **kwargs)

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 details

quiet : bool

If false, print out some messages about the fitting

logp(xarr, data, error, pars=None)

Return the log probability of the model. If the parameter is out of range, return -inf

mpfitfun(x, y, err=None)

Wrapper function to compute the fit residuals in an mpfit-friendly format

n_modelfunc(pars=None, debug=False, **kwargs)

Simple wrapper to deal with N independent peaks for a given spectral model

parse_3par_guesses(guesses)

Try to convert a set of interactive guesses (peak, center, width) into guesses appropriate to the model.

slope(xinp)

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)

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

Module API

[edit this page on github]

[edit this page on bitbucket]