Line list based Gaussian spectral model

class PyAstronomy.modelSuite.LLGauss(lineList, uniSig=None, modelBinsize=0.005, useFastRB=True, verbose=False, onlyAbs=True)

A spectral model based on Gaussian profiles and a line list.

This class provides a simple spectral model based on a number of Gaussian lines, whose strength may be fitted individually. Note that the EW of the lines is given by: A{n}`*`lineScale, where A{n} is the area of the n-th Gaussian. The scale parameter does not influence the EW of the Gaussians.

Note

The unit of the EWs given in the lineList needs to be the same as the wavelength units.

Fit parameters:

  • lineScale - A common scaling of the area of all lines.

  • scale - A scaling of the entire spectrum.

  • eps - Linear limb-darkening coefficient.

  • vrad - Radial velocity [km/s].

  • vsini - Projected stellar rotational velocity [km/s].

  • A{n} - The amplitudes (area) parameters of the individual Gaussians.

  • sig{n} - The standard deviations of the individual Gaussian.

  • mu{n} - The position of the individual Gaussians.

Parameters
onlyAbsboolean, optional

If True (default), restrictions will be applied, which prevent emission lines in the spectrum.

lineListarray

An array with either two or three columns. The first column given the position of the lines, the second gives the EW of the lines, and the third—if present—gives the depth of the lines. The depth is the maximal depression of the continuum, e.g., a value of 0.96 means that the center of the line of 4% below the continuum. If the depth is given, the width of the individual Gaussians is determined from it, unless the uniSig is specified.

uniSigfloat, optional

Use “unified sigma”, i.e., the same width for all lines. Note that this flag overrules the “depth” column in the lineList, if it has been specified.

modelBinsizefloat, optional

Internally, the model should be calculated on a finer grid than the actual spectrum. This parameter specifies the used bin size, which is 0.005 by default.

useFastRBboolean, optional

Use the “fast” rotational broadening algorithm. This algorithm uses a wavelength-independent broadening kernel, which is considerably faster than considering the wavelength dependence. Setting this flag to False is necessary if you use very long wavelength ranges; by default it is True.

verboseboolean, optional

If True, the class print the current parameters during the evaluation.

Methods

MCMCautoParameters(ranges[, picky, ...])

Convenience function to generate parameters for MCMC fit.

addConditionalRestriction(*args)

Define a conditional restriction.

assignValue(specval)

Assign new values to variables.

assignValues(specval)

Assign new values to variables.

autoFitMCMC(x, y, ranges[, picky, stepsize, ...])

Convenience function to using auto-generated sampling parameters in MCMC.

availableParameters()

Provides a list of existing parameters.

delRestriction(parName)

Delete restriction

description([parenthesis])

Returns a description of the model based on the names of the individual components.

errorConfInterval(par[, dstat, statTol, ...])

Calculate confidence interval for a parameter.

evaluate(x)

Calculates the model for current parameters.

fit(x, y[, yerr, X0, minAlgo, mAA, ...])

Carries out a fit.

fitEMCEE([x, y, yerr, nwalker, priors, ...])

MCMC sampling using emcee package.

fitMCMC(x, y, X0, Lims, Steps[, yerr, ...])

Carry out MCMC fit/error estimation.

freeParamNames()

Get the names of the free parameters.

freeParameters()

Get names and values of free parameters.

freeze(specifiers)

Consider variables free to float.

frozenParameters()

Get names and values of frozen parameters.

getRelationsOf(specifier)

Return relations of a variable.

getRestrictions()

Get all restrictions.

hasVariable(specifier)

Determine whether the variable exists.

numberOfFreeParams()

Get number of free parameters.

numberOfLines()

Get the number of lines in the model.

parameterSummary([toScreen, prefix, sorting])

Writes a summary of the parameters in text form.

parameters()

Obtain parameter names and values.

relate(dependentVar, independentVars[, func])

Define a relation.

removeConditionalRestriction(*args)

Remove an existing conditional constraint.

renameVariable(oldName, newName)

Change name of variable.

restoreState(resource)

Restores parameter values from file or dictionary.

saveState(*args, **kwargs)

Save the state of the fitting object.

setObjectiveFunction([miniFunc])

Define the objective function.

setPenaltyFactor(penalFac)

Change the penalty factor.

setRestriction(restricts)

Define restrictions.

setRootName(root[, rename])

Define the root name of the model.

