added page for DREAM

Author: alexhroom

source/algorithms/DREAM.rst

@@ -4,19 +4,87 @@
 Bayesian Analysis - DREAM
 =========================
 
-The main Bayesian Sampler bundled with RAT is an implementation of the DREAM (DiffeRential Evolution Adaptive Metropolis)
-algorithm first described by `Vrugt`_. 
+The main Bayesian Sampler in RAT is a MATLAB implementation of the DREAM (DiffeRential Evolution Adaptive Metropolis)
+algorithm of Vrugt (2016), [#vrugt2016]_ modified to be compilable to C++ via MATLAB Coder.
 
+As the name suggests, DREAM combines :ref:`differential evolution<DE>` with a variant on the 
+`Metropolis-Hastings optimisation algorithm <https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm>`_.
 
+The basic idea of the Metropolis-Hastings algorithm is that it randomly "walks" around the parameter space, by taking
+a "step" of a fixed distance and either accepting the step or rejecting it (i.e. going back to the previous step and starting again)
+based on the following process:
 
+#. Let :math:`x` be the current point, :math:`x'` be the proposed new step, and :math:`f` the function that calculates chi-squared at any point.
+#. Calculate the ratio
 
+   .. math:: \alpha = f(x')/f(x).
 
+#. Generate a random number :math:`u` between 0 and 1.
+#. If :math:`u \leq \alpha`, accept the candidate. Otherwise, reject it.
 
+This means that while steps into points with better chi-squared are more likely, they are not *guaranteed*, and the algorithm is also able to step
+to points of worse fit (but this is less likely to be accepted). The sequence of steps is an example of a Markov chain. After a large number of steps
+we can infer various things about the shape of our probability distribution for each parameter. The algorithm can be adapted to, for example, adapt
+the size of the steps to shrink proportionally to the number of steps that have been accepted, with the assumption that the algorithm has found a minimum
+after many accepted steps.
 
+DREAM modifies this idea and combines it with differential evolution. It creates multiple Markov chains ("walks") in parallel,
+and the new proposed points are generated by using a differential evolution process between all the chains:
 
+#. Let :math:`v_i` be the vector representing chain :math:`i`; 
+   :math:`v_i = (x^i_1, x^i_2, x^i_3,...)` 
+   where :math:`x^i_j` is the :math:`j`'th point on chain :math:`i`.
+#. Let there be :math:`N` chains, and say we are on iteration :math:`t`.
+   Use differential evolution over all points on all chains to propose a new point for each chain:
+   :math:`x^1_t, x^2_t, ..., x^N_t`.
+#. For each point, either accept or reject using the same criteria as in Metropolis-Hastings. If
+   accepted, add it to its respective chain.
 
+This means that over time, DREAM creates :math:`N` parallel Markov chains that explore the parameter space,
+and each chain can 'learn' from better-performing chains via the differential evolution process. After some iterations,
+we start 'cutting off' the older end of the chain so that the new proposals are generated from a sample set of better points.
 
+It also performs several auxilliary improvements such as subspace sampling (i.e. generating proposals from only a subset of the parameters
+at a time), and a correction algorithm for chains that are outliers. It has been found to outperform other adaptive
+MCMC approaches, and in some cases even outperforms classical optimisation algorithms at parameter estimation. [#laloy2012]_
 
 
+Algorithm control parameters
+----------------------------
+The following parameters in the :ref:`Controls object<controlsInfo>` are specific to DREAM:
 
-.. _Vrugt: https://www.sciencedirect.com/science/article/abs/pii/S1364815215300396
+- ``nSamples``: The number of samples in the initial population for each chain.
+
+- ``nChains``: The number of Markov chains to use in the algorithm.
+
+- ``jumpProbability``: A probability range for the size of jumps when performing subspace sampling. Larger values mean
+  more random variation in jump sizes for new proposed values.
+
+- ``pUnitGamma``: The 'probability of unit jump' rate. New proposals are scaled down by a factor :math:`\gamma = 2.38/2d`,
+  where :math:`d` is the current length of the chain. This means that as the chain gets longer and 'better', jumps to new
+  points aren't as far from the rest of the chain (as we are likely close to a minimum of the space).
+  ``pUnitGamma`` is the probability that this will be ignored for an iteration
+  and instead :math:`\gamma = 1` for that iteration (hence the name: "p unit gamma", the probability that gamma
+  will be unital for an iteration). This allows chains to 'jump' out of their local area and more easily explore other
+  parts of the parameter space.
+
+- ``boundHandling``: How steps across the boundary of the space should be handled. Can be:
+  - ``"off"``: allow steps to pass the boundary. 
+  - ``"reflect"``: if a step passes a boundary, reflect it back across the boundary.
+  - ``"bound"``: if a step passes a boundary, set it equal to the nearest point on the boundary.
+  - ``"fold"``: treat the boundary as periodic and 'wrap the step around' to the other side of the space.
+
+- ``adaptPCR``: whether crossover probabilities for differential evolution should be adapted by the algorithm as it runs.
+
+
+.. [#vrugt2016] 
+   Vrugt, Jasper (2016),
+   "Markov chain Monte Carlo simulation using the DREAM software package: Theory, concepts, and MATLAB implementation".
+   DOI: 10.1016/j.envsoft.2015.08.013,
+   URL: https://www.sciencedirect.com/science/article/abs/pii/S1364815215300396
+
+.. [#laloy2012]
+   Laloy, Eric; Vrugt, Jasper (2012),
+   "High-dimensional posterior exploration of hydrologic models using multiple-try DREAM(ZS) and high-performance computing".
+   DOI: 10.1029/2011WR010608,
+   URL: https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2011WR010608