# Optimizer Selection¶

Bumps has a number of different optimizers available, each with its own control parameters:

Random Lines [experimental]

Particle Swarm [experimental]

Parallel Tempering [experimental]

In general there is a trade-off between convergence rate and robustness, with the fastest algorithms most likely to find a local minimum rather than a global minimum. The gradient descent algorithms (Levenberg-Marquardt, Quasi-Newton BFGS) tend to be fast but they will find local minima only, while the population algorithms (DREAM, Differential Evolution) are more robust and likely slower. Nelder-Mead Simplex is somewhere between, with a small population keeping the search local but more robust than the gradient descent algorithms.

Each algorithm has its own set of control parameters for adjusting the search process and the stopping conditions. The same option may mean slightly different things to different optimizers. The Bumps package provides a dialog box for selecting the optimizer and its options when running the fit wx application. This only includes the common options for the most useful optimizers. For full control, the fit will need to be run from the command line interface or through a python script.

For parameter uncertainty, most algorithms use the covariance matrix at
the optimum to estimate an uncertainty ellipse. This is okay for a
preliminary analysis, but only works reliably for weakly correlated parameters.
For full uncertainty analysis, DREAM uses a random walk to explore
the parameter space near the minimum, showing pair-wise correlations
amongst the parameter values. In order for DREAM to return the
correct uncertainty, the function to be optimized should be a conditional
probability density, with *nllf* as the negative log likelihood function
of seeing point \(x\) in the parameter space. Other functions
can be fitted, but uncertainty estimates will be meaningless.

Most algorithms have been adapted to run in parallel at least to some degree. The implementation is not heavily tuned, either in terms of minimizing the overhead per function evaluation or for distributing the problem across multiple processors. If the theory function is implemented in parallel, then the optimizer should be run in serial. Mixed mode is also possible when running on a cluster with a multi-threaded theory function. In this case, only one theory function will be evaluated on each cluster node, but the optimizer will distribute the parameters values to the cluster nodes in parallel. Do not run serial algorithms (Levenberg-Marquardt, Quasi-Newton BFGS) on a cluster.

We have included a number of optimizers in Bumps that did not perform particularly well on our problem sets. However, they may be perfect for your problem, so we have left them in the package for you to explore. They are not available in the GUI selection.

## Levenberg-Marquardt¶

The Levenberg-Marquardt [1] [2] algorithm has been the standard method for non-linear data fitting. As a gradient descent trust region method, it starts at the initial value of the function and steps in the direction of the derivative until it reaches the minimum. Set up as an explicit minimization of the sum of square differences between theory and model, it uses a numerical approximation of the Jacobian matrix to set the step direction and an adaptive algorithm to set the size of the trust region.

### When to use¶

Use this method when you have a reasonable fit near the minimum, and
you want to get the best possible value. This can then be used as the starting
point for uncertainty analysis using DREAM. This method requires
that the problem definition includes a *residuals* method, but this should
always be true when fitting data.

When modeling the results of an experiment, the best fit value is an accident of the measurement. Redo the same measurement, and the slightly different values you measure will lead to a different best fit. The important quantity to report is the credible interval covering 68% (1-\(\sigma\)) or 95% (2-\(\sigma\)) of the range of parameter values that are somewhat consistent with the data.

This method uses *lmfit* from *scipy*, and does not run in parallel.

### Options¶

*Steps* is the number of gradient steps to take. Each step requires
a calculation of the Jacobian matrix to determine the direction. This
needs \(2 m n\) function evaluations, where \(n\) is the number of parameters and
each function is evaluated and \(m\) data points (assuming center point
formula for finite difference estimate of the derivative). The resulting
linear equation is then solved, but for small \(n\) and expensive function
evaluation this overhead can be ignored. Use `--steps=n`

from
the command line.

