R/fitByDream.R
fitByDream.RdFit 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.
fitByDream(
MODEL,
loglik = hydromad.getOption("loglik"),
control = hydromad.getOption("dream.control"),
vcov = TRUE,
save = NULL
)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.
log-likelihood function (log of the posterior probability
density), given as a function(Q, X, ...). See
objFunVal.
settings for the DREAM algorithm. See
dream.
if vcov = TRUE, the parameter variance-covariance matrix
will be estimated from the last half of the sequences. It can be extract
using vcov.
Optional function(pars,objective value,model) that will
be called for every model evaluation, for example to save every model run.
the best model from those sampled, according to the given
loglik function. Also, these extra elements are inserted:
the result from dream.
the loglik function used.
total number of evaluations of the model simulation function.
timing vector as returned by system.time.
if (requireNamespace("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")
}
#> Warning: Outliers detected outside burn-in period, returning to burn in
#> Warning: Outliers detected outside burn-in period, returning to burn in
#> Warning: Outliers detected outside burn-in period, returning to burn in
#> Warning: Outliers detected outside burn-in period, returning to burn in
#> Warning: Outliers detected outside burn-in period, returning to burn in
#> List of 15
#> $ call : language dream::dream(FUN = do_dream, func.type = "logposterior.density", pars = parlist, control = control)
#> ..- attr(*, ".Environment")=<environment: 0x7f830fc84528>
#> $ control :List of 19
#> ..$ nCR : num 3
#> ..$ gamma : num 0
#> ..$ steps : num 10
#> ..$ eps : num 0.05
#> ..$ outlierTest : chr "IQR_test"
#> ..$ pCR.Update : logi TRUE
#> ..$ boundHandling: chr "reflect"
#> ..$ burnin.length: num 50
#> ..$ ndraw : num 500
#> ..$ maxtime : num Inf
#> ..$ Rthres : num 1.01
#> ..$ thin.t : logi NA
#> ..$ REPORT : num 0
#> ..$ parallel : chr "none"
#> ..$ ndim : int 5
#> ..$ DEpairs : num 2
#> ..$ nseq : int 5
#> ..$ Cb : logi NA
#> ..$ Wb : logi NA
#> $ in.burnin : logi TRUE
#> $ EXITFLAG : logi NA
#> $ AR : num [1:100, 1:2] NA 5 10 15 20 25 30 35 40 45 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : NULL
#> .. ..$ : chr [1:2] "fun.evals" "AR"
#> $ R.stat : num [1, 1:6] 5 -2 -2 -2 -2 -2
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : NULL
#> .. ..$ : chr [1:6] "fun.evals" "a" "b" "c" ...
#> $ CR : num [1:11, 1:4] 5 50 100 150 200 250 300 350 400 450 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : NULL
#> .. ..$ : chr [1:4] "fun.evals" "CR1" "CR2" "CR3"
#> $ outlier : int [1, 1:2] 250 2
#> $ EXITMSG : chr "Maximum function evaluations reached"
#> $ Sequences :List of 5
#> ..$ : 'mcmc' num [1:100, 1:5] 31.8 44.1 44.1 48.2 48.2 ...
#> .. ..- attr(*, "dimnames")=List of 2
#> .. .. ..$ : NULL
#> .. .. ..$ : chr [1:5] "a" "b" "c" "tw" ...
#> .. ..- attr(*, "mcpar")= num [1:3] 1 100 1
#> ..$ : 'mcmc' num [1:100, 1:5] 19.5 19.5 19.5 19.5 19.5 ...
#> .. ..- attr(*, "dimnames")=List of 2
#> .. .. ..$ : NULL
#> .. .. ..$ : chr [1:5] "a" "b" "c" "tw" ...
#> .. ..- attr(*, "mcpar")= num [1:3] 1 100 1
#> ..$ : 'mcmc' num [1:100, 1:5] 59.7 59.7 59.7 59.7 59.7 ...
#> .. ..- attr(*, "dimnames")=List of 2
#> .. .. ..$ : NULL
#> .. .. ..$ : chr [1:5] "a" "b" "c" "tw" ...
#> .. ..- attr(*, "mcpar")= num [1:3] 1 100 1
#> ..$ : 'mcmc' num [1:100, 1:5] 3.89 3.89 3.89 3.89 3.89 ...
#> .. ..- attr(*, "dimnames")=List of 2
#> .. .. ..$ : NULL
#> .. .. ..$ : chr [1:5] "a" "b" "c" "tw" ...
#> .. ..- attr(*, "mcpar")= num [1:3] 1 100 1
#> ..$ : 'mcmc' num [1:100, 1:5] 47.3 32.4 32.4 32.4 32.4 ...
#> .. ..- attr(*, "dimnames")=List of 2
#> .. .. ..$ : NULL
#> .. .. ..$ : chr [1:5] "a" "b" "c" "tw" ...
#> .. ..- attr(*, "mcpar")= num [1:3] 1 100 1
#> ..- attr(*, "class")= chr "mcmc.list"
#> $ X : num [1:5, 1:7] 9.06 6.72 7.61 8.03 4.76 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : NULL
#> .. ..$ : chr [1:7] "a" "b" "c" "tw" ...
#> $ hist.logp : num [1:100, 1:5] -230 -192 -192 -181 -181 ...
#> $ iterations: num 12
#> $ fun.evals : int 500
#> $ time : num 15.7
#> - attr(*, "class")= chr [1:2] "dream" "list"
#> [1] coef plot print simulate summary window
#> see '?methods' for accessing help and source code