DREAM ( D ) : an adaptive markov chain monte carlo simulation algorithm to solve discrete , noncontinuous , posterior parameter estimation problems

Introduction Conclusions References


Introduction
Formal and informal Bayesian methods have found widespread application and use to summarize parameter and model predictive uncertainty in hydrologic modeling.These parameters generally represent model dynamics, but could also include rainfall multipliers (Kavetski et al., 2006;Kuczera et al., 2006;Vrugt et al., 2008), error model variables (Smith et al., 2008;Schoups and Vrugt, 2010), and calibration data measurement errors (Sorooshian and Dracup, 1980;Schaefli et al., 2007;Vrugt et al., 2008).Monte Carlo methods are admirably suited to generate samples from the posterior parameter distribution, but generally inefficient when confronted with complex, multimodal, Figures and high-dimensional model-data synthesis problems.This has stimulated the development of Markov Chain Monte Carlo (MCMC) methods that generate a random walk through the search (parameter space) and iteratively visit solutions with stable frequencies stemming from an invariant probability distribution.If well designed, such MCMC methods should be more efficient than brute force Monte Carlo or importance sampling methods.
To visit configurations with a stable frequency, an MCMC algorithm generates trial moves from the current position of the Markov chain x t−1 to a new state z.The earliest and most general MCMC approach is the random walk Metropolis (RWM) algorithm (Metropolis et al., 1953).Assume that we have already sampled points {x 0 ,...,x t−1 } this algorithms proceeds in the following three steps.First, a candidate point z is sampled from a proposal distribution q that depends on the present location, x t−1 and is symmetric, q(x t−1 ,z) = q(z,x t−1 ).Next, the candidate point is either accepted or rejected using the Metropolis acceptance probability: where π(•) denotes the probability density function (pdf) of the target distribution.Finally, if the proposal is accepted the chain moves to z, otherwise the chain remains at its current location x t−1 .The result is a Markov chain which under some regularity conditions has a unique stationary distribution with pdf π(•).
The standard RWM algorithm has been designed to maintain detailed balance with respect to π(•) at each step in the chain: where π(x t−1 ) (π(z)) denotes the probability of finding the system in state x t−1 (z), and p(x t−1 → z)(p(z → x t−1 )) denotes the conditional probability to perform a trial move from x t−1 to z (z to x t−1 ).Hastings (Hastings, 1970) extended Eq. ( 1) to include nonsymmetrical proposal distributions, i.e. q(x t−1 ,z) = q(z,x t−1 ) in which a proposal jump Figures to z and the reverse jump do not have equal probability.This extension is called the Metropolis Hastings algorithm (MH), and has become the basic building block of many existing MCMC sampling schemes.
Existing theory and experiments prove convergence of well-constructed MCMC schemes to the appropriate limiting distribution under a variety of different conditions.
In practice, this convergence is often observed to be impractically slow.This deficiency is frequently caused by an inappropriate selection of q(•) used to generate trial moves in the Markov Chain.This inspired Vrugt et al. (2008Vrugt et al. ( , 2009) ) to develop a simple adaptive RWM algorithm called Differential Evolution Adaptive Metropolis (DREAM) that runs multiple chains simultaneously for global exploration, and automatically tunes the scale and orientation of the proposal distribution during the evolution to the posterior distribution.This scheme is an adaptation of the Shuffled Complex Evolution Metropolis (Vrugt et al., 2003) global optimization algorithm and has the advantage of maintaining detailed balance and ergodicity while showing excellent efficiency on complex, highly nonlinear, and multimodal target distributions (Vrugt et al., 2008(Vrugt et al., , 2009)).
In DREAM, N different Markov Chains are run simultaneously in parallel.If the state of a single chain is given by a single d-dimensional vector x, then at each generation the N chains in DREAM define a population X, which corresponds to an N x d matrix, with each chain as a row.Jumps in each chain i = {1,...,N} are generated by taking a fixed multiple of the difference of the states of randomly chosen pairs of chains of X: where δ signifies the number of pairs used to generate the proposal, x r 1 (j ) and x r 2 (n) are randomly selected without replacement from the population X d , the number of dimensions that will be updated jointly.By comparison with RWM, a good choice for γ = 2.4/ √ 2δd (Roberts and Rosenthal, 2001;Ter Braak, 2006).This choice is expected, for Gaussian and Student target distributions, to yield an acceptance probability of 0.44 for d = 1, 0.28 for d = 5 and 0.23 for large d .Every 5th generation γ = 1.0 to facilitate jumping between disconnected posterior modes (Vrugt et al., 2008).
The difference vector in Eq. ( 3) contains the desired information about the scale and orientation of the target distribution, π(x|•).By accepting each jump with the Metropo- , a Markov chain is obtained, the stationary or limiting distribution of which is the posterior distribution.The proof of this is given in Ter Braak and Vrugt (2008) and Vrugt et al. (2008Vrugt et al. ( , 2009)).Because the joint pdf of the N chains factorizes to π(x 1 |•) × ... × π(x N |•), the states x 1 ...x N of the individual chains are independent at any generation after DREAM has become independent of its initial value.After this burn-in period, the convergence of DREAM can thus be monitored with the R-statistic of Gelman and Rubin (1992).This convergence diagnostic compares the within and in-between variances of the N different chains.
Various recent contributions have shown the utility of MCMC simulation with DREAM to treat different error sources, and help quantify and analyze parameter, model structural, forcing, and calibration data uncertainty (Vrugt et al., 2008(Vrugt et al., , 2009b;;Dekker et al., 2010;He et al., 2010;Huisman et al., 2010;Keating et al., 2010;Minasny et al., 2011;Schoups and Vrugt, 2010;Vrugt et al., 2011).These, and most other model-data synthesis studies in the earth sciences, typically involve continuous variables (parameters and probability density functions) that can take on any numerical value within their prior ranges defined by the user.Yet, relatively little attention has been given to posterior sampling problems involving discrete variables.The existing DREAM framework has been developed based on the assumption of continuity of the parameter space, lim in many fields of study, and therefore of considerable theoretical and practical interest.For instance, in hydrology there is increasing interest in using optimization approaches to help find optimal experimental design strategies that attempt to minimize cost, parameter and model predictive uncertainty, or combinations thereof.This involves selecting one or multiple different measurement locations amongst a discrete set of possibilities.The solution of this problem lies in an adaptation of DREAM, and its various extensions including DREAM (ZS) (Vrugt et al., 2011), and MT-DREAM (ZS) (Laloy and Vrugt, 2011).
In this paper we present a discrete implementation of DREAM, that is especially designed to solve noncontinuous search and optimization problems.This new code, DREAM (D) uses DREAM as its main building block, yet implements integer search to facilitate solving discrete sampling problems.The DREAM (D) algorithm maintains detailed balance and ergodicity, and achieves excellent performance across a range of noncontinuous posterior sampling problems.The DREAM (D) code is perhaps the only method available to solve discontinuous parameter estimation problems, and simultaneously provide an estimate of uncertainty.
The remainder of this paper is organized as follows.Section 2 presents a short introduction to MCMC, followed by a detailed description of DREAM (D) in Sect.3. Section 4 demonstrates the performance of DREAM (D) using two different case studies involving a simple 52-dimensional sudoku puzzle, and a rainfall-runoff model calibration problem.
These results are illustrated in detail, and used to highlight further possible improvements.Finally, in Sect. 5 we summarize the methodology and discuss the results.

Nonlinear optimization involving discrete variables
Discrete optimization problems are abundant in many fields of study, and have begun to appear in the hydrologic literature (Furman et al., 2004;Harmancioglu et al., 2004;Perrin et al., 2008;Kleidorfer et al., 2008;Neuman et al., 2011).To illustrate a noncontinuous parameter estimation problem, please consider Other problems that require integer search involve finding optimal experimental design strategies.Usually this consists of finding the best measurement locations among a prior defined and often restricted set of possibilities.This requires integer optimization in a similar way as done in the tile puzzle.

DREAM (D) ⇒ differential evolution adaptive metropolis with discrete sampling
We now describe our new code, entitled DREAM (D) , which uses DREAM as main building block.
Let X be a N × d -matrix defining the N initialized starting positions, x i , i = 1,...,N of the different Markov chains by drawing samples from p d (x), the prior distribution.
Similarly, let Z be an external archive of points that periodically appends the elements of X at regular intervals.The initial population [X t ;t = 0] is translated into a sample from the posterior target distribution using the following pseudo code:

Conclusions References
Tables Figures

Back Close
Full a. Generate a candidate point, z i in chain i , where the function • d rounds each element j = 1,...,d of the jump vector to the nearest integer.
b. Replace each element (j = 1,...,d ) of the proposal z i j with z i j using a binomial scheme with probability 1 − CR, where CR denotes the crossover probability, and U ∈ [0,1] is a draw from a uniform distribution.
c. Compute π(z i ) and accept the candidate points with Metropolis acceptance probability, α(x i ,z i ), 5. Summarize the posterior pdf using Z after discarding the initial and burn-in samples.
The DREAM (D) algorithm presented herein is similar as DREAM, but especially designed to solve discrete search and optimization problems.The function • d in Eq. ( 4) forces the jumps to maintain integer values, and maintains detailed balance (as demonstrated later).Compared to the original DREAM sampling scheme, DREAM (D) contains two additional algorithmic parameters.These are T max , the maximum number of generations, and K , the thinning rate used to periodically add samples to Z. Based on recommendations in (Vrugt et al., 2011), we set K = 10, and assign a large value for T max , thus DREAM (D) automatically stops after convergence has been achieved.To speed up convergence to the target distribution, DREAM estimates a distribution of CR values during burn-in that favors large jumps over smaller ones in each of the N chains.Details of this can be found in (Vrugt et al., 2008(Vrugt et al., , 2009(Vrugt et al., , 2011)).
The transition kernel in DREAM (D) is especially designed to handle discrete search and optimization problems, yet only samples integer values.This is by no means a limitation, but a simple computational trick is therefore required to solve non-integer problems.For such problems, we discretize the range of each individual parameter into equidistant intervals, and rewrite the actual problem using integers only.The posterior distribution of these integers is subsequently derived with DREAM (D) .For example, consider a discrete parameter that can only take on the values x 1 ∈ [0,0.2,0.4,...,5].We can rewrite this problem as: x 1 = min(x 1 ) + 0.2j and sample j ∈ [0,1,...,25] using DREAM (D) .This approach, if deemed necessary, allows for a different discretization of each individual parameter.Also, the simultaneous use of DREAM and DREAM (D) , enables joint sampling of continuous and discrete parameters.We are now left with a proof that DREAM (D) yields an invariant distribution that is identical to the posterior distribution.For this we need to demonstrate that the transition kernel in DREAM maintains detailed balance, and thus results in a reversible Markov chain.This essentially means that the forward (p(x i → z i )) and backward (p jump have equal probability at every single step in the chain.This is easy to proof for standard RWM algorithms that use a fixed proposal distribution.Yet, the jumping distribution used in DREAM (D) is adaptive, and continuously changes scale and orientation en route to the posterior distribution.This significantly enhances efficiency, but it is not immediately clear whether Eq. ( 4) also satisfies the detailed balance condition.
In previous papers, we have given formal proofs of convergence of DREAM (Vrugt et al., 2008), DREAM (ZS) (Vrugt et al., 2011), and MT-DREAM (ZS) (Laloy and Vrugt, 2011) to the appropriate limiting distribution.For simplicity, we use a hypothetical example to illustrate that DREAM (D) maintains detailed balance.Lets consider a twodimensional, (d = 2) discrete sampling problem with N = 5 different chains, δ = 1, and The initial states of these chains are color coded and depicted in Fig. 2, and used to generate candidate points.Lets assume that the jump in the first chain (purple) uses the states of the 2nd (r 1 : blue), and 4th (r 2 : green) chain, and that e j = j = 0; j = 1,...,2.Following Eq. ( 4), the forward jump and proposal point then becomes: This point is indicated in Fig. 2 with the red square.Lets assume that we accept this candidate point, and chain 1 transitions to this new state.
We are now left with studying the probability of the backward jump.This requires the green chain to be r 1 , and the blue chain to be r 2 , thus selected in opposite order respective chain pairs for the candidate point.Hence, the chance to select r 1 as chain 2, and r 2 as chain 4 is equal to drawing r 2 = 4, and r 2 = 2. Reversibility is thus ensured, yet is is not directly obvious that this proof also holds for DREAM (D) that uses an explicit integer rounding function in the transition kernel.If we proceed with our new selection of r 1 and r 2 then the next candidate point becomes: Indeed, the reverse jump results in a proposal point that is identical to the initial state of chain 1 (purple point), and reversibility is ensured.Apparently, the function, • d , used to round the d -dimensional jumping vector to the closest integers does not violate detailed balance.This proof also holds when e d , and d are drawn from their respective symmetric probability distributions, and when δ > 1 and d > 2. We leave this up to the reader.This concludes the proof of detailed balance.
A final remark is appropriate.In theory it is possible that at least one of the d arguments, y d of the rounding function, y d has a fractional part of .5.Or in mathematical notation, ∃y j ∈ Y j : y j − y j = 0.5;j = 1,...,d , where Y j denotes the feasible jump space of the j th dimension, generated with Eqs. 4 Case studies I now present two different case studies with increasing complexity.The first study consists of a typical integer estimation problem, and involves a sudoku puzzle.This puzzle has become quite popular in the past 10 years, and many newspapers, journals, and magazines around the world publish sudokus for entertainment.This synthetic study illustrates the ability of DREAM (D) to help solve a relatively difficult integer optimization problem.The second study is concerned with estimating parameters in a mildly complex lumped watershed model using observed daily discharge data from the Guadalupe River in Texas.The parameter estimation problem is posed in discrete form, and solved using DREAM (D) .This results in a discrete The posterior parameter distribution is compared against a classical continuous formulation of the model calibration inverse problem separately inferred using DREAM (ZS) , and used to highlight the advantages of DREAM (D) and discrete sampling.

The daily sudoku
The first case study considers a synthetic integer parameter estimation problem, involving a sudoku puzzle.The objective is to fill a 9 × 9 grid with digits so that each column, each row, and each of the nine 3 × 3 sub-grids that make up the total square contains values of 1 to 9. The same single integer may not appear twice in the same 9 × 9 playing board row or column or in any of the nine 3 × 3 subregions.The puzzle setter provides a partially completed grid, which typically has a single (unique) solution.
The puzzle was popularized in 1986 by the Japanese puzzle company Nikoli, under the name Sudoku, meaning single number (Hayes, 2006).Nowadays, Sudoku puzzles are very popular, and widely practised by many millions of people throughout the world.I consider the Sudoku puzzle in Fig. 3, taken from Wikipedia (http://en.wikipedia.org/wiki/Sudoku).The initial grid is depicted at the left-hand side, whereas the final solution is presented at the right-hand side.In this particular puzzle, the solution at estimated.These parameters can take on values between 1 to 9. Figure 4 illustrates the sampled values of DREAM (D) at various stages during the integer search.A total of N = 25 different chains were used to search the parameter space, and the initial sample was created using Latin Hypercube Sampling.I did not impose any constraints on this initial solution, and thus each value of 1 to 9 can appear in any number in the grid.The log-likelihood function measures the constraint violation, details of which are outside the scope of this publication.
The first grid, at the left-hand side of Fig. 4 presents an example of a possible starting solution used in one of the chains.Obviously, this solution is rather bad, and violates each of the three different constraints discussed previously (1 to 9 in each column and row, and each 3 × 3 subregion).The second grid, Fig. 4b, shows considerable improvements, and better resemblances the final solution.The third panel in Fig. 4c is a further refinement, but with a few noticeable deviations from the true solution.Finally, after about 1 million sudoku function evaluations, the puzzle is successfully solved (Fig. 4d).
Obviously, it is inspiring that DREAM (D) is able to successfully solve a Sudoku puzzle, but the required number of function evaluations (computational time) to do so is rather large.Indeed, a branch and bound optimization approach would be about 3 orders of magnitude more efficient.Such an approach shuffles the d = 52 integer values until all constraint are met.We could implement this procedure in DREAM (D) , and change the proposal distribution in such a way that jumps in each individual chain are no longer created by taking the difference of the states of two (or more) other chains but simply generated by shuffling selected parameter values of the same chain.Yet, it is not directly obvious how to make such modification and ensure reversibility of the Markov chain.This is an interesting problem, and deserves additional investigation, and theoretical development.

Watershed model calibration using discrete parameter estimation
The final case study involves flood forecasting, and consists of the calibration of a mildly complex lumped watershed model using historical data from the Guadalupe River at Spring Branch, Texas.This is the driest of the 12 MOPEX basins described in the study of Duan et al. (2006).The model structure and hydrologic process representations are found in (Schoups and Vrugt, 2010).The model transforms rainfall into runoff at the watershed outlet using explicit process descriptions of interception, throughfall, evaporation, runoff generation, percolation, and surface and subsurface routing.Table 1 summarizes the seven different model parameters, and their prior uncertainty ranges.Each parameter is discretized equidistantly in 250 intervals with respective step size listed in the last column at the right hand side.This gridding is necessary to create a non-continuous, discrete, parameter estimation problem.Unlike the previous case study in which integer values are sampled only, this particular study (mostly) involves non-integer values.Daily discharge, mean areal precipitation, and mean areal potential evapotranspiration were derived from (Duan et al., 2006) and used for model calibration.Details about the basin, experimental data, and likelihood function can be found there, and will not be discussed herein.The same model and data was used in a previous study (Schoups and Vrugt, 2010), and used to introduce a generalized likelihood function for heteroscedastic, non-Gaussian, and correlated (streamflow) prediction errors.
Figure 5 presents histograms of the marginal distribution of a selected few hydrologic model parameters using five years of observed daily discharge data.The top panel presents the results of DREAM (D) , whereas the bottom panel presents the results for a continuous parameter space.These histograms were derived by separately running DREAM for the same data set and model.Notice the close agreement between the histograms derived with both MCMC methods.This is a testament to the ability of DREAM (D) to successfully solve discrete posterior parameter estimation problems.The influence of gridding is hardly noticeable, but becomes apparent if we use at least 25 bins to represent the marginal density (not shown herein).To better illustrate the effect of discretization, please consider Figure 6 that presents two-dimensional scatter plots of the DREAM (D) derived posterior samples for a few selected parameter pairs.The bottom panel shows similar plots but then assuming continuity of the parameter space.The effect of gridding is immediately apparent.Whereas the original bivariate scatter plots sample the parameter space in a (bloodstain) spatter pattern, two-dimensional plots of the posterior samples derived with DREAM (D) exhibit an obvious grid pattern with horizontally and vertically aligned points.The posterior samples take on discrete values with a distance between subsequent points that is similar to the intervals listed in Table 1.Despite this difference in sampling pattern, the shape of the DREAM and DREAM (D) derived bivariate scatter plots are very similar, commensurate with the covariance structure of the posterior distribution.The results presented in Fig. 5 inspire confidence in the ability of DREAM (D) to solve noncontinuous posterior sampling problems.
The excellent correspondence of the posterior parameter distributions derived with DREAM and DREAM (D) results in marginal differences of the resulting streamflow predictions.I therefore do not show any times series plots of model predictions, and corresponding 95% uncertainty ranges.Such plots can be found in (Schoups and Vrugt, 2010) and details can be found in that publication.It is interesting to observe that the maximum log-likelihood value of 543 found with DREAM (D) is somewhat larger than its counterpart estimated with DREAM (540).This difference was consistently observed for multiple different trials with both MCMC algorithms.
To provide better insights into the efficiency of DREAM (D) , Fig. 7a plots the evolution of the R-statistic of Gelman and Rubin (Gelman and Rubin, 1992) for each of the different parameters.To benchmark these results, the bottom panel (Fig. 7b feasible parameter space, this has no immediate effect on sampling efficiency.An approximately similar number of model runs is required to converge to and explore the posterior target distribution.Yet, this by no means is a universal finding.My initial results to date for more complex models with much higher parameter dimensionality demonstrate considerable enhancements in efficiency (sometimes dramatically) when solving the calibration problem in the discretized rather than continuous space.Such transformation is particularly useful for insensitive parameters as gridding significantly reduces the space of feasible solutions.Further research on this topic is warranted.

Conclusions
In the past decade much progress has been made in the development of sampling algorithms for statistical inference of the posterior parameter distribution.The typical assumption in this work is that the parameters are continuous and can take on any value within their upper and lower bounds.Unfortunately, such problems typically do not work for discrete parameter estimation problems.Such problems are abundant in many fields of study, and therefore of considerable theoretical and practical interest.
Examples include selecting among different measurement locations in the design of optimal experimental strategies, finding the best members of an ensemble of predictors, and more generally discrete model calibration problems.Here, I have introduced a discrete MCMC simulation algorithm that is especially designed to solve non-continuous posterior sampling problems.This method, entitled DREAM (D) uses DREAM as its main building block, yet uses a modified proposal distribution to facilitate solve discrete sampling problems.The DREAM (D) algorithm maintains detailed balance and ergodicity, and receives good performance across a range of problems involving a sudoku puzzle, and discrete rainfall -runoff model calibration problem.
The theory developed herein is easily implemented in DREAM (ZS) (Vrugt et al., 2011) and MT-DREAM (ZS) (Laloy and Vrugt, 2011), which provides a venue to further increase the efficiency of MCMC simulation.The DREAM (D) code is written in MATLAB and is available upon request (jasper@uci.edu).Figures

Back Close
Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | without x i t−1 ); r 1 (j ),r 2 (n) ∈ {1,...,N} and r 1 (j ) = r 2 (n).The values of e d and d are drawn from U d (−b,b) and N d (0,b Discussion Paper | Discussion Paper | Discussion Paper | x→c π(x) = π(c), and thus will exhibit some practical problems when confronted with discrete variables.Such non-continuous parameter estimation problems are abundant Discussion Paper | Discussion Paper | Discussion Paper |

Fig
Discussion Paper | Discussion Paper | Discussion Paper |puzzle consisting of sixteen different surfaces.Each tile contains a different letter from the alphabet.The goal is to get the letters in the appropriate order.The solution to this problem is immediately obvious to a human, but not immediately clear to a computer.A search algorithm is therefore required to solve this problem.If we assign numbers to each letter, a = 1, b = 2 and so forth, we could measure the distance from our initial guess to the actual solution.This constitutes an integer optimization problem.The tiles can only take on integer values, between 1 and 16.This results in a d = 16 dimensional search problem with each dimension x ∈ [1,2,...,16].Figure1aillustrates an initial guess, and in a series of panels (B ⇒ D) it is shown how DREAM (D) translates this guess into the final solution.
Discussion Paper | Discussion Paper | Discussion Paper | d.If accepted, move the chain to the candidate point, x i = z i , otherwise remain at the old location, x i .END FOR (CHAIN EVOLUTION) END FOR (POPULATION EVOLUTION) 2. Append X to Z. 3. Compute the Gelman-Rubin convergence diagnostic, Rj (Gelman and Rubin, 1992) for each dimension j = 1,...,d using the last 50% of the samples of Z. Discussion Paper | Discussion Paper | Discussion Paper | 4. If Rj ≤ 1.2 for j = 1,...,d or T > T max , stop and go to step 5, otherwise set T ← T +1 and go to POPULATION EVOLUTION.
from the previous jump.The chance of this backward jump is equal to the chance of the forward jump because of the uniform random number generator used to select the Discussion Paper | Discussion Paper | Discussion Paper | (3) or (4).If this (mathematical) statement is true, then the respective argument(s) is (are) rounded down (floor) to the closest integer.This directed rounding violates detailed balance and introduces a possible bias in the proposal jump.A simple one-dimensional example will immediately illustrate this, but the chance this bias happens in practice is virtually zero.If nothing else, the stochastic nature of e d ∼ U d (−b,b), and d ∼ N d (0,b * ) will eliminate this possibility.In the rare event of y j − y j = 0.5 we implement stochastic rounding, and choice among y j − 0.5 and y j + 0.5 with equal probability.This ensures reversibility of the Markov chains generated with DREAM (D) .Note that the size of Y j essentially depends on the choice of the prior distribution.Discussion Paper | Discussion Paper | Discussion Paper | 29 different cells is known, leaving us with d = 81 − 29 = 52 parameter values to be Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ) illustrates the convergence results of DREAM assuming continuity of the parameter space.The presented traces represent an average of 25 different MCMC trials.The convergence speed for both algorithms is strikingly similar.Both MCMC methods need about 30 000 model evaluations to converge to a stationary (limiting) distribution (indicated with the dashed blue horizontal line).Although gridding significantly reduces the size of the Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 1 .Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .
Fig. 1.A 4 × 4 square tile puzzle with 16 different letters.The goal is to get the letters in order of the alphabet.Each letter is assigned a different integer value, and the resulting inverse problem is solved using DREAM (D) .The panels (B)-(D) show the various steps of DREAM (D) to the target solution.

Fig. 7 .
Fig. 7. Simulated traces of the R-statistic of Gelman and Rubin (Gelman and Rubin, 1992) using the (A) DREAM (D) (top panel)(discrete sampling), and (B) DREAM (bottom panel) (continuous sampling) MCMC algorithms.Each parameter is coded with a different color.A similar convergence speed is observed for both algorithms.

Table 1 .
Prior Uncertainty Ranges of Hydrologic and Error Model Parameters.