state - Sampling history for MCMC

MCMCDraw	Initialization parameters:
Draw	A set of sample points from an MCMC draw with restrictions applied.
dream_load	Load the saved DREAM state.
h5load
h5dump

Sampling history for MCMC.

MCMC keeps track of a number of things during sampling.

The results may be queried as follows:

draws, generation, thinning, total_generations
sample(condition) returns draws, points, logp
gen_logp()        returns genid, logp
acceptance_rate() returns draws, AR
traces()          returns genid, chains
CR_weight()       returns draws, CR_weight
best()            returns best_x, best_logp
outliers()        returns outliers
show()/save(file)/load(file)

Data is stored in circular arrays, which keeps the last N generations and throws the rest away.

draws is the total number of draws from the sampler.

generation is the total number of generations. Due to tests for partially full circular buffers throughout the code (state.generation < state.Ngen) we are resetting generation to the size of the stored history on resume, and setting the generation_offset to the start of the history. If you need the number of generations across resume then use total_generations.

thinning is the number of generations per stored sample.

draws[i] is the number of draws including those required to produce the information in the corresponding return vector. Note that draw numbers need not be linearly spaced, since techniques like delayed rejection will result in a varying number of samples per generation.

logp[i] is the set of log likelihoods, one for each member of the population. The logp() method returns the complete set, and the sample() method returns a thinned set, with on element of logp[i] for each vector point[i, :].

AR[i] is the acceptance rate at generation i, showing the proportion of proposed points which are accepted into the population.

chains[i, :, :] is the set of points in the differential evolution population at thinned generation i. Ideally, the thinning rate of the MCMC process is chosen so that thinned generations i and i+1 are independent samples from the posterior distribution, though there is a chance that this may not be the case, and indeed, some points in generation i+1 may be identical to those in generation i. Actual generation number is i*thinning.

points[i, :] is the ith point in a returned sample. The i is just a place holder; there is no inherent ordering to the sample once they have been extracted from the chains. Note that the sample may be from a marginal distribution.

R[i] is the Gelman R statistic measuring convergence of the Markov chain.

CR_weight[i] is the set of weights used for selecting between the crossover ratios available to the candidate generation process of differential evolution. These will be fixed early in the sampling, even when adaptive differential evolution is selected.

outliers[i] is a vector containing the thinned generation number at which an outlier chain was removed, the id of the chain that was removed and the id of the chain that replaced it. We leave it to the reader to decide if the cloned samples, point[:generation, :, removed_id], should be included in further analysis.

best_logp is the highest log likelihood observed during the analysis and best_x is the corresponding point at which it was observed.

generation is the last generation number

class bumps.dream.state.Draw(state: MCMCDraw, portion: float = 1.0, vars: List[str | int] | None = None, exclude: List[str | int] | None = None, selection: Dict[int | str, Tuple[float, float]] | None = None, thin: int = 1, outliers: bool = False, derived: Callable | None = None)

Bases: object

A set of sample points from an MCMC draw with restrictions applied.

Inputs

state is the MCMCDraw containing all points.

portion is the fraction of points to use for the sample, starting from the ends of the Markov chains.

vars is the list of vars that are active. Variables can be 0-origin integers or strings matching a variable label. Match uses standard glob rules, with * matching any sequence, ? matching any character, [abc] matching one of abc, and [!abc] not matching one of abc.

exclude all labels that match a pattern in the list. Matching rules are the same as for vars. If a variable matches both vars and exclude then it is included.

selection restricts the points returned to a specific region of the parameter space using {var: (low, high)}. The var can be a 0-origin integer or label as for vars. It can also be logp if restricting by log likelihood. If a label matches multiple patterns then choose the first restriction (this may change to all that match in future).

thin is the amount of thinning to apply, using every nth point in the chain. This reduces autocorrelation in the chains and reduces artificial spikes in the marginal distributions. (The distributions will still be spiky, but the spikiness will be appropriate for the number of points in the bins.)

Attributes

points[n,k] are the set of sampled points, where n is the number of points and k is the number of parameters. The earlier points come from earlier generations, but this is an accident of the implementation and should not be relied on. n is usually the number of chains times the number of saved generations, but this will change depending on portion, selection, thin and outliers parameters. k is usually the number of fitted parameters, but the derived and active parameter lists will adjust this.

Nvar is the number of dimensions per point in the draw after including derived parameters and excluding nuisance parameters.

