Extract galaxy stellar kinematics (V, sigma, h3, h4, h5, h6,...) or the stellar population and gas emission by fitting a template to an observed spectrum in pixel space, using the Penalized PiXel-Fitting (pPXF) method originally described in
and significantly upgraded in
The following key optional features are also available:
from ppxf.ppxf import ppxf
pp = ppxf(templates, galaxy, noise, velscale, start,
bias=None, bounds=None, clean=False, component=0, degree=4,
fixed=None, fraction=None, ftol=1e-4, gas_component=None,
gas_names=None, gas_reddening=None, goodpixels=None, lam=None,
linear=False, mask=None, method='capfit', mdegree=0, moments=2,
plot=False, quiet=False, reddening=None, reddening_func=None,
reg_ord=2, reg_dim=None, regul=0, sigma_diff=0, sky=None,
templates_rfft=None, tied=None, trig=False, velscale_ratio=1,
vsyst=0)
print(pp.sol) # print best-fitting kinematics (V, sigma, h3, h4)
pp.plot() # Plot best fit and gas linesVector containing the spectrum of a single template star or more commonly an array of dimensions templates[nPixels, nTemplates] containing different templates to be optimized during the fit of the kinematics. It has to be nPixels >= galaxy.size.
To apply linear regularization to the weights via the keyword regul, templates should be an array of two templates[nPixels, nAge], three templates[nPixels, nAge, nMetal] or four templates[nPixels, nAge, nMetal, nAlpha] dimensions, depending on the number of population variables one wants to study. This can be useful to try to attach a physical meaning to the output weights, in term of the galaxy star formation history and chemical composition distribution. In that case the templates may represent single stellar population SSP models and should be arranged in sequence of increasing age, metallicity or alpha along the second, third or fourth dimension of the array respectively.
Vector containing the spectrum of the galaxy to be measured. The star and the galaxy spectra have to be logarithmically rebinned but the continuum should not be subtracted. The rebinning may be performed with the log_rebin routine that is distributed with pPXF.
For high redshift galaxies, one should bring the spectra close to the restframe wavelength, before doing the pPXF fit. This can be done by dividing the observed wavelength by (1 + z), where z is a rough estimate of the galaxy redshift, before the logarithmic rebinning. See Section 2.4 of Cappellari (2017) for details.
galaxy can also be an array of dimensions galaxy[nGalPixels, 2] containing two spectra to be fitted, at the same time, with a reflection-symmetric LOSVD. This is useful for spectra taken at point-symmetric spatial positions with respect to the center of an equilibrium stellar system. For a discussion of the usefulness of this two-sided fitting see e.g. Section 3.6 of Rix & White (1992).
IMPORTANT: (1) For the two-sided fitting the vsyst keyword has to be used. (2) Make sure the spectra are rescaled to be not too many order of magnitude different from unity, to avoid over or underflow problems in the calculation. E.g. units of erg/(s cm^2 A) may cause problems!
Vector containing the 1*sigma error (per pixel) in the galaxy spectrum, or covariance matrix describing the correlated errors in the galaxy spectrum. Of course this vector/matrix must have the same units as the galaxy spectrum.
If galaxy is a Nx2 array, noise has to be an array with the same dimensions.
When noise has dimensions NxN it is assumed to contain the covariance matrix with elements sigma(i, j). When the errors in the spectrum are uncorrelated it is mathematically equivalent to input in pPXF an error vector noise=errvec or a NxN diagonal matrix noise=np.diag(errvec**2) (note squared!).
IMPORTANT: the penalty term of the pPXF method is based on the relative change of the fit residuals. For this reason the penalty will work as expected even if no reliable estimate of the noise is available (see Cappellari & Emsellem [2004] for details). If no reliable noise is available this keyword can just be set to:
noise = np.ones_like(galaxy) # Same weight for all pixelsVelocity scale of the spectra in km/s per pixel. It has to be the same for both the galaxy and the template spectra. An exception is when the velscale_ratio keyword is used, in which case one can input templates with smaller velscale than galaxy.
velscale is defined in pPXF by velscale = c*Delta[np.log(lambda)], which is approximately velscale ~ c*Delta(lambda)/lambda. See Section 2.3 of Cappellari (2017) for details.
Vector, or list/array of vectors [start1, start2, ...], with the initial estimate for the LOSVD parameters.
When LOSVD parameters are not held fixed, each vector only needs to contain start = [velStart, sigmaStart] the initial guess for the velocity and the velocity dispersion in km/s. The starting values for h3-h6 (if they are fitted) are all set to zero by default. In other words, when moments=4:
start = [velStart, sigmaStart]is interpreted as:
start = [velStart, sigmaStart, 0, 0]When the LOSVD for some kinematic components is held fixed (see fixed keyword), all values for [Vel, Sigma, h3, h4,...] can be provided.
Unless a good initial guess is available, it is recommended to set the starting sigma >= 3*velscale in km/s (i.e. 3 pixels). In fact when the sigma is very low, and far from the true solution, the chi^2 of the fit becomes weakly sensitive to small variations in sigma (see pPXF paper). In some instances the near-constancy of chi^2 may cause premature convergence of the optimization.
In the case of two-sided fitting a good starting value for the velocity is velStart = 0.0 (in this case vsyst will generally be nonzero). Alternatively on should keep in mind that velStart refers to the first input galaxy spectrum, while the second will have velocity -velStart.
With multiple kinematic components start must be a list of starting values, one for each different component.
EXAMPLE: We want to fit two kinematic components. We fit 4 moments for the first component and 2 moments for the second one as follows:
component = [0, 0, ... 0, 1, 1, ... 1]
moments = [4, 2]
start = [[V1, sigma1], [V2, sigma2]]This parameter biases the (h3, h4, ...) measurements towards zero (Gaussian LOSVD) unless their inclusion significantly decreases the error in the fit. Set this to bias=0 not to bias the fit: the solution (including [V, sigma]) will be noisier in that case. The default bias should provide acceptable results in most cases, but it would be safe to test it with Monte Carlo simulations. This keyword precisely corresponds to the parameter lambda in the Cappellari & Emsellem (2004) paper. Note that the penalty depends on the relative change of the fit residuals, so it is insensitive to proper scaling of the noise vector. A nonzero bias can be safely used even without a reliable noise spectrum, or with equal weighting for all pixels.
Lower and upper bounds for every kinematic parameter. This is an array, or list of arrays, with the same dimensions as start, except for the last one, which is two. In practice, for every elements of start one needs to specify a pair of values [lower, upper].
EXAMPLE: We want to fit two kinematic components, with 4 moments for the first component and 2 for the second (e.g. stars and gas). In this case:
moments = [4, 2]
start_stars = [V1, sigma1, 0, 0]
start_gas = [V2, sigma2]
start = [start_stars, start_gas]then we can specify boundaries for each kinematic parameter as:
bounds_stars = [[V1_lo, V1_up], [sigma1_lo, sigma1_up],
[-0.3, 0.3], [-0.3, 0.3]]
bounds_gas = [[V2_lo, V2_up], [sigma2_lo, sigma2_up]]
bounds = [bounds_stars, bounds_gas]When fitting more than one kinematic component, this keyword should contain the component number of each input template. In principle every template can belong to a different kinematic component.
EXAMPLE: We want to fit the first 50 templates to component 0 and the last 10 templates to component 1. In this case:
component = [0]*50 + [1]*10which, in Python syntax, is equivalent to:
component = [0, 0, ... 0, 1, 1, ... 1]This keyword is especially useful when fitting both emission (gas) and absorption (stars) templates simultaneously (see example for moments keyword).
Set this keyword to use the iterative sigma clipping method described in Section 2.1 of Cappellari et al. (2002). This is useful to remove from the fit unmasked bad pixels, residual gas emissions or cosmic rays.
IMPORTANT: This is recommended only if a reliable estimate of the noise spectrum is available. See also note below for .chi2.
Degree of the additive Legendre polynomial used to correct the template continuum shape during the fit (default: 4). Set degree=-1 not to include any additive polynomial.
Boolean vector set to True where a given kinematic parameter has to be held fixed with the value given in start. This is an array, or list, with the same dimensions as start.
EXAMPLE: We want to fit two kinematic components, with 4 moments for the first component and 2 for the second. In this case:
moments = [4, 2]
start = [[V1, sigma1, h3, h4], [V2, sigma2]]then we can held fixed e.g. the sigma (only) of both components using:
fixed = [[0, 1, 0, 0], [0, 1]]NOTE: Setting a negative moments for a kinematic component is entirely equivalent to setting fixed = 1 for all parameters of the given kinematic component. In other words:
moments = [-4, 2]is equivalent to:
moments = [4, 2]
fixed = [[1, 1, 1, 1], [0, 0]]This keyword allows one to fix the ratio between the first two kinematic components. This is a scalar defined as follows:
fraction = np.sum(weights[component == 0])
/ np.sum(weights[component < 2])This is useful e.g. to try to kinematically decompose bulge and disk.
IMPORTANT: The templates and galaxy spectra should be normalized with mean ~ 1 (as order of magnitude) for the fraction keyword to work as expected. A warning is printed if this is not the case and the resulting output fraction is inaccurate.
The remaining kinematic components (component > 1) are left free, and this allows, for example, to still include gas emission line components.
Fractional tolerance for stopping the non-linear minimization.
Boolean vector, of the same size as component, set to True where the given component describes a gas emission line. If given, pPXF provides the pp.gas_flux and pp.gas_flux_error in output.
EXAMPLE: In the common situation where component = 0 are stellar templates and the rest are gas emission lines, one will set:
gas_component = component > 0This keyword is also used to plot the gas lines with a different color.
String array specifying the names of the emission lines (e.g. gas_names=["Hbeta", "[OIII]",...], one per gas line. The length of this vector must match the number of nonzero elements in gas_component. This vector is only used by pPXF to print the line names on the console.
Set this keyword to an initial estimate of the gas reddening E(B-V) >= 0 to fit a positive reddening together with the kinematics and the templates. The fit assumes by default the extinction curve of Calzetti et al. (2000) but any other prescription can be passed via the reddening_func keyword.
Integer vector containing the indices of the good pixels in the galaxy spectrum (in increasing order). Only these spectral pixels are included in the fit.
IMPORTANT: in all likely situations this keyword has to be specified.
When the keyword reddening or gas_reddening are used, the user has to pass in this keyword a vector with the same dimensions of galaxy, giving the restframe wavelength in Angstrom of every pixel in the input galaxy spectrum. If one uses my log_rebin routine to rebin the spectrum before the pPXF fit:
from ppxf.ppxf_util import log_rebin
specNew, logLam, velscale = log_rebin(lamRange, galaxy)the wavelength can be obtained as lam = np.exp(logLam).
When lam is given, the wavelength is shown in the best-fitting plot, instead of the pixels.
Set to True to keep all nonlinear parameters fixed and only perform a linear fit for the templates and additive polynomials weights. The output solution is a copy of the input one and the errors are zero.
Boolean vector of length galaxy.size specifying with True the pixels that should be included in the fit. This keyword is just an alternative way of specifying the goodpixels.
Algorithm to perform the non-linear minimization step. The default 'capfit' is a novel trust-region implementation of the Levenberg-Marquardt method, generalized to deal with box constraints in a rigorous manner, while also allowing for tied or fixed variables. See documentation of scipy.optimize.least_squares for the other methods.
Degree of the multiplicative Legendre polynomial (with a mean of 1) used to correct the continuum shape during the fit (default: 0). The zero degree multiplicative polynomial is always included in the fit as it corresponds to the weights assigned to the templates. Note that the computation time is longer with multiplicative polynomials than with the same number of additive polynomials.
IMPORTANT: Multiplicative polynomials cannot be used when the reddening keyword is set, as they are degenerate with the reddening.
Order of the Gauss-Hermite moments to fit. Set this keyword to 4 to fit [h3, h4] and to 6 to fit [h3, h4, h5, h6]. Note that in all cases the G-H moments are fitted (non-linearly) together with [V, sigma].
If moments=2 or moments is not set then only [V, sigma] are fitted.
If moments is negative then the kinematics of the given component are kept fixed to the input values. NOTE: Setting a negative moments for a kinematic component is entirely equivalent to setting fixed = 1 for all parameters of the given kinematic component.
EXAMPLE: We want to keep fixed component = 0, which has a LOSVD described by [V, sigma, h3, h4] and is modelled with 100 spectral templates; At the same time, we fit [V, sigma] for component = 1, which is described by 5 templates (this situation may arise when fitting stellar templates with pre-determined stellar kinematics, while fitting the gas emission). We should give in input to pPXF the following parameters:
component = [0]*100 + [1]*5 # --> [0, 0, ... 0, 1, 1, 1, 1, 1]
moments = [-4, 2]
start = [[V, sigma, h3, h4], [V, sigma]]Integer. Gives the integer ratio >= 1 between the velscale of the galaxy and the templates. When this keyword is used, the templates are convolved by the LOSVD at their native resolution, and only subsequently are integrated over the pixels and fitted to galaxy. This keyword is generally unnecessary and mostly useful for testing.
Note that in realistic situations the uncertainty in the knowledge and variations of the intrinsic line-spread function become the limiting factor in recovering the LOSVD well below velscale.
Set this keyword to plot the best fitting solution and the residuals at the end of the fit.
One can also call the class function pp.plot() after the call to pp = ppxf(...).
Set this keyword to suppress verbose output of the best fitting parameters at the end of the fit.
Set this keyword to an initial estimate of the reddening E(B-V) >= 0 to fit a positive reddening together with the kinematics and the templates. The fit assumes by default the extinction curve of Calzetti et al. (2000, ApJ, 533, 682) but any other prescription can be passed via the reddening_func keyword.
IMPORTANT: The mdegree keyword cannot be used when reddening is set.
If this keyword is nonzero, the program applies first or second order linear regularization to the weights during the pPXF fit. Regularization is done in one, two or three dimensions depending on whether the array of templates has two, three or four dimensions respectively. Large regul values correspond to smoother weights output. When this keyword is nonzero the solution will be a trade-off between smoothness of weights and goodness of fit.
Section 3.5 of Cappellari (2017) gives a description of regularization.
When fitting multiple kinematic component the regularization is applied only to the first component = 0, while additional components are not regularized. This is useful when fitting stellar population together with gas emission lines. In that case, the SSP spectral templates must be given first and the gas emission templates are given last. In this situation one has to use the reg_dim keyword (below), to give pPXF the dimensions of the population parameters (e.g. nAge, nMetal, nAlpha). An usage example is given in the file ppxf_example_population_gas_sdss.py.
The effect of the regularization scheme is the following:
reg_ord=1 it enforces the numerical first derivatives between neighbouring weights (in the 1-dim case) to be equal to w[j] - w[j+1] = 0 with an error Delta = 1/regul.reg_ord=2 it enforces the numerical second derivatives between neighboring weights (in the 1-dim case) to be equal to w[j-1] - 2*w[j] + w[j+1] = 0 with an error Delta = 1/regul.It may be helpful to define regul = 1/Delta and think of Delta as the regularization error.
IMPORTANT: Delta needs to be smaller but of the same order of magnitude of the typical weights to play an effect on the regularization. One quick way to achieve this is:
Divide the full templates array by a scalar in such a way that the typical template has a median of one:
templates /= np.median(templates)Do the same for the input galaxy spectrum:
galaxy /= np.median(galaxy)Delta will be a few percent (e.g. 0.01 --> regul=100).Alternatively, for a more rigorous definition of the parameter regul:
Perform an un-regularized fit (regul=0) and then rescale the input noise spectrum so that:
Chi^2/DOF = Chi^2/goodPixels.size = 1.This is achieved by rescaling the input noise spectrum as:
noise = noise*sqrt(Chi**2/DOF) = noise*sqrt(pp.chi2);regul and iteratively redo the pPXF fit until the Chi^2 increases from the unregularized value Chi^2 = goodPixels.size value by DeltaChi^2 = sqrt(2*goodPixels.size).The derived regularization corresponds to the maximum one still consistent with the observations and the derived star formation history will be the smoothest (minimum curvature or minimum variation) that is still consistent with the observations.
Order of the derivative that is minimized by the regularization. The following two rotationally-symmetric estimators are supported:
reg_ord=1: minimizes the integral over the weights of the squared gradient:
Grad[w] @ Grad[w].reg_ord=2: minimizes the integral over the weights of the squared curvature:
Laplacian[w]**2.When using regularization with more than one kinematic component (using the component keyword), the regularization is only applied to the first one (component=0). This is useful to fit the stellar population and gas emission together.
In this situation, one has to use the reg_dim keyword, to give pPXF the dimensions of the population parameters (e.g. nAge, nMetal, nAlpha). One should creates the initial array of population templates like e.g. templates[nPixels, nAge, nMetal, nAlpha] and define:
reg_dim = templates.shape[1:] # = [nAge, nMetal, nAlpha]The array of stellar templates is then reshaped into a 2-dim array as:
templates = templates.reshape(templates.shape[0], -1)and the gas emission templates are appended as extra columns at the end. An usage example is given in ppxf_example_population_gas_sdss.py.
When using regularization with a single component (the component keyword is not used, or contains identical values), the number of population templates along different dimensions (e.g. nAge, nMetal, nAlpha) is inferred from the dimensions of the templates array and this keyword is not necessary.
Quadratic difference in km/s defined as:
sigma_diff**2 = sigma_inst**2 - sigma_temp**2between the instrumental dispersion of the galaxy spectrum and the instrumental dispersion of the template spectra.
This keyword is useful when the templates have higher resolution than the galaxy and they were not convolved to match the instrumental dispersion of the galaxy spectrum. In this situation, the convolution is done by pPXF with increased accuracy, using an analytic Fourier Transform.
vector containing the spectrum of the sky to be included in the fit, or array of dimensions sky[nPixels, nSky] containing different sky spectra to add to the model of the observed galaxy spectrum. The sky has to be log-rebinned as the galaxy spectrum and needs to have the same number of pixels.
The sky is generally subtracted from the data before the pPXF fit. However, for observations very heavily dominated by the sky spectrum, where a very accurate sky subtraction is critical, it may be useful not to subtract the sky from the spectrum, but to include it in the fit using this keyword.
When calling pPXF many times with an identical set of templates, one can use this keyword to pass the real FFT of the templates, computed in a previous pPXF call, stored in the pp.templates_rfft attribute. This keyword mainly exists to show that there is no need for it...
IMPORTANT: Use this keyword only if you understand what you are doing!
A list of string expressions. Each expression "ties" the parameter to other free or fixed parameters. Any expression involving constants and the parameter array p[j] are permitted. Since they are totally constrained, tied parameters are considered to be fixed; no errors are computed for them.
This is an array, or list of arrays, with the same dimensions as start. In practice, for every element of start one needs to specify either an empty string '' implying that the parameter is free, or a string expression involving some of the variables p[j], where the index j represents the index of the flattened list of kinematic parameters.
EXAMPLE: We want to fit three kinematic components, with 4 moments for the first component and 2 moments for the second and third (e.g. stars and two gas components). In this case:
moments = [4, 2, 2]
start = [[V1, sigma1, 0, 0], [V2, sigma2], [V3, sigma3]]then we can force the equality constraint V2 = V3 as follows:
tied = [['', '', '', ''], ['', ''], ['p[4]', '']]or we can force the equality constraint sigma2 = sigma3 as follows:
tied = [['', '', '', ''], ['', ''], ['', 'p[5]']]NOTE: One could in principle use the tied keyword to completely tie the LOSVD of two kinematic components. However this same effect is more efficient achieved by assigning them to the same kinematic component using the component keyword.
Set trig=True to use trigonometric series as an alternative to Legendre polynomials, for both the additive and multiplicative polynomials. When trig=True the fitted series below has N = degree/2 or N = mdegree/2:
poly = A_0 + sum_{n=1}^{N} [A_n*cos(n*th) + B_n*sin(n*th)]IMPORTANT: The trigonometric series has periodic boundary conditions. This is sometimes a desirable property, but this expansion is not as flexible as the Legendre polynomials.
Reference velocity in km/s (default 0). The input initial guess and the output velocities are measured with respect to this velocity. This keyword is generally used to account for the difference in the starting wavelength of the templates and the galaxy spectrum as follows:
vsyst = c*np.log(wave_temp[0]/wave_gal[0])The value assigned to this keyword is crucial for the two-sided fitting. In this case vsyst can be determined from a previous normal one-sided fit to the galaxy velocity profile. After that initial fit, vsyst can be defined as the measured velocity at the galaxy center. More accurately vsyst is the value which has to be subtracted to obtain a nearly anti-symmetric velocity profile at the two opposite sides of the galaxy nucleus.
IMPORTANT: this value is generally different from the systemic velocity one can get from the literature. Do not try to use that!
Stored as attributes of the pPXF class:
Vector with the best fitting additive polynomial.
Vector with the best fitting model for the galaxy spectrum. This is a linear combination of the templates, convolved with the best fitting LOSVD, multiplied by the multiplicative polynomials and with subsequently added polynomial continuum terms or sky components.
The reduced chi^2 (namely chi^2/DOF) of the fit (where DOF ~ pp.goodpixels.size).
IMPORTANT: if Chi^2/DOF is not ~1 it means that the errors are not properly estimated, or that the template is bad and it is not safe to set the clean keyword.
If gas_component is not None, this attribute returns the best-fitting gas spectrum alone. The stellar spectrum alone can be computed as stellar_spectrum = pp.bestfit - pp.gas_bestfit
Vector with the integrated flux (in counts) of all lines set as True in the input gas_component keyword. If a line is composed of a doublet, the flux is that of both lines. If the Balmer series is input as a single template, this is the flux of the entire series.
IMPORTANT: pPXF makes no assumptions about the input flux units: The returned .gas_flux has the same units as the value one would obtain by just summing the values of the spectral pixels of the gas emission. This implies that, if the spectrum is in units of erg/(cm^2 s A), the pPXF value should be multiplied by the pixel size in Angstrom at the line wavelength to obtain the integrated line flux in units of erg/(cm^2 s).
NOTE: If there is no gas reddening and each input gas templates was normalized to sum = 1, then pp.gas_flux = pp.weights[pp.gas_component].
When a gas template is identically zero within the fitted region, then pp.gas_flux = pp.gas_flux_error = np.nan. The corresponding components of pp.gas_zero_template are set to True. These np.nan values are set at the end of the calculation to flag the undefined values. They do not indicate numerical issues with the actual pPXF calculation, and the rest of the pPXF output is reliable.
formal uncertainty (1*sigma) for the quantity pp.gas_flux.
This error is approximate as it ignores the covariance between the gas flux and any non-linear parameter. Bootstrapping can be used for more accurate errors.
These errors are meaningless unless Chi^2/DOF ~ 1. However if one assumes that the fit is good, a corrected estimate of the errors is:
gas_flux_error_corr = gas_flux_error*sqrt(chi^2/DOF)
= pp.gas_flux_error*sqrt(pp.chi2).vector with the best-fitting gas reddening curve.
Best fitting E(B-V) value if the gas_reddening keyword is set.
vector of size gas_component.sum() set to True where the gas template was identically zero within the fitted region. For those gas components pp.gas_flux = pp.gas_flux_error = np.nan.
Integer vector containing the indices of the good pixels in the fit. This vector is a copy of the input goodpixels if clean = False otherwise it will be updated by removing the detected outliers.
This variable contains a vector of formal uncertainty (1*sigma) for the fitted parameters in the output vector sol. This option can be used when speed is essential, to obtain an order of magnitude estimate of the uncertainties, but we strongly recommend to run bootstrapping simulations to obtain more reliable errors. In fact these errors can be severely underestimated in the region where the penalty effect is most important (sigma < 2*velscale).
These errors are meaningless unless Chi^2/DOF ~ 1. However if one assumes that the fit is good, a corrected estimate of the errors is:
error_corr = error*sqrt(chi^2/DOF) = pp.error*sqrt(pp.chi2).IMPORTANT: when running Monte Carlo simulations to determine the error, the penalty (bias) should be set to zero, or better to a very small value. See Section 3.4 of Cappellari & Emsellem (2004) for an explanation.
This is largely superseded by the .apoly attribute above.
When degree >= 0 contains the weights of the additive Legendre polynomials of order 0, 1,... degree. The best fitting additive polynomial can be explicitly evaluated as:
from numpy.polynomial import legendre
x = np.linspace(-1, 1, len(galaxy))
apoly = legendre.legval(x, pp.polyweights)When trig=True the polynomial is evaluated as:
apoly = pp.trigval(x, pp.polyweights)When doing a two-sided fitting (see help for galaxy parameter), the additive polynomials are allowed to be different for the left and right spectrum. In that case the output weights of the additive polynomials alternate between the first (left) spectrum and the second (right) spectrum.
Prediction matrix[nPixels, degree+nTemplates] of the linear system.
pp.matrix[nPixels, :degree] contains the additive polynomials if degree >= 0.
pp.matrix[nPixels, degree:] contains the templates convolved by the LOSVD, and multiplied by the multiplicative polynomials if mdegree > 0.
pp.matrix[nPixels, -nGas:] contains the nGas emission line templates if given. In the latter case the best fitting gas emission line spectrum is:
lines = pp.matrix[:, -nGas:] @ pp.weights[-nGas:]Best fitting multiplicative polynomial (or reddening curve when reddening is set).
This is largely superseded by the .mpoly attribute above.
When mdegree > 0 this contains in output the coefficients of the multiplicative Legendre polynomials of order 1, 2,... mdegree. The polynomial can be explicitly evaluated as:
from numpy.polynomial import legendre
x = np.linspace(-1, 1, len(galaxy))
mpoly = legendre.legval(x, np.append(1, pp.mpolyweights))When trig = True the polynomial is evaluated as:
mpoly = pp.trigval(x, np.append(1, pp.mpolyweights))Best fitting E(B-V) value if the reddening keyword is set.
Vector containing in output the parameters of the kinematics.
moments=2 this contains [Vel, Sigma]moments=4 this contains [Vel, Sigma, h3, h4]moments=N this contains [Vel, Sigma, h3,... hN]When fitting multiple kinematic component, pp.sol contains a list with the solution for all different components, one after the other, sorted by component: pp.sol = [sol1, sol2,...].
Vel is the velocity, Sigma is the velocity dispersion, h3 - h6 are the Gauss-Hermite coefficients. The model parameters are fitted simultaneously.
IMPORTANT: The precise relation between the output pPXF velocity and redshift is Vel = c*np.log(1 + z). See Section 2.3 of Cappellari (2017) for a detailed explanation.
These are the default safety limits on the fitting parameters (they can be changed using the bounds keyword):
velscale/100 < Sigma < 1000 km/s-0.3 < [h3, h4, ...] < 0.3 (extreme value for real galaxies)In the case of two-sided LOSVD fitting the output values refer to the first input galaxy spectrum, while the second spectrum will have by construction kinematics parameters [-Vel, Sigma, -h3, h4, -h5, h6]. If vsyst is nonzero (as required for two-sided fitting), then the output velocity is measured with respect to vsyst.
Contains the output status of the optimization. Positive values generally represent success (the status is defined as in scipy.optimize.least_squares).
Receives the value of the weights by which each template was multiplied to best fit the galaxy spectrum. The optimal template can be computed with an array-vector multiplication:
bestemp = templates @ weightsThese weights do not include the weights of the additive polynomials which are separately stored in pp.polyweights.
When the sky keyword is used weights[:nTemplates] contains the weights for the templates, while weights[nTemplates:] gives the ones for the sky. In that case the best fitting galaxy template and sky are given by:
bestemp = templates @ weights[:nTemplates]
bestsky = sky @ weights[nTemplates:]When doing a two-sided fitting (see help for galaxy parameter) together with the sky keyword, the sky weights are allowed to be different for the left and right spectrum. In that case the output sky weights alternate between the first (left) spectrum and the second (right) spectrum.
The pPXF routine can give sensible quick results with the default bias parameter, however, like in any penalized/filtered/regularized method, the optimal amount of penalization generally depends on the problem under study.
The general rule here is that the penalty should leave the line-of-sight velocity-distribution (LOSVD) virtually unaffected, when it is well sampled and the signal-to-noise ratio (S/N) is sufficiently high.
EXAMPLE: If you expect an LOSVD with up to a high h4 ~ 0.2 and your adopted penalty (bias) biases the solution towards a much lower h4 ~ 0.1, even when the measured sigma > 3*velscale and the S/N is high, then you are misusing the pPXF method!
THE RECIPE: The following is a simple practical recipe for a sensible determination of the penalty in pPXF:
(S/N)_min level for your kinematics extraction and spatially bin your data so that there are no spectra below (S/N)_min;bias=0). The solution will be noisy and may be affected by spurious solutions, however this step will allow you to check the expected mean ranges in the Gauss-Hermite parameters [h3, h4] for the galaxy under study;ppxf_example_simulation.py routine. Adopt as S/N in the simulation the chosen value (S/N)_min and as input [h3, h4] the maximum representative values measured in the non-penalized pPXF fit of the previous step;bias) the largest value such that, for sigma > 3*velscale, the mean difference delta between the output [h3, h4] and the input [h3, h4] is well within (e.g. delta ~ rms/3) the rms scatter of the simulated values (see an example in Fig. 2 of Emsellem et al. 2004).Common problems with your first pPXF fit are caused by incorrect wavelength ranges or different velocity scales between galaxy and templates. To quickly detect these problems try to overplot the (log rebinned) galaxy and the template just before calling the pPXF procedure.
You can use something like the following Python lines while adjusting the smoothing window and the pixels shift. If you cannot get a rough match by eye it means something is wrong and it is unlikely that pPXF (or any other program) will find a good match:
import numpy as np
import matplotlib.pyplot as plt
import scipy.ndimage
sigma = 2 # Velocity dispersion in pixels
shift = -20 # Velocity shift in pixels
template = np.roll(ndimage.gaussian_filter1d(template, sigma), shift)
plt.plot(galaxy, 'k')
plt.plot(template*np.median(galaxy)/np.median(template), 'r')reddening_cal00(). It only affected NIR lam > 1000 nm.start is not longer than its moments.ftol keyword..gas_reddening when mdegree > 0.misc.factorial, deprecated in Scipy 1.0.gas_flux_error normalization, when input not normalized..gas_bestfit, .gas_mpoly, .mpoly and .apoly attributes.gas_reddening (useful with tied Balmer emission lines).axvspan to visualize masked regions in plot.linear keyword.reddening_func keyword.galaxy and check all bounds have two elements.start to be either a list or an array or vectors..gas_flux output units. Thanks to Xihan Ji (Tsinghua University) for the feedback.gas_component and gas_names keywords.MPFIT with CAPFIT for a Levenberg-Marquardt method with fixed or tied variables, which rigorously accounts for box constraints.method keyword to use Scipy's least_squares() as alternative to MPFIT.start and bounds.linear_fit() does not return unused status any more, for consistency with the correspinding change to cap_mpfit.tied keyword to tie parameters during fitting..weights._templates_rfft(). Many thanks to Jesus Falcon-Barroso for a very clear and useful bug report!reg_ord keyword to allow for both first and second order regularization.dim > 1.trig keyword to use a trigonometric series as alternative to Legendre polynomials.next_fast_len() for optimal rfft() zero padding.gas_component in the .plot() method, to distinguish gas emission lines in best-fitting plots.linear_fit() and nonlinear_fit() functions to better clarify the code structure. Included templates_rfft keyword.linear keyword. Added .status attribute. Changes suggested by Kyle B. Westfall (Santa Cruz).linear keyword to only perform a linear fit and skip the non-linear optimization.Chi**2/DOF instead of Biweight estimate.moments to be an arbitrary integer.moments with multiple kinematic components.sigma < velscale.oversample keyword, which is now unnecessary.start to be smaller than moments. After feedback by Masato Onodera (NAOJ).bounds and fixed.velscale_ratio keyword to pass a set of templates with higher resolution than the galaxy spectrum.oversample keyword to require integers not Booleans.bounds, fixed and fraction keywords.moments.oversample.reddening keyword. Thanks to Masatao Onodera for reporting the problem.matrix in output. Modified plotting colours.