Code used to conduct research at the University of Hawaii at Manoa given by Dr. Ahmed Elshall
The code presented herein is a Markov Chain Monte Carlo algorithm that runs multiple chains in parallel for efficient posterior exploration. The algorithm, entitled DREAM_(ZS) is based on the original DREAM sampling scheme, but uses sampling from an archive of past states to generate candidate points in each individual chain. Theoy and numerical examples of DREAM_(ZS) have been presented in Vrugt et al. (2009). Details can also be found in Ter Braak and Vrugt (2008)
Sampling from past has three main advantages:
(1) Circumvents the requirement of using N = d for posterior exploration. This will speed-up convergence to a limiting distribution, especially for high-dimensional problems (large d).
(2) Outlier chains do not need explicit consideration. By sampling historical states, aberrant trajectories an jump directly to the modal region at any time during the simulation. The N path ways simulated with DREAM_(ZS) therefore maintain detailed balance at every singe step in the chain.
(3) The transition kernel defining the jumps in each of the chains does not require information about the current states of the chains. This is of great advantage in a multi-processor environment where the N candidate points can be generated simultaneously so that each chain can evolve most efficiently on a different computer. Details of this will be given in a later publication, which should be ready within the next few months.
DREAM_(ZS) also contains a snooker updater to maximize the diversity of candidate points and generate jumps beyond parallel direction updates. Finally, DREAM_(ZS) contains subspace learning in a similar way as DREAM, to maximize the squared jumping distance between two subsequent points in each chain. This idea has been presented in Vrugt et al. (2008) and shown to significantly increase the efficiency of posterior exploration. All these options can be activated from the input file.
DREAM_(ZS) developed by Jasper A. Vrugt and Cajo ter Braak
SYNOPSIS:
[Sequences,X,Z,output] = dream_zs(MCMCPar,ModelName)
[Sequences,X,Z,output] = dream_zs(MCMCPar,ModelName,Extra)
[Sequences,X,Z,output] = dream_zs(MCMCPar,ModelName,Extra,ParRange)
[Sequences,X,Z,output] = dream_zs(MCMCPar,ModelName,Extra,ParRange,Measurement)Input:
MCMCPar % structure with DREAM parameters
ModelName % name of the function
Extra % optional structure with arguments to be passed to function
ParRange % optional structure with parameter ranges
Measurement % optional structure with measurement informationOutput:
Sequences % 3D array with Markov chain evolution
X % final position of chains and correponding density values
Z % matrix with thinned sample history
output % structure with convergence properties, acceptance rate, CR values, etc.The directory \PostProcessing contains a script “PostProcMCMC” that will compute various posterior statistics (MEAN, STD, MAP, CORR) and create create various plots including, marginal posterior parameter distributions, R_stat convergence diagnostic, two-dimensional correlation plots of the posterior parameter samples, chain convergence plots, and parameter and total posterior simulation uncertainty ranges (interval can be specified)
Laloy, E., and J.A. Vrugt, High-dimensional posterior exploration of hydrologic models using multiple-try DREAM_(ZS) and high-performance computing, Water Resources Research, 48, W01526, doi:10.1029/2011WR010608, 2012.
ter Braak, C.J.F., and J.A. Vrugt, Differential Evolution Markov Chain with snooker updater and fewer chains, Statistics and Computing, 10.1007/s11222-008-9104-9, 2008
Vrugt, J.A., E. Laloy, and C.J.F. ter Braak, DiffeRential Evolution Adaptive Metropolis with Sampling from the Past and Subspace Updating, SIAM journal on Optimization
For more information please read:
Vrugt J.A., H.V. Gupta, W. Bouten and S. Sorooshian, A Shuffled Complex Evolution Metropolis algorithm for optimization and uncertainty assessment of hydrologic model parameters, Water Resour. Res., 39 (8), 1201, doi:10.1029/2002WR001642, 2003.
ter Braak, C.J.F., A Markov Chain Monte Carlo version of the genetic algorithm Differential Evolution: easy Bayesian computing for real parameter spaces, Stat. Comput., 16, 239 - 249, doi:10.1007/s11222-006-8769-1, 2006.
Vrugt, J.A., C.J.F. ter Braak, M.P. Clark, J.M. Hyman, and B.A. Robinson, Treatment of input uncertainty in hydrologic modeling: Doing hydrology backward using Markov chain Monte Carlo, Water Resour. Res., 44, W00B09, doi:10.1029/2007WR006720, 2008.
Vrugt, J.A., C.J.F. ter Braak, C.G.H. Diks, D. Higdon, B.A. Robinson, and J.M. Hyman, Accelerating Markov chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling, International Journal of Nonlinear Sciences and Numerical Simulation}, 10(3), 273-290, 2009.
This program is free software: you can modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
Version 0.5: January 2009
Version 1.0: April 2011 Maintenance update, explicit treatment of prior distribution
Version 1.1: August 2011 Whittle likelihood function (SPECTRAL ANALYSIS !!)
Version 1.2: April 2012 Simplified code (removed variables) + graphical interface
Version 1.3: June 2012 Simulations stored, new example, and updated likelihood func.
Version 1.4: January 2013 Simplification of metrop.m and dream_zs.m
Example 1: N-dimensional banana shaped Gaussian distribution
Example 2: N-dimensional Gaussian distribution
Example 3: N-dimensional multimodal mixture distribution
Example 4: Real-world example using hymod rainfall - runoff model (HYMOD code in MATLAB)
Example 5: Real-world example using hymod rainfall - runoff model (HYMOD code in FORTRAN)
Example 6: Rainfall-runoff model with generalized log-likelihood function
Example 7: HYDRUS-1D soil hydraulic model: using prior information on soil hydraulic parameters
Example 8: Multivariate student t distribution
Example 9: Rainfall-runoff model with Whittle’s likelihood function
Example 10: The use of prior information in a multimodel mixture distrbibution
Example 11: Multivariate student t distribution
Example 12: Pedometrics problem involving variogram fitting
Example 13: Nash-Cascade example –> heteroscedastic errors
Example 14: ABC inference for hydrologic model
Example 15: ABC inference using 10 bivariate normal distributions
Example 16: Hydrogeophysics example
Example 17: 81 Parameter Ground water PDE