weights are the weights associated with each point. Since the chains are run at a single temperature the weights will all be 1.

logp[n] is the negative log likelihuud for each sampled point.

labels are the labels for the variables remaining after applying derived variables and restricting to the requested vars list.

integers is an array of flags, one per parameter* indicating whether the parameter is float (0) or int (1=floor). The values in the chain are managed as floats, and it is the responsibility of the likelihood function to transform them to integers. Models which use round, trunc or ceil rather than float are not yet supported.[1]

state is the original MCMCDraw.

[1] Note that discrete parameter handling is likely to change in the future. The derivative based optimizers in particular will need to know that the minimum step size for discrete parameters is one when computing the partial derivative. If bumps converts the parameters to integers before calling the likelihood function then we will only need a discrete flag since we will always use the floor() function.

Nvar: int

get_argsort_indices(var: int)

Sort var by value. Unlike argsort this caches the results for reuse.

integers: ndarray[tuple[Any, ...], dtype[_ScalarT]]

labels: List[str]

logp: ndarray[tuple[Any, ...], dtype[_ScalarT]]

points: ndarray[tuple[Any, ...], dtype[_ScalarT]]

portion: float

selection: Dict[int | str, Tuple[float, float]] | None

state: MCMCDraw

thin: int

title: str | None

vars: List[str | int] | None

weights: ndarray[tuple[Any, ...], dtype[_ScalarT]]

class bumps.dream.state.MCMCDraw(Ngen: int, Nthin: int, Nupdate: int, Nvar: int, Npop: int, Ncr: int, thinning: int)

Bases: object

Initialization parameters:

Ngen is the number of generations to store. These are retrievable through gen_logp() and acceptance_rate() methods. Note that this is not the same as the length of the saved chains.

Nthin is the number of thinned generations to store. These are retrievable through the draw() and traces() methods.

Nupdate is the number of crossover ratio sets to save. The DREAM algorithm runs in batches of steps, with various checks such as convergence and crossover ratio updates after each batch. The results are retrievable through the CR_weight() method.

Nvar is the number of parameters stored in each point.

Npop is the number of running chains. Note that unlike the –pop command line option, this is not the number of chains per parameter.

Ncr is the number of crossover ratios. DREAM maintains a fixed set of crossover ratios, choosing amongst them at each DE step. Depending on the success rate of the steps, it will adapt the weights after every batch of updates.

thinning is the number of update generations to skip between saves to the MC chains. Additional thinning can be done when drawing samples from the saved points.

Attributes and properties:

Ngen, Nthin, Nupdate, Nvar, Npop, Ncr gives the size of the buffers stored in state. If the buffers are not yet full, the returned sizes for the traces(), gen_logp(), acceptance_rate(), CR_weight() and draw() methods will be smaller.

generation is the number of generations in the current fit run. This is reset to zero each time the fit is resumed.

draws is the number of draws in the current fit run. This is reset to zero each time the fit is resumed. Use total_generations time Npop if you want the number of draws including all previous runs.

total_generations is the total number of generations for the fit across all fit runs.

title a title for the plot

labels a list of names for each parameter

thinning amount of thinning between the gen_logp buffer and the traces buffer. DREAM stores every nth generation starting at generation n (1-origin). It is possible to change thinning on resume, which will change the short range correlation in the traces, however these changes are not recorded. Generation number on the trace plot will be incorrect until DREAM burns through the entire buffer.

portion gives the fraction of each chain to use for statistical plots. A value can be automatically determined by calling state.portion = state.trim_portion(), or supplied by the user to various functions that work with the chains. The value should be between 0.0 and 1.0, where 1.0 means the entire chain is used for statistical plots.

CR_weight()

Return the crossover ratio weights to be used in the next generation.

For example, to see if the adaptive CR is stable use:

draw, weight = state.CR_weight()
plot(draw, weight)

See crossover for details.

property Ncr

property Ngen

property Npop

property Nsamples

property Nthin

property Nupdate

property Nvar

Number of parameters in the fit

acceptance_rate(portion: float | None = None)

Return the iteration number and the acceptance rate for that iteration.

For example, to plot the acceptance rate over time:

genid, AR = state.acceptance_rate()
plot(genid, AR)

V1.0.2: Now returns (genid, AR) rather than (num_draws, AR)

best()

