# Bumps documentation

## state - Sampling history for MCMC

`MCMCDraw` |
`load_state` | filename is the path to the saved MCMC state up to the final -chain.mc, etc.
`save_state` |

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
    logp()            returns draws, logp
    acceptance_rate() returns draws, AR
    chains()          returns draws, chains, logp
    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.MCMCDraw(Ngen: int, Nthin: int, Nupdate: int, Nvar: int, Npop: int, Ncr: int, thinning: int, portion: float = 1.0)

  Bases: `object`

  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()

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

    For example, to plot the acceptance rate over time:

        draw, AR = state.acceptance_rate()
        plot(draw, 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, selection: Dict[int | Literal['logp'], Tuple[float, float]] | None = None, thin: int = 1)

    Return a sample from the posterior distribution.

    portion is the portion of each chain to use

    vars is a list of variables to return for each point

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

    thin takes every nth item.

    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 | Literal['logp'], 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()

    Compute the Gelman and Rubin R-statistic for the Markov 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
  logp(full: 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:

        draw, logp = state.logp()
        plot(draw, logp)

    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.

  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)
  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
  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.load_state(filename, skip=0, report=0, derived_vars=0)

  filename is the path to the saved MCMC state up to the final -chain.mc, etc.

  derived_vars is the number of columns added to each point, derived from other columns in that point. The newer set_derived_vars interface generates the derived variables on demand rather than storing them in the state object and so it will be zero always.

bumps.dream.state.save_state(state: MCMCDraw, filename: str)