*f(x) tolerance* and *x tolerance* are used to determine when
the fit has reached the point where no significant improvement is expected.
If the function value does not improve significantly within the step, or
the step is too short, then the fit will terminate. Use `--ftol=v`

and
`--xtol=v`

from the command line.

From the command line, `--starts=n`

will automatically restart the algorithm
after it has converged so that a slightly better value can be found. If
`--keep_best`

is included then restart will use a value near the minimum,
otherwise it will restart the fit from a random point in the parameter space.

Use `--fit=lm`

to select the Levenberg-Marquardt fitter from the command line.

### Notes¶

v0.8.2 Changed from scipy.leastsq to mpfit for better bounds handling. Use
`--fit=scipy.leastsq`

to restore the previous behaviour.

### References¶

## Nelder-Mead Simplex¶

The Nelder-Mead [3] downhill simplex algorithm is a robust optimizer which does not require the function to be continuous or differentiable.

It uses the relative values of the function at the corners of a
simplex (an n-dimensional triangle) to decide which points of the simplex
to update. It will take the worst value and try moving it inward or
outward, or reflect it through the centroid of the remaining values
stopping if it finds a better value. If none of these values are
better, then it will shrink the simplex and start again. The name
amoeba comes from the book *Numerical Recipes* [4] wherein they
describe the search as acting like an amoeba, squeezing through narrow valleys
as it makes its way down to the minimum.

### When to use¶

Use this method as a first fit to your model. If your fitting function is well behaved with few local minima this will give a quick estimate of the model, and help you decide if the model needs to be refined. If your function is poorly behaved, you will need to select a good initial value before fitting, or use a more robust method such as Differential Evolution or DREAM.

The uncertainty reported comes from a numerical derivative estimate at the minimum.

This method requires a series of function updates, and does not benefit much from running in parallel.

### Options¶

*Steps* is the simplex update iterations to perform. Most updates
require one or two function evaluations, but shrinking the simplex evaluates
every value in the simplex. Use `--steps=n`

from the command line.

*Starts* tells the optimizer to restart a given number of times.
Each time it restarts it uses a random starting point. Use
`--starts=n`

from the command line.

*Simplex radius* is the initial size of the simplex, as a portion of
the bounds defining the parameter space. If a parameter is unbounded, then
the radius will be treated as a portion of the parameter value. Use
`--radius=n`

from the command line.

*x tolerance* and *f(x) tolerance* are used to determine when the
fit has reached the point where no significant improvement is expected.
If the simplex is tiny (that is, the corners are close to each other) and
flat (that is, the values at the corners are close to each other),
then the fit will terminate. Use `--xtol=v`

and `--ftol=v`

from
the command line.

From the command line, use `--keep_best`

so that restarts are centered on a
value near the minimum rather than restarting from a random point within the
parameter bounds.

Use `--fit=amoeba`

to select the Nelder-Mead simplex fitter from the
command line.

### References¶

Nelder, J. A.; Mead, R.
*The Computer Journal*
1965, 7 (4), 308–313.
DOI: 10.1093/comjnl/7.4.308

Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T.
In *Numerical Recipes in C: The Art of Scientific Computing, Second Edition*;
Cambridge University Press: Cambridge; New York, 1992; pp 408–412.

## Quasi-Newton BFGS¶

Broyden-Fletcher-Goldfarb-Shanno [5] is a gradient descent method which uses the gradient to determine the step direction and an approximation of the Hessian matrix to estimate the curvature and guess a step size. The step is further refined with a one-dimensional search in the direction of the gradient.

### When to use¶

Like Levenberg-Marquardt, this method converges quickly to the minimum. It does
not assume that the problem is in the form of a sum of squares and does not
require a *residuals* method.

The \(n\) partial derivatives are computed in parallel.

### Options¶