Return the best point seen and its log likelihood.

chains()

Returns the observed Markov chains and the corresponding likelihoods.

The return value is a tuple (draws, chains, logp).

draws is the number of samples taken up to and including the samples for the current generation.

chains is a three dimensional array of generations X chains X vars giving the set of points observed for each chain in every generation. Only the thinned samples are returned.

logp is a two dimensional array of generation X population giving the log likelihood of observing the set of variable values given in chains.

derive_vars(fn: Callable[[ndarray[tuple[Any, ...], dtype[_ScalarT]]], ndarray[tuple[Any, ...], dtype[_ScalarT]]], labels: List[str] | None = None)

* DEPRECATED *

Like set_derived_vars but operating in place, modifying the points in the history.

draw(portion: float | None = None, vars: List[int] | None = None, exclude: List[int] | None = None, selection: Dict[int | str, Tuple[float, float]] | None = None, thin: int = 1, outliers: bool = False, derived: Callable | None = None)

Return a sample from the posterior distribution.

portion is the portion of each chain to use

vars is the list of vars that are active. Variables can be 0-origin integers or strings matching a variable label. Match uses standard glob rules, with * matching any sequence, ? matching any character, [abc] matching one of abc, and [!abc] not matching one of abc.

selection sets the range each parameter in the returned distribution, using {var: (low, high)}. If var matches multiple labels, then use the first restriction only. Missing variables use the full range.

thin takes every nth item.

outliers is True if outlier chains should be included (default False).

To plot the distribution for parameter p1:

draw = state.draw()
hist(draw.points[:, 0])

To plot the interdependence of p1 and p2:

draw = state.sample()
plot(draw.points[:, 0], draw.points[:, 1], '.')

entropy(vars: List[int] | None = None, portion: float | None = None, selection: Dict[int | str, Tuple[float, float]] | None = None, n_est: int = 10000, thin: int | None = None, method: str | None = None)

Return entropy estimate and uncertainty from an MCMC draw.

portion is the portion of each chain to use (uses self.portion if None).

vars is the set of variables to marginalize over. It is None for the visible variables, or a list of variables.

vars is the list of variables to use for marginalization.

selection sets the range each parameter in the returned distribution, using {variable: (low, high)}. Missing variables use the full range.

n_est is the number of points to use from the draw when estimating the entropy (default=10000).

thin is the amount of thinning to use when selecting points from the draw.

method determines which entropy calculation to use:

• gmm: fit sample to a gaussian mixture model (GMM) with \(5 \sqrt{d}\) components where \(d\) is the number fitted parameters and estimate entropy by sampling from the GMM.

• llf: estimates likelihood scale factor from ratio of density estimate to model likelihood, then computes Monte Carlo entropy from sample; this does not work for marginal likelihood estimates. DOI:10.1109/CCA.2010.5611198

• mvn: fit sample to a multi-variate Gaussian and return the entropy of the best fit gaussian; uses bootstrap to estimate uncertainty.

• wnn: estimate entropy from nearest-neighbor distances in sample. DOI:10.1214/18-AOS1688

gelman(portion: float | None = None)

Compute the Gelman and Rubin R-statistic for the Markov chains.

gen_logp(portion: float | None = None, outliers: bool = False)

Return the iteration number and the log likelihood for each point in the individual sequences in that iteration.

For example, to plot the convergence of each sequence:

genid, logp = state.gen_logp()
plot(genid, logp)

portion is the amount to trim from the head of the chain, or None to use the stored burn point.

If outliers is True, then return all chains, not just good chains.

keep_best()

Place the best point at the end of the last good chain.

Good chains are defined by mark_outliers.

Because the Markov chain is designed to wander the parameter space, the best individual seen during the random walk may have been observed during the burn-in period, and may no longer be present in the chain. If this is the case, replace the final point with the best, otherwise swap the positions of the final and the best.

property labels

Labels associated with each parameter in a point. This doesn’t include the labels for derived quantities. For these you will need to query state.draw().labels.

logp(full: bool = False)

Return the cumulative number of draws and the log likelihoods for each chain.

Note that draw[i] represents the total number of samples taken, including those for the samples in logp[i].

If full is True, then return all chains, not just good chains.

** Deprecated ** Use state.gen_logp() instead.

logp_slice(n: int)

Return a slice of the logp chains, either the first n if n > 0 or the last n if n < 0. Avoids unrolling the circular buffer if possible.

