These routines represent an implementation of an adaptive Monte Carlo Markov Chain method to a hydrological modelling example (using the same model and data as the earlier GLUE Routines). It makes use of a statistical likelihood function based on an error model that includes a mean bias, constant variance and autocorrelation. It allows for the Box-Cox transformation of residuals to more closely match the assumptions of the error model on which the likelihood function is derived.
These routines require the R libraries CODA, MASS and LATTICE
There is also a general stand-alone R library called MCMC.
Download sample code and data file for MCMC routines [Here]
Stages in the formal Bayes Analysis (programs modified from originals written by Paul Smith, Lancaster University)
- Initialise the process: load the MCMCSetup.Rsource
This will run set-up code automatically. If you open it for editing you will see that you can change the minimum, maximum and initial values of the PDM model parameters (something to play with later, use the defaults for now)
- Fit an appropriate error transformation: load the BoxCox.Rsource.
This program carries out an initial (trial and error) estimate of the form of error model to use based on Box-Cox transformation to a zero-mean Gaussian autoregressive model. It first uses a least square optimisation of the model parameters to provide a “best-fit” model for calculating the residual series (which can take a couple of minutes). The routine provides graphics of the autocorrelation function (ACF) and quantile-quantile plots before and after transformation to illustrate how well the (optimised) transformed model errors conform to the Gaussian assumption.
You can try this with different values of the Box-Cox parameters by using the call:
Default values of the parameters are (1,0) which gives no transformation.
- Run the MCMC Identification algorithm: load the MCMCrun.Rsource
This will run automatically using the most recent values of lambda and alpha in the Box-Cox transformation but will take some time!! It is running multiple iterations of the chain in moving from the prior distribution to the posterior distribution of the parameters.
This routine provides graphics of the way in each the parameters change after each MCMC iteration, together with a plot of the 95% prediction quantiles (from both parameter and residual error) for discharge (these calculations also take some time – quantiles are being calculated over all posterior model samples at every time step).
Things to play with in these programs
Open MCMCsetup.R for editing. There you will see that the vectors of minimum, maximum and initial values of the parameters are set. The default values will give reasonable results, but try changing values to see whether it makes a big impact on the resulting error model and posterior distributions. In general, the initial parameter values will affect the resulting optimisation, so it is always worth trying multiple starting points.
You can also try out different Box-Cox transforms before running MCMCrun.R. Note, however, in doing so that the best error model at the parameter set found by a least squares optimisation might not be the best throughout the posterior density space.