R: Fit a hydromad model using the DREAM (DiffeRential Evolution...

fitByDream {hydromad}

R Documentation

Fit a hydromad model using the DREAM (DiffeRential Evolution Adaptive Metropolis) algorithm.

Description

Fit a hydromad model using the DREAM (DiffeRential Evolution Adaptive Metropolis) algorithm. This is a Markov Chain Monte Carlo algorithm which gives estimates of the joint probability distribution of parameters according to a likelihood function. The fitting function returns the maximum likelihood model, but the full MCMC results are also available as component $fit.result. The result can also be used to define a feasible parameter set.

Usage

fitByDream(MODEL, loglik = hydromad.getOption("loglik"),
           control = hydromad.getOption("dream.control"),
           vcov = TRUE,save=NULL)

Arguments

`MODEL`	a model specification created by `hydromad`. It should not be fully specified, i.e one or more parameters should be defined by ranges of values rather than exact values.
`loglik`	log-likelihood function (log of the posterior probability density), given as a `function(Q, X, ...)`. See `objFunVal`.
`control`	settings for the DREAM algorithm. See `dream`.
`vcov`	if `vcov = TRUE`, the parameter variance-covariance matrix will be estimated from the last half of the sequences. It can be extract using `vcov`.
`save`	Optional `function(pars,objective value,model)` that will be called for every model evaluation, for example to save every model run.

Value

the best model from those sampled, according to the given loglik function. Also, these extra elements are inserted:

`fit.result`	the result from `dream`.
`objective`	the `loglik` function used.
`funevals`	total number of evaluations of the model simulation function.
`timing`	timing vector as returned by `system.time`.

Author(s)

Felix Andrews felix@nfrac.org

Examples


if (require("dream", quietly = TRUE)) {

data(Cotter)
x <- Cotter[1:1000]

## IHACRES CWI model with power law unit hydrograph
modx <- hydromad(x, sma = "cwi", routing = "powuh")
modx

## a very short run! just to demonstrate methods
foo <- fitByDream(modx, control = list(ndraw = 500))

summary(foo)

## parameter correlation matrix with symbols
symnum(cov2cor(vcov(foo)))

## return value from dream:
str(foo$fit.result)

## plot log-likelihood value convergence over time
xyplot(window(optimtrace(foo, raw = TRUE), start = 50),
  superpose = TRUE, auto.key = FALSE,
  xlab = "function evaluations", ylab = "neg. log likelihood")

## calculate corresponding objective function values over time.
xyplot(optimtrace(foo, objective = ~ - hmadstat("r.squared")(Q,X)),
  xlab = "function evaluations", ylab = "negative R Squared")

## MCMC diagnostics and more are available:
methods(class = "dream")

}

[Package hydromad version 0.9-18 Index]