*Steps* is the number of gradient steps to take. Each step requires
a calculation of the Jacobian matrix to determine the direction. This
needs \(2 m n\) function evaluations, where \(n\) is the number of parameters and
each function is evaluated and \(m\) data points (assuming center point
formula for finite difference estimate of the derivative). The resulting
linear equation is then solved, but for small \(n\) and expensive function
evaluation this overhead can be ignored.
Use `--steps=n`

from the command line.

*Starts* tells the optimizer to restart a given number of times.
Each time it restarts it uses a random starting point.
Use `--starts=n`

from the command line.

*f(x) tolerance* and *x tolerance* are used to determine when
the fit has reached the point where no significant improvement is expected.
If the function is small or the step is too short then the fit
will terminate. Use `--ftol=v`

and `--xtol=v`

from the command line.

From the command line, `--keep_best`

uses a value near the previous minimum
when restarting instead of using a random value within the parameter bounds.

Use `--fit=newton`

to select BFGS from the command line.

### References¶

Dennis, J. E.; Schnabel, R. B.
*Numerical Methods for Unconstrained Optimization and Nonlinear Equations*;
Society for Industrial and Applied Mathematics: Philadelphia, 1987.

## Differential Evolution¶

Differential evolution [6] is a population based algorithm which uses
differences between points as a guide to selecting new points. For each member
of the population a pair of points is chosen at random, and a difference vector
is computed. This vector is scaled, and a random subset of its components are
added to the current point based on crossover ratio. This new point is
evaluated, and if its value is lower than the current point, it replaces
it in the population. There are many variations available within DE that
have not been exposed in Bumps. Interested users can modify
`bumps.fitters.DEFit`

and experiment with different crossover and
mutation algorithms, and perhaps add them as command line options.

Differential evolution is a robust directed search strategy. Early in the search, when the population is disperse, the difference vectors are large and the search remains broad. As the search progresses, more of the population goes into the valleys and eventually all the points end up in local minima. Now the differences between random pairs will often be small and the search will become more localized.

The population is initialized according to the prior probability distribution for each each parameter. That is, if the parameter is bounded, it will use a uniform random number generate within the bounds. If it is unbounded, it will use a uniform value in [0,1]. If the parameter corresponds to the result of a previous measurement with mean \(\mu\) and standard deviation \(\sigma\), then the initial values will be pulled from a Gaussian random number generator.

### When to use¶

Convergence with differential evolution will be slower, but more robust.

Each update will evaluate \(k\) points in parallel, where \(k\) is the size of the population.

### Options¶

*Steps* is the number of iterations. Each step updates each member
of the population. The population size scales with the number of fitted
parameters. Use `--steps=n`

from the command line.

*Population* determines the size of the population. The number of
individuals, \(k\), is equal to the number of fitted parameters times the
population scale factor. Use `--pop=k`

from the command line.

*Crossover ratio* determines what proportion of the dimensions to update
at each step. Smaller values will likely lead to slower convergence, but
more robust results. Values must be between 0 and 1. Use `--CR=v`

from
the command line.

*Scale* determines how much to scale each difference vector before adding
it to the candidate point. The selected mutation algorithm chooses a scale
factor uniformly in \([0,F]\). Use `--F=v`

from the command line.

*f(x) tolerance* and *x tolerance* are used to determine when the
fit has reached the point where no significant improvement is expected.
If the population is flat (that is, the minimum and maximum values are
within tolerance) and tiny (that is, all the points are close to each
other) then the fit will terminate. Use `ftol=v`

and `xtol=v`

from the
command line.

Use `--fit=de`

to select differential evolution from the command line.

### References¶

Storn, R.; Price, K.
*Journal of Global Optimization*
1997, 11 (4), 341–359.
DOI: 10.1023/A:1008202821328

## DREAM¶

DREAM [7] is a population based algorithm like differential evolution, but instead of only keeping individuals which improve each generation, it will sometimes keep individuals which get worse. Although it is not fast and does not give the very best value for the function, we have found it to be a robust fitting engine which will give a good value given enough time.