mark_outliers(test: str = 'IQR', portion: float | None = None)

Mark some chains as outliers but don’t remove them. This can happen after drawing is complete, so that chains that did not converge are not included in the statistics.

test is ‘IQR’, ‘mahal’, ‘grubbs’, or ‘none’.

portion indicates what portion of the samples should be included in the outlier test. If None, then the stored portion is used.

min_slice(n: int)

Return the minimum logp for n slices, from the head if positive or the tail if negative.

This is a specialized function so it can be fast. Convergence can be quickly rejected if the min in a short head is smaller than the min in a long tail. Unfortunately, if the data is wrapped, then the max function will cost extra.

outliers()

Return a list of outlier removal operations.

Each outlier operation is a tuple giving the thinned generation in which it occurred, the old chain id and the new chain id.

The chains themselves have already been updated to reflect the removal.

Curiously, it is possible for the maximum likelihood seen so far to be removed by this operation.

remove_outliers(x: ndarray[tuple[Any, ...], dtype[_ScalarT]], logp: ndarray[tuple[Any, ...], dtype[_ScalarT]], test: str = 'IQR')

Replace outlier chains with clones of good ones. This should happen early in the sampling processes so the clones have an opportunity to evolve their own identity. Only the head of the chain is modified.

state contains the chains, with log likelihood for each point.

x, logp are the current population and the corresponding log likelihoods; these are updated with cloned chain values.

test is the name of the test to use (one of IQR, Grubbs, Mahal or none). See outliers.identify_outliers() for details.

Updates state, x and logp to reflect the changes.

Returns a list of the outliers that were removed.

resize(Ngen: int, Nthin: int, Nupdate: int, Nvar: int, Npop: int, Ncr: int, thinning: int)

sample(**kw)

Return a sample from the posterior distribution.

Deprecated use draw() instead.

save(filename: Path | str)

set_derived_vars(fn: Callable[[ndarray[tuple[Any, ...], dtype[_ScalarT]]], ndarray[tuple[Any, ...], dtype[_ScalarT]]], labels: List[str])

Define derived variables from the sample. When calling draw() it will add columns for the derived variables to each sample.

fn is a function taking points p[:, k] for k in 0 … samples and returning a set of derived variables pj[k] for each sample k. The variables can be returned as any kind of sequence including an array or a tuple with one entry per variable. The caller uses asarray to convert the returned variables into a vars X samples array. For convenience, a single variable can be returned by itself.

labels are the labels to use for the derived variables.

The following example adds the new variable x+y = P[0] + P[1]:

state.derive_vars(lambda p: p[0]+p[1], labels=["x+y"])

set_integer_vars(labels: List[str])

Indicate tha variables should be considered integer variables when computing statistics.

set_visible_vars(labels: List[str])

show(portion: float | None = None, figfile: str | Path | None = None)

show_labels()

List of parameters in the state, with the parameter number for each. Use this to help with the inputs to state.draw().

stable_best()

Return True if the at least one full cycle of the circular buffer has passed since the best logp was first observed.

title = None

property total_generations

traces(portion: float | None = None, thin: int = 1, outliers: bool = False)

Grab traces for trace plot.

portion gives the starting point of the trace, skipping over the initial frames that may not have converged.

thin is the amount of additional thinning to apply on the traces.

outliers (default False) is True if bad chains should be included in the trace.

Returns generation [Ngen] and chains [Ngen x Nchain x Nvar]

trim_index(portion: float | None = None, generation: int | None = None)

Returns the generation index corresponding to the trim portion. The returned generation is relative to the start of the sampling, even if resume has been called multiple times to extend the burn or the sample size.

The optional generation parameter is needed in case we have a state object that is out of date with respect to the number of generations seen so far. This can happen when the state is transferred to the user interface thread much less frequently than the step monitor.

trim_portion()

Estimate the point at which the Markov chains have reached equilibrium, returning it as a portion of the currently stored length.

If not converged, then trim the first half of the samples.

bumps.dream.state.dream_load(store)

Load the saved DREAM state.

store is the path to the stored state, either as an HDF5 history file or as a bumps export directory. If the directory contains multiple exports then use the path to the .par file within the directory.

See also: h5load, load_state

bumps.dream.state.h5dump(group: Group, state: MCMCDraw)

bumps.dream.state.h5load(group: Group)