Line Profiles & MCMC Fitting

Analytical line profile functions and the MCMC fitting routines used by the collapse_analytical family of methods.

Line Profiles

bettermoments.profiles.gaussian(x, *params)

Gaussian function with Doppler width.

Parameters:

  • x (arr) – Velocity axis in [m/s].

  • params (tuple) – The line center in [m/s], the line Doppler width in [m/s] and the line peak in [Jy/beam].

Returns:

Model spectrum in [Jy/beam].

Return type:

model (arr)

bettermoments.profiles.gaussthick(x, *params)

Gaussian profile with non-negligible optical depth,

I(v) = I_{\nu} \big(1 - \exp(\mathcal{G}(v, v0, \Delta V, \tau))\big) \, / \, (1 - \exp(-\tau))

where \mathcal{G} is a Gaussian function. Note that \tau is forced to be non-negative, but negative values will be clipped to 0.

Parameters:

  • x (arr) – Velocity axis in [m/s].

  • params (tuple) – The line center in [m/s], the line Doppler width in [m/s], the line peak in [Jy/beam] and the optical depth.

Returns:

Model spectrum in [Jy/beam].

Return type:

model (arr)

bettermoments.profiles.gausshermite(x, *params)

Gauss-Hermite expansion with the Doppler width. This allows for a flexible line profile that purely a Gaussian, where the h3 and h4 terms quantify the skewness and kurtosis of the line as in van der Marel & Franx (1993).

Parameters:

  • x (arr) – Velocity axis in [m/s].

  • params (tuple) – The line center in [m/s], the Doppler width in [m/s], the line peak in [Jy/beam], the h3 skewness coefficient and the h4 kurtosis coefficient.

Returns:

Model spectrum in [Jy/beam].

Return type:

model (arr)

bettermoments.profiles.doublegauss(x, *params)

Two gaussian components with individual widths.

Parameters:

  • x (arr) – Velocity axis in [m/s].

  • params (tuple) – The line center in [m/s], the line Doppler width in [m/s] and the line peak in [Jy/beam] for each of the two components.

Returns:

Model spectrum in [Jy/beam].

Return type:

model (arr)

bettermoments.profiles.build_cube(x, moments, method)

From a list of (N, M) moment maps, construct a model data cube on the coordinate systems of the provided cube.

Parameters:

  • x (arr) – Velocity axis in [m/s].

  • moments (array) – Array of best-fit parameters.

  • method (str) – Method used to decompose the data.

Returns:

3D data cube.

Return type:

cube (array)

bettermoments.profiles.free_params(model_function)

Return the number of free parameters for the given model function.

Parameters:

model_function (str) – Name of the model function, e.g. 'gaussian', 'gaussthick', 'gausshermite', 'doublegauss', or any of their _cont variants.

Returns:

Number of free parameters in the model.

Return type:

int

MCMC Fitting

bettermoments.mcmc_sampling.fit_cube(velax, data, rms, model_function, indices=None, ncpu=1, acf=None, **kwargs)

Fit each spectrum in indices using model_function. Spectra with fewer finite values than twice the number of free parameters are skipped.

Parallelism is handled at the per-pixel level via multiprocessing.Pool so work is evenly distributed across workers. The shared arrays (velax, data, etc.) are copied once per worker process via the pool initializer rather than being re-pickled for every task.

For more information on kwargs, see the fit_spectrum documentation.

Parameters:

  • velax (ndarray) – Velocity axis of the cube.

  • data (ndarray) – Intensity or brightness temperature array. The first axis must be the velocity axis.

  • rms (float) – Noise per pixel in same units as data.

  • model_function (str) – Name of the model function to fit to the data. Must be a function within profiles.py.

  • indices (list) – A list of pixels described by (y_idx, x_idx) tuples to fit. If none are provided, will fit all pixels.

  • ncpu (Optional[int]) – Number of worker processes. Defaults to 1 (no pool overhead).

Returns:

A (Npix, Ndim, 2) shaped array of the fits and

associated uncertainties. The uncertainties will be interleaved with the best-fit values.

Return type:

fits (ndarray)

bettermoments.mcmc_sampling.fit_spectrum(x, y, dy, model_function, p0=None, priors=None, nwalkers=None, nburnin=500, nsteps=500, mcmc='emcee', scatter=0.001, niter=1, returns='default', plots=False, cov=None, cov_chol=None, **kwargs)

Fit the provided spectrum with model_function. If mcmc is not specified, the results of the scipy.optimize.curve_fit optimization will be returned instead. Using plots=True is only recommended for debugging and when this function is not called as part of fit_cube.

Parameters:

  • x (array) – Velocity axis.

  • y (array) – Intensity axis.

  • dy (array) – Uncertainties on the intensity.

  • model_function (str) – Name of the model to fit to the spectrum. Must be a function defined in profiles.py.

  • p0 (Optional[array]) – An array of starting positions.

  • priors (Optioinal[list]) – User-defined priors.

  • nwalkers (Optional[int]) – Number of walkers for the MCMC.

  • nburnin (Optional[int]) – Number of steps to discard as burnin.

  • nsteps (Optional[int]) – Number of steps to take beyond burnin to sample the posterior distribution.

  • mcmc (Optional[str/None]) – The MCMC package to import EnsembleSampler from: 'emcee' or 'zeus'. If None, will skip the MCMC sampling and return the scipy.optimize.curve_fit results.

  • scatter (Optional[float]) – Scatter to apply to p0 values for the walkers.

  • niter (Optional[int]) – Number of MCMC iterations to run, each time using the median of the posterior samples as the new p0. Between each iteration scipy.optimize.curve_fit is not called.

  • returns (Optional[str]) – What the function returns. 'default' will return (mu, sig) for each parameter, 'percentiles' will return the 16th, 50th and 84th percentiles for each marginalized posterior distribution, 'samples' will return all posterior samples, while 'sampler' will return the EnsembleSampler.

  • plots (Optioanl[bool]) – If True, make diagnost plots.

  • free_params (Optional[int]) – The number of free parameters expected.

Returns:

Various depending on the value of returns.