The progress of each individual in the population from generation to generation can considered a Markov chain, whose transition probability is equal to the probability of taking the step times the probability that it keeps the step based on the difference in value between the points. By including a purely random stepper with some probability, the detailed balance condition is preserved, and the Markov chain converges onto the underlying equilibrium distribution. If the theory function represents the conditional probability of selecting each point in the parameter space, then the resulting chain is a random draw from the posterior distribution.

This means that the DREAM algorithm can be used to determine the parameter uncertainties. Unlike the hessian estimate at the minimum that is used to report uncertainties from the other fitters, the resulting uncertainty need not Gaussian. Indeed, the resulting distribution can even be multi-modal. Fits to measured data using theory functions that have symmetric solutions have shown all equivalent solutions with approximately equal probability.

### When to use¶

Use DREAM when you need a robust fitting algorithm. It takes longer but it does an excellent job of exploring different minima and getting close to the global optimum.

Use DREAM when you want a detailed analysis of the parameter uncertainty.

Like differential evolution, DREAM will evaluate \(k\) points in parallel, where \(k\) is the size of the population.

### Options¶

*Samples* is the number of points to be drawn from the Markov chain.
To estimate the 68% interval to two digits of precision, at least
1e5 (or 100,000) samples are needed. For the 95% interval, 1e6
(or 1,000,000) samples are needed. The default 1e4 samples
gives a rough approximation of the uncertainty relatively quickly.
Use `--samples=n`

from the command line.

*Burn-in Steps* is the number of iterations to required for the Markov
chain to converge to the equilibrium distribution. If the fit ends
early, the tail of the burn will be saved to the start of the steps.
Use `--burn=n`

from the command line.

*Population* determines the size of the population. The number of
individuals, \(k\), is equal to the number of fitted parameters times the
population scale factor. Use `--pop=k`

from the command line.

*Initializer* determines how the population will be initialized.
The options are as follows:

eps(epsilon ball), in which the entire initial population is chosen at random from within a tiny hypersphere centered about the initial point

lhs(latin hypersquare), which chops the bounds within each dimension in \(k\) equal sized chunks where \(k\) is the size of the population and makes sure that each parameter has at least one value within each chunk across the population.

cov(covariance matrix), in which the uncertainty is estimated using the covariance matrix at the initial point, and points are selected at random from the corresponding Gaussian ellipsoid

random(uniform random), in which the points are selected at random within the bounds of the parameters

Use `--init=type`

from the command line.

*Thinning* is the amount of thinning to use when collecting the
population. If the fit is somewhat stuck, with most steps not improving
the fit, then you will need to thin the population to get proper
statistics. Use `--thin=k`

from the command line.

*Convergence* gives a cutoff value \(\alpha\) for determining when
the Markov chain has converged. The default is `--alpha=0.00`

for no
convergence tests. Various tests are used, such as comparing the distribution
of points in the first part of the chain to the last part and looking for
trends in the log-likelihood values. You may need to use smaller \(\alpha\) for
shorter sequences (samples over variables times population) since the test
statistics will have higher variance. Convergence is tested every \(n\) steps.

*Outliers* is the test to use to check for outlier chains. Default is
`--outliers=none`

for no outlier test. Options are *iqr*, which uses
the inter-quartile range on the likelihoods, *grubbs*, which uses a t-test
on the likelihoods, and *mahal* which looks at the distance from the best
chain in parameter space. Outlier removal occurs every \(2n\) steps where
\(n\) is #samples/(#pars #pop), or when the convergence test indicates the
chains are stable. Outliers are replaced by non-outlier chains at random.
These new chains need at least \(n\) steps to mix before being used. If the
MCMC exploration stops due to time, some of the chains may not be properly
mixed.

*Burn-in trim* is used to clear spurious samples from the Markov chains.
If `--trim=true`

