emcee

SeparableModelResult.emcee(params=None, steps=1000, nwalkers=100, burn=0, thin=1, ntemps=1, pos=None, reuse_sampler=False, workers=1, float_behavior='posterior', is_weighted=True, seed=None, progress=True)

Bayesian sampling of the posterior distribution using emcee.

Bayesian sampling of the posterior distribution for the parameters using the emcee Markov Chain Monte Carlo package. The method assumes that the prior is Uniform. You need to have emcee installed to use this method.

Parameters:

• params (Parameters, optional) – Parameters to use as starting point. If this is not specified then the Parameters used to initialize the Minimizer object are used.
• steps (int, optional) – How many samples you would like to draw from the posterior distribution for each of the walkers?
• nwalkers (int, optional) – Should be set so \(nwalkers >> nvarys\), where nvarys are the number of parameters being varied during the fit. “Walkers are the members of the ensemble. They are almost like separate Metropolis-Hastings chains but, of course, the proposal distribution for a given walker depends on the positions of all the other walkers in the ensemble.” - from the emcee webpage.
• burn (int, optional) – Discard this many samples from the start of the sampling regime.
• thin (int, optional) – Only accept 1 in every thin samples.
• ntemps (int, optional) – If ntemps > 1 perform a Parallel Tempering.
• pos (numpy.ndarray, optional) – Specify the initial positions for the sampler. If ntemps == 1 then pos.shape should be (nwalkers, nvarys). Otherwise, (ntemps, nwalkers, nvarys). You can also initialise using a previous chain that had the same ntemps, nwalkers and nvarys. Note that nvarys may be one larger than you expect it to be if your userfcn returns an array and is_weighted is False.
• reuse_sampler (bool, optional) – If you have already run emcee on a given Minimizer object then it possesses an internal sampler attribute. You can continue to draw from the same sampler (retaining the chain history) if you set this option to True. Otherwise a new sampler is created. The nwalkers, ntemps, pos, and params keywords are ignored with this option. Important: the Parameters used to create the sampler must not change in-between calls to emcee. Alteration of Parameters would include changed min, max, vary and expr attributes. This may happen, for example, if you use an altered Parameters object and call the minimize method in-between calls to emcee.
• workers (Pool-like or int, optional) – For parallelization of sampling. It can be any Pool-like object with a map method that follows the same calling sequence as the built-in map function. If int is given as the argument, then a multiprocessing-based pool is spawned internally with the corresponding number of parallel processes. ‘mpi4py’-based parallelization and ‘joblib’-based parallelization pools can also be used here. Note: because of multiprocessing overhead it may only be worth parallelising if the objective function is expensive to calculate, or if there are a large number of objective evaluations per step (ntemps * nwalkers * nvarys).
• float_behavior (str, optional) –

  Specifies meaning of the objective function output if it returns a float. One of:

  • ’posterior’ - objective function returns a log-posterior probability
  • ’chi2’ - objective function returns \(\chi^2\)

  See Notes for further details.

• is_weighted (bool, optional) – Has your objective function been weighted by measurement uncertainties? If is_weighted is True then your objective function is assumed to return residuals that have been divided by the true measurement uncertainty (data - model) / sigma. If is_weighted is False then the objective function is assumed to return unweighted residuals, data - model. In this case emcee will employ a positive measurement uncertainty during the sampling. This measurement uncertainty will be present in the output params and output chain with the name __lnsigma. A side effect of this is that you cannot use this parameter name yourself. Important this parameter only has any effect if your objective function returns an array. If your objective function returns a float, then this parameter is ignored. See Notes for more details.
• seed (int or numpy.random.RandomState, optional) – If seed is an int, a new numpy.random.RandomState instance is used, seeded with seed. If seed is already a numpy.random.RandomState instance, then that numpy.random.RandomState instance is used. Specify seed for repeatable minimizations.

Returns:

MinimizerResult object containing updated params, statistics, etc. The updated params represent the median (50th percentile) of all the samples, whilst the parameter uncertainties are half of the difference between the 15.87 and 84.13 percentiles. The MinimizerResult also contains the chain, flatchain and lnprob attributes. The chain and flatchain attributes contain the samples and have the shape (nwalkers, (steps - burn) // thin, nvarys) or (ntemps, nwalkers, (steps - burn) // thin, nvarys), depending on whether Parallel tempering was used or not. nvarys is the number of parameters that are allowed to vary. The flatchain attribute is a pandas.DataFrame of the flattened chain, chain.reshape(-1, nvarys). To access flattened chain values for a particular parameter use result.flatchain[parname]. The lnprob attribute contains the log probability for each sample in chain. The sample with the highest probability corresponds to the maximum likelihood estimate.

Return type:

MinimizerResult

Notes

This method samples the posterior distribution of the parameters using Markov Chain Monte Carlo. To do so it needs to calculate the log-posterior probability of the model parameters, F, given the data, D, \(\ln p(F_{true} | D)\). This ‘posterior probability’ is calculated as:

\[\ln p(F_{true} | D) \propto \ln p(D | F_{true}) + \ln p(F_{true})\]

where \(\ln p(D | F_{true})\) is the ‘log-likelihood’ and \(\ln p(F_{true})\) is the ‘log-prior’. The default log-prior encodes prior information already known about the model. This method assumes that the log-prior probability is -numpy.inf (impossible) if the one of the parameters is outside its limits. The log-prior probability term is zero if all the parameters are inside their bounds (known as a uniform prior). The log-likelihood function is given by [1]:

\[\ln p(D|F_{true}) = -\frac{1}{2}\sum_n \left[\frac{(g_n(F_{true}) - D_n)^2}{s_n^2}+\ln (2\pi s_n^2)\right]\]

The first summand in the square brackets represents the residual for a given datapoint (\(g\) being the generative model, \(D_n\) the data and \(s_n\) the standard deviation, or measurement uncertainty, of the datapoint). This term represents \(\chi^2\) when summed over all data points. Ideally the objective function used to create lmfit.Minimizer should return the log-posterior probability, \(\ln p(F_{true} | D)\). However, since the in-built log-prior term is zero, the objective function can also just return the log-likelihood, unless you wish to create a non-uniform prior.

If a float value is returned by the objective function then this value is assumed by default to be the log-posterior probability, i.e. float_behavior is ‘posterior’. If your objective function returns \(\chi^2\), then you should use a value of ‘chi2’ for float_behavior. emcee will then multiply your \(\chi^2\) value by -0.5 to obtain the posterior probability.

However, the default behaviour of many objective functions is to return a vector of (possibly weighted) residuals. Therefore, if your objective function returns a vector, res, then the vector is assumed to contain the residuals. If is_weighted is True then your residuals are assumed to be correctly weighted by the standard deviation (measurement uncertainty) of the data points (res = (data - model) / sigma) and the log-likelihood (and log-posterior probability) is calculated as: -0.5 * numpy.sum(res**2). This ignores the second summand in the square brackets. Consequently, in order to calculate a fully correct log-posterior probability value your objective function should return a single value. If is_weighted is False then the data uncertainty, s_n, will be treated as a nuisance parameter and will be marginalized out. This is achieved by employing a strictly positive uncertainty (homoscedasticity) for each data point, \(s_n = \exp(\_\_lnsigma)\). __lnsigma will be present in MinimizerResult.params, as well as Minimizer.chain, nvarys will also be increased by one.

References

[1]	http://dan.iel.fm/emcee/current/user/line/