ModifiedBlackbody (MBB)

ModifiedBlackbody (MBB)#

Modified blackbody model for SED fitting. See Modeling options for possible models that can be used.

class mbb.ModifiedBlackbody(L, T, beta, z, alpha=2.0, l0=200.0, opthin=True, pl=False, pl_piecewise=False)#

Class to represent a modified blackbody (or MBB).

This class can be used to encapsulate a single MBB model, or to perform an SED fit to photometry. The results can be easily plotted or updated as needed, and various parameters/statistics can be extracted. The models are based off of Casey (2012).

Parameters:

  • L (float) – log10 of luminosity in solar units. If fitting data, this will set the initial guess for the fit.

  • T (float) – dust temperature in K. If fitting data, this will set the initial guess for the fit.

  • beta (float) – dust emissivity spectral index. If fitting data, this will set the initial guess for the fit.

  • z (float) – Redshift of this galaxy.

  • alpha (float, optional) – mid-IR power-law slope.

  • l0 (float,optional) – wavelength where opacity equals unity, in microns.

  • opthin (bool) – Whether or not the model should assume optically thin dust emission.

  • pl (pool) – Whether or not the model should include a MIR power law (as in Casey+ 2012)

  • pl_piecewise (bool) – if the powerlaw should be attached piecewise (as in Drew&Casey 2022) or with exponential fall-off (as in Casey+ 2012). Currently, if set to True, the powerlaw is joined at the wavelength where the blackbody slope equals alpha, otherwise it smoothly turns over at 3/4 times this wavelength.

Note: By default, ModifiedBlackbody assumes a flat \(\Lambda\mathrm{CDM}\) cosmology with \(\Omega_m = 0.3\) and \(\Omega_\Lambda = 0.7\). If you wish to change this, the code allows you to set the cosmo attribute of the ModifiedBlackbody to an instance of astropy.cosmology.Cosmology after it is created.

eval(wl, z=None)#

Evaluate MBB at a given wavelength

Return evaluation of this MBB’s model if observed at the given wavelengths wl shifted to redshift z, in Jy. Set z=0 to get rest-frame evaluation. Default is to give observed-frame flux.

Note that this function is not luminosity-preserving…i.e., it does not change the underlying MBB parameters and will simply shift the SED to the redshift without applying any dimming. To get a constant luminosity at different redshifts, use MBB.update(z=new_z) to update the MBB parameters to be consistent with the new redshift, and then use MBB.eval(wl) to get the SED.

Parameters:

  • wl (float) – wavelength(s) in micron

  • z (float) – redshift at which the model should be evaluated.

Returns:

value of mbb at the wavelength wl

Return type:

float

fit(phot, uplims=None, params=['L', 'T', 'beta'], priors=None, restframe=False, nwalkers=400, nburn=1000, niter=1000, ncores=0, stepsize=1e-07, pool=None)#

Fit model to photometry

Fit a modified blackbody to photometry. Updates the parameters of this MBB model to the resulting parameters of the fit, and populates the fit_result attribute of the ModifiedBlackbody with the fit results. The parameter values of the MBB model will be set to the medians of the posterior distributions. Percentiles of the posterior (e.g., to measure the 68% confidence interval) can be obtained with the post_percentile() method.

The fit_result attribute is a dictionary containing the following:

  • sampler: an emcee.EnsembleSampler representing the chain of walker values from the fit. The parameter values are ordered in the same order as the params argument passed to fit().

  • chi2: the raw chi-squared value at the end of the fitting process.

  • n_params: the number of fitted parameters.

  • n_bands: the number of bands in the fit (the length of phot)

Example:

from mbb import ModifiedBlackbody as MBB
m = MBB(L=12, beta=1.8, T=35, z=2.5, alpha=2.0, opthin=True, pl = True)
phot = ([450, 850],[0.005, 0.0021],[0.0006,0.00032]) #wl, flux, error

# fit for redshift
result = m.fit(phot=phot,niter=100,params=['L','z'],
               restframe=False, priors = {'z':dict(mu=2.5,sigma=1.0)})
# could equivalently use 'm.fit_result' instead of 'result'