then Bumps finds the “burn point” after which the
chains appear to have converged. Samples before this point are ignored
when computed statistics and making plots. The trimmed samples are still
written to the MCMC output files so they will be available when the fit
is resumed.

*Calculate entropy*, if true, computes the entropy for the fit. This is
an estimate of the amount of information in the data. Use `--entropy=method`

from the command line, where method is one of *llf* (default), *gmm*, *mvn*
or *wnn*. See below for details.

*Steps*, if not zero, determines the number of iterations to use for
drawing samples after burn in. Each iteration updates the full population,
which is (population x number of fitted parameters) points. This option
is available for compatibility; it is more useful to set the number of
samples directly. Use `--steps=n`

from the command line.

Use `--fit=dream`

to select DREAM from the command line. Consider using
`--parallel`

and `--checkpoint`

as well. When running in a batch queue,
add `--batch`

and use `--mpi`

rather than `--parallel`

.

### Output¶

DREAM produces a number of different outputs, and there are a number of
things to check before using its reported uncertainty values. The main
goal of selecting `--burn=n`

is to wait long enough to reach the
equilibrium distribution.

For each parameter in the fit, DREAM finds the mean, median and best value, as well as the 68% and 95% credible intervals. The mean value is defined as \(\int x P(x) dx\), which is just the expected value of the probability distribution for the parameter. The median value is the 50% point in the probability distribution, and the best value is the maximum likelihood value seen in the random walk. The credible intervals are the central intervals which capture 68% and 95% of the parameter values respectively. You need approximately 100,000 samples to get two digits of precision on the 68% interval, and 1,000,000 samples for the 95% interval. [9]

# |
Parameter |
mean |
median |
best |
[ 68% interval] |
[ 95% interval] |
---|---|---|---|---|---|---|

1 |
M1.background |
0.059925(41) |
0.059924 |
0.059922 |
[0.05988 0.05997] |
[0.05985 0.06000] |

2 |
M1.radius |
2345.3(15) |
2345.234 |
2345.174 |
[2343.83 2346.74] |
[2342.36 2348.29] |

3 |
M1.radius.width |
0.00775(41) |
0.00774 |
0.00777 |
[ 0.0074 0.0081] |
[ 0.0070 0.0086] |

4 |
M1.scale |
0.21722(20) |
0.217218 |
0.217244 |
[0.21702 0.21743] |
[0.21681 0.21761] |

The *Convergence* plot shows the range of \(\chi^2\) values in the population
for each iteration. The band shows the 68% of values around the median, and
the solid line shows the minimum value. If the distribution has reached
equilibrium, then convergence graph should be roughly flat, with little
change in the minimum value throughout the graph. If there is no convergence,
then the remaining plots don’t mean much.

The *Correlations* plot shows cross correlation between each pair of
parameters. If the parameters are completely uncorrelated then the boxes
should contain circles. Diagonals indicate strong correlation. Square
blocks indicate that the fit is not sensitive to one of the parameters.
The range plotted on the correlation plot is determined by the 95% interval
of the data. The individual correlation plots are too small to show the
range of values for the parameters. These can instead be read from the
*Uncertainty* plot for each parameter, which covers the same range of values
and indicates 68% and 95% intervals. If there are some chains that are
wandering around away from the minimum, then the plot will look fuzzy, and
not have a nice blob in the center. If a correlation plot has multiple blobs,
then there are multiple minima in your problem space, usually because there
are symmetries in the problem definition. For example, a model fitting
\(x + a^2\) will have identical solutions for \(\pm\,a\).