showConditionalRestrictions(**kwargs)

Show conditional restrictions.

steppar(pars, ranges[, extractFctVal, quiet])

Allows to step a parameter through a specified range.

thaw(specifiers)

Consider variables fixed.

thawLineStrengths([wlmin, wlmax])

Thaw line strengths.

thawLineWidths([wlmin, wlmax])

Thaw line widths.

untie(parName[, forceFree])

Remove all relations of parameter parName, i.e., the parameter is not dependend on other parameters.

updateModel()

Recalculate the model using current settings.
evaluate(x)

Calculates the model for current parameters.

The “model” is calculated on a wavelength axis with binning specified by the modelBinsize parameter.

The line positions are Doppler shifted and the resulting model is rotationally broadened. Finally, the entire model is multiplied by the scale parameter to account for a global normalization.

Parameters
xarray

The wavelengths at which to calculate the model.

Returns
modelarray

The model evaluated at the specified positions.

numberOfLines()

Get the number of lines in the model.

Returns
Number of linesint

Number of Gaussian in the model.

thawLineStrengths(wlmin=None, wlmax=None)

Thaw line strengths.

Thaws parameters of the from A{n}, where n is the number of the Gaussian component. By default all such parameters will be thawed. The selection may, however, be influenced by specifying wlmin and wlmax.

Parameters
wlminfloat, optional

If specified, only the strength of lines at wavelengths larger than this limits will be thawed.

wlmaxfloat, optional

If specified, only the strength of lines at wavelengths below this limit will be thawed.

Returns
Thawed parameterslist

A list of thawed parameter names.

thawLineWidths(wlmin=None, wlmax=None)

Thaw line widths.

Thaws parameters of the from sig{n}, where n is the number of the Gaussian component. By default all such parameters will be thawed. The selection may, however, be influenced by specifying wlmin and wlmax.

Parameters
wlminfloat, optional

If specified, only the strength of lines at wavelengths larger than this limits will be thawed.

wlmaxfloat, optional

If specified, only the strength of lines at wavelengths below this limit will be thawed.

Returns
Thawed parameterslist

A list of thawed parameter names.

Example: Evaluation and fitting

from PyAstronomy import modelSuite as ms
import numpy as np
import matplotlib.pylab as plt

# Create our line list with 4 line
lineList = np.zeros((4, 3))
# Assign wavelengths (in A)
lineList[0, 0] = 5002.37
lineList[1, 0] = 5005.9
lineList[2, 0] = 5007.52
lineList[3, 0] = 5007.64
# Assign EWs (in A)
lineList[0, 1] = 0.01
lineList[1, 1] = 0.05
lineList[2, 1] = 0.009
lineList[3, 1] = 0.12
# Assign depths (0-1)
lineList[0, 2] = 0.97
lineList[1, 2] = 0.9
lineList[2, 2] = 0.99
lineList[3, 2] = 0.35

wvl = np.arange(5000., 5010., 0.01)

# Get an instance of the LLGauss class
llg = ms.LLGauss(lineList)
# Have a look at the model parameters
llg.parameterSummary()
# Evaluate the model
m1 = llg.evaluate(wvl)
# Now apply rotational broadening [km/s]
# with limb-darkening of 0.6
llg["vsini"] = 61.0
llg["eps"] = 0.6
# and evaluate again
mvsini = llg.evaluate(wvl)
# Next, apply a Doppler shift [km/s]
llg["vrad"] = -32.7
# and evaluate
mvrad = llg.evaluate(wvl)
# Plot the results
plt.subplot(2, 1, 1)
plt.plot(wvl, m1, 'b.-')
plt.plot(wvl, mvsini, 'g.-')
plt.plot(wvl, mvrad, 'y.-')

# Now use the model for fitting
# We need "data" ...
data = llg.evaluate(wvl)
# ... with noise
data += np.random.normal(0.0, 0.01, len(data))
# Lets modify the strengths of the Gaussians
# and get it back.
for i in range(llg.numberOfLines()):
    llg["A"+str(i+1)] += np.random.normal(0.0, 0.1)
# Use all line strengths for fitting
llg.thawLineStrengths()
# and fit
llg.fit(wvl, data)
# Plot the result
plt.subplot(2, 1, 2)
plt.errorbar(wvl, data, yerr=np.ones(len(wvl))*0.01, fmt='bp')
plt.plot(wvl, llg.evaluate(wvl), 'r--')
plt.show()