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

See Also

dream, objFunVal

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]