The *Uncertainty* plot shows histograms for each fitted parameter generated
from the values for that parameter across all chains. Within each histogram
bar the values are sorted and displayed as a gradient from black to copper,
with black values having the lowest \(\chi^2\) and copper values having the
highest. The resulting histogram should be dark brown, with a black hump
in the center and light brown tips. If there are large lumps of light brown,
or excessive black then its likely that the optimizer did not converge. The
green line over the histogram shows the best value seen within each
histogram bin (the maximum likelihood given \(p_k == x\)).
With enough samples and proper convergence, it should roughly follow the
outline of the histogram. The yellow band in the center of the plot
represents the 68% interval for the data. The histogram cuts off at 95%.
These values along with the median are shown as labels along the x axis.
The green asterisk represents the best value, the green *E* the mean value
and the vertical green line the median value. If the fit is not sensitive
to a parameter, or if two parameters are strongly correlated, the parameter
histogram will show a box rather than a hump. Spiky shapes (either in the
histogram or the maximum likelihood line) indicate lack of convergence or
maybe not enough steps. A chopped histograms indicates that the range for
that parameter is too small.

The *Parameter Trace* plot is diagnostic for models which have poor mixing.
In this cases no matter how the parameter values are changing, they are
landing on much worse values for the \(\chi^2\). This can happen if the
problem is highly constrained with many tight and twisty values.

The *Data and Theory* plot should show theory and data lining up pretty well,
with the theory overlaying about 2/3 of the error bars on the data
(1-\(\sigma\) = 68%). The *Residuals* plot shows the difference between
theory and data divided by uncertainty. The residuals should be 2/3 within
[-1, 1], They should not show any structure, such as humps where the theory
misses the data for long stretches. This indicates some feature missing
from the model, or a lack of convergence to the best model.

If entropy is requested, then Bumps will show the total number of bits of information in the fit, where entropy is defined as:

Since we already have a sample from the posterior distribution \(p(\Theta)\) the Monte Carlo integral should be \(S \approx \sum_{k} \log_2 p(\theta_k)\). However, we do not know \(p(\theta_k)\), especially when we are integrating over nuisance parameters and only computing entropy for the parameters of interest. There are numerous methods in the literature for performing this calculation, and we have implemented the following:

*gmm*fits the MCMC sample to a Gaussian mixture model (GMM) and then estimates the entropy of the GMM through Monte Carlo integration.*llf*finds the average ratio between the unnormalized negative log likelihood (NLLF) and a kernel density estimate (sklearn*KernelDensity*with default options), then estimates the entropy from the normalized likelihood through Monte Carlo integration. [8] This technique will not work for marginal likelihood estimates.*mvn*fits the MCMC sample to a multivariate Gaussian and returns the entropy of that Gaussian. This is fast and accurate when the sample is well behaved (i.e., the uncertainty distribution is approximately Gaussian).*wnn*estimates entropy from nearest-neighbour distances in the sample. [10]

Using entropy and simulation we hope to be able to make experiment planning decisions in a way that maximizes information, by estimating whether it is better to measure more precisely or to measure different but related values and fit them with shared parameters.

### References¶

Vrugt, J. A.; Ter Braak, C. J. F.; Diks, C. G. H.; Robinson, B. A.; Hyman, J. M.; Higdon, D. International Journal of Nonlinear Sciences and Numerical Simulation, 2009, 10 (3), 273–290. DOI: 10.1515/IJNSNS.2009.10.3.273

Kramer, A.; Hasenauer, J.; Allgower, F.; Radde, N.
*In 2010 IEEE International Conference on Control Applications (CCA)*
2010; pp 493–498.
DOI: 10.1109/CCA.2010.5611198

JCGM.
*Evaluation of measurement data — Supplement 1 to the “Guide to the
expression of uncertainty in measurement” — Propagation of distributions
using a Monte Carlo method*; Joint Committee for Guides in Metrology,
JCGM 101:2008; Geneva, Switzerland, 2008; p 90.
http://www.bipm.org/utils/common/documents/jcgm/JCGM_101_2008_E.pdf

Berrett, T. B.; Samworth, R.J.; Yuan, M.;
*Efficient multivariate entropy estimation via k-nearest neighbour distances.*
Annals of Statistics 2019, 47 (1), 288-318.
DOI: 10.1214/18-AOS1688