reduc_chi2 = result['chi2']/(result['n_bands']-result['n_params']) # reduced chi-squared
Parameters:

  • phot (array-like, 3 x N) – wavelengths and photometry, arranged as a 3 x N array (wavelength, flux, error).

  • True. (Wavelengths are treated as observed-frame unless restframe is set to)

  • uplims (bool, array-like of length N) – whether each band in phot is a upper limit (limit equal to flux if > error, otherwise limit = error)

  • nwalkers (int) – how many walkers should be used in the MCMC fit.

  • nburn (int) – number of MCMC steps to discard (default is 1000)

  • niter (int) – how many iterations to run in the fit.

  • ncores (int) – how many CPU cores to use in multiprocessing. If pool is not None, this is ignored. Set to 1 to not use multiprocessing. Default is number of available CPUs - 2.

  • stepsize (float) – stepsize used to randomize the initial walker values.

  • params (list) – list of parameter names, e.g., [L, beta, T, z, alpha, l0] to vary in the fit. The rest will be fixed. L should generally be included since it represents the normalization of the model.

  • priors (dict) – Priors to use in the Bayesian fitting. This should be a dictionary with keys corresponding to the elements of params. For each key, the corresponding value can be either (1) a function that takes the value of the parameter and returns the relative probability density (in linear space), or (2) a dictionary with keys ‘mu’ and ‘sigma’, in which case the code will use these to generate a Gaussian prior on that parameter. If None, or for any params not included in priors, flat (uniform) priors will be assumed. Currently, T is constrained to be between 5 K and 120 K, beta is between 0.1 and 5.0, and z is between 0.1 and 12.

  • restframe (bool) – whether wavelengths in phot are given in the rest frame (default is observed frame)

  • pool (multiprocessing.pool.Pool) – an optional pool to pass to the sampler for multiprocessing; otherwise fit() will generate one internally. Can be faster to use an external Pool if many fits are being performed, also may be useful depending on the OS/kernel being used.

Returns:

the fit results, which are also stored in the fit_result attribute of this ModifiedBlackbody. See above for details on the contents of this dictionary.

Return type:

dict

classmethod from_file(file)#

load in saved ModifiedBlackbody from file

get_luminosity(wllimits=(8, 1000))#

get integrated luminosity for the current MBB state between wavelength limits in microns.

Parameters:

wllimits (tuple) – rest-frame wavelength limits in microns (lo, hi) between which to integrate

Returns:

the luminosity integrated between rest-frame wavelength limits given by wllimits

Return type:

float

get_peak_wavelength()#

Get the peak (rest-frame) wavelength of this ModifiedBlackbody, in microns.

property phot#

Photometry (observed frame) used for fitting. Note that wavelengths are stored as rest-frame internally, so a conversion is applied when returning these values in observed frame.

property phot_rf#

Photometry (rest frame) used for fitting.

plot_corner(**kwargs)#

Plot a corner plot showing the results from the MCMC fit to the data. All kwargs will be passed on to corner.corner().

Returns:

the current figure

Return type:

matplotlib.pyplot.figure

plot_sed(obs_frame=False, ax=None)#

plot the rest-frame form of this mbb just for basic visualization. It is recommended to use a separate, more detailed plotting function for figures.

Parameters:

  • obs_frame (bool) – whether to plot against observed-frame wavelengths (default is rest frame).

  • ax (matplotlib.pyplot.Axes) – axes to plot the model on.

Returns:

the current figure matplotlib.axes.Axes: the current axes

Return type:

matplotlib.pyplot.figure

post_percentile(param, q=[16, 50, 84], sample_by=10)#

Determine the posterior percentile values of a given fit parameter.

Example:

m.post_percentile('beta', q = [16,50,84]) # get median and 16th--84th percentile interval
Parameters:

  • param (str) – name of parameter, element of the ‘params’ argument passed to ModifiedBlackbody.fit()

  • q (array-like of float) – percentile or sequence of percentiles to compute from the posterior distribution

  • sample_by (int) – sample every sample_by values in the posterior chain, helps with speed (especially for L, which must be computed).

Returns:

the percentiles of the posterior distribution for parameter ‘param’

Return type:

float or array

posterior(param, sample_by=10)#

Determine the posterior chain of a given fit parameter.

Parameters:

  • param (str) – name of parameter, element of the ‘params’ argument passed to ModifiedBlackbody.fit()

  • sample_by (int) – sample every sample_by values in the posterior chain, helps with speed (especially for L, which must be computed).

Returns:

the chain of values in the posterior distribution for parameter ‘param’

Return type:

array

reset()#

Clear the fit results.

Clear the ModifiedBlackbody fit results, priors, and photometry. Current values of parameters (L, beta, etc) will remain unchanged.

save(file, protocol=None)#

save out ModifiedBlackbody (using pickle) to file in binary format, including fit and sampler

update(L=None, T=None, beta=None, z=None, alpha=None, l0=None)#

update specified modified blackbody parameters (not the underlying model), keeping others constant.