## Particle Swarm¶

Inspired by bird flocking behaviour, the particle swarm [11] algorithm is a population-based method which updates an individual according to its momentum and a force toward the current best fit parameter values. We did not explore variations of this algorithm in any detail.

### When to use¶

Particle swarm performed well enough in our low dimensional test problems, but made little progress when more fit parameters were added.

The population updates can run in parallel, but the tiny population size limits the amount of parallelism.

### Options¶

`--steps=n`

is the number of iterations. Each step updates each member
of the population. The population size scales with the number of fitted
parameters.

`--pop=k`

determines the size of the population. The number of
individuals, \(k\), is equal to the number of fitted parameters times the
population scale factor. The default scale factor is 1.

Use `--fit=ps`

to select particle swarm from the command line.

Add a few more lines

### References¶

Kennedy, J.; Eberhart, R.
Particle Swarm Optimization
*Proceedings of IEEE International Conference on Neural Networks. IV.*
1995; pp 1942–1948.
DOI: 10.1109/ICNN.1995.48896

## Random Lines¶

Most of the population based algorithms ignore the value of the function when choosing the points in the next iteration. Random lines [12] is a new style of algorithm which fits a quadratic model to a selection from the population, and uses that model to propose a new point in the next generation of the population. The hope is that the method will inherit the robustness of the population based algorithms as well as the rapid convergence of the newton descent algorithms.

### When to use¶

Random lines works very well for some of our test problems, showing rapid convergence to the optimum, but on other problems it makes very little progress.

The population updates can run in parallel.

### Options¶

`--steps=n`

is the number of iterations. Each step updates each member
of the population. The population size scales with the number of fitted
parameters.

`--pop=k`

determines the size of the population. The number of
individuals, \(k\), is equal to the number of fitted parameters times the
population scale factor. The default scale factor is 0.5.

`--CR=v`

is the crossover ratio, determining what proportion of the
dimensions to update at each step. Values must be between 0 and 1.

`--starts=n`

tells the optimizer to restart a given number of times.
Each time it restarts it uses a random starting point.

`--keep_best`

uses a value near the previous minimum when restarting
instead of using a random value within the parameter bounds. This option is
not available in the options dialog.

Use `--fit=rl`

to select random lines from the command line.

### References¶

Sahin, I.
*An International Journal of Optimization and Control: Theories & Applications (IJOCTA)*
2013, 3 (2), 111–119.

## Parallel Tempering¶

Parallel tempering [13] is an MCMC algorithm for uncertainty analysis. This version runs at multiple temperatures simultaneously, with chains at high temperature able to more easily jump between minima and chains at low temperature to fully explore the minima. Like DREAM it has a differential evolution stepper, but this version uses the chain history as the population rather than maintaining a population at each temperature.

This is an experimental algorithm which does not yet perform well.

### When to use¶

When complete, parallel tempering should be used for problems with widely spaced local minima which dream cannot fit.

### Options¶

`--steps=n`

is the number of iterations to include in the Markov
chain. Each iteration updates the full population. The population size
scales with the number of fitted parameters.

`--burn=n`

is the number of iterations to required for the Markov
chain to converge to the equilibrium distribution. If the fit ends
early, the tail of the burn will be saved to the start of the steps.

`--CR=v`

is the differential evolution crossover ratio to use when
computing step size and direction. Use a small value to step through the
dimensions one at a time, or a large value to step through all at once.

`-nT=k`

, `-Tmin=v`

and `--Tmax=v`

specify a log-spaced initial
distribution of temperatures. The default is 25 points between
0.1 and 10. DREAM runs at a fixed temperature of 1.0.

Use `--fit=pt`

to select parallel tempering from the command line.

### References¶

Swendsen, R. H.; Wang J. S.
Replica Monte Carlo simulation of spin glasses
*Physical Review Letters*
1986, 57, 2607-2609