# \donttest{
# use of \donttest{} because execution time exceeds 5 seconds
data(sim_wrc_hcc)
# define number of cores for parallel computations
if(interactive()) ncpu <- parallel::detectCores() - 1L else ncpu <- 1L
# estimate parameters for single sample ...
# ... from wrc and hcc data
plot(rfit_wrc_hcc <- fit_wrc_hcc(
wrc_formula = wc ~ head, hcc_formula = hc ~ head,
data = sim_wrc_hcc, subset = id == 1))
coef(rfit_wrc_hcc, gof = TRUE)
# ... from wrc data
plot(rfit_wrc <- fit_wrc_hcc(
wrc_formula = wc ~ head,
data = sim_wrc_hcc, subset = id == 1))
coef(rfit_wrc, gof = TRUE)
# ... from hcc data
plot(rfit_hcc <- fit_wrc_hcc(
hcc_formula = hc ~ head,
data = sim_wrc_hcc, subset = id == 1))
coef(rfit_hcc, gof = TRUE)
# ... from wrc and hcc data
# keeping some parameters fixed
plot(rfit_wrc_hcc_fixed <- fit_wrc_hcc(
wrc_formula = wc ~ head, hcc_formula = hc ~ head,
data = sim_wrc_hcc, subset = id == 1,
param = c(alpha = 2.1, thetar = 0.1),
fit_param = default_fit_param(alpha = FALSE, thetar = FALSE)),
y = rfit_wrc_hcc)
coef(rfit_wrc_hcc, gof = TRUE)
coef(rfit_wrc_hcc_fixed, gof = TRUE)
# estimate parameters for 3 samples simultaneously by ...
# ... unconstrained, global optimisation algorithm NLOPT_GN_MLSL (default)
rfit_uglob <- fit_wrc_hcc(
wrc_formula = wc ~ head | id, hcc_formula = hc ~ head | id,
data = sim_wrc_hcc,
control = control_fit_wrc_hcc(pcmp = control_pcmp(ncores = ncpu)))
summary(rfit_uglob)
op <- par(mfrow = c(3, 2))
plot(rfit_uglob)
on.exit(par(op))
# ... unconstrained, global optimisation algorithm SCEoptim
rfit_sce <- update(
rfit_uglob,
control = control_fit_wrc_hcc(
settings = "sce", pcmp = control_pcmp(ncores = ncpu)))
coef(rfit_sce, gof = TRUE, lc = TRUE)
convergence_message(2, sce = TRUE)
op <- par(mfrow = c(3, 2))
plot(rfit_sce, y = rfit_uglob)
on.exit(par(op))
# ... unconstrained, local optimisation algorithm NLOPT_LN_BOBYQA,
# logging iteration results to console
rfit_uloc <- update(
rfit_uglob,
param = as.matrix(coef(rfit_uglob, what = "nonlinear")),
control = control_fit_wrc_hcc(
settings = "ulocal", pcmp = control_pcmp(ncores = 1L)),
verbose = 2)
coef(rfit_uloc, gof = TRUE, lc = TRUE)
# ... constrained, global optimisation algorithm NLOPT_GN_ISRES
rfit_cglob <- update(
rfit_uglob,
ratio_lc_lt_bound = c(lower = 0.8, upper = 1.2),
control = control_fit_wrc_hcc(
settings = "cglobal", nloptr = control_nloptr(ranseed = 1),
pcmp = control_pcmp(ncores = ncpu)))
coef(rfit_cglob, gof = TRUE, lc = TRUE)
# ... constrained, local optimisation algorithm NLOPT_LD_CCSAQ
# starting from unconstrained, locally fitted initial values
rfit_cloc_1 <- update(
rfit_uglob,
param = coef(rfit_uloc, what = "nonlinear"),
ratio_lc_lt_bound = c(lower = 0.8, upper = 1.2),
control = control_fit_wrc_hcc(
settings = "clocal", pcmp = control_pcmp(ncores = ncpu)))
coef(rfit_cloc_1, gof = TRUE, lc = TRUE)
op <- par(mfrow = c(3, 2))
plot(x = rfit_uloc, y = rfit_cloc_1)
on.exit(par(op))
# ... constrained, local optimisation algorithm NLOPT_LD_CCSAQ
# starting from constrained, globally fitted initial values,
# using sample-specific evaporation rates and limits for ratio lc/lt
rfit_cloc_2 <- update(
rfit_uglob,
param = as.matrix(coef(rfit_cglob, what = "nonlinear")),
e0 = c("1" = 0.002, "2" = 0.0025, "3" = 0.003),
ratio_lc_lt_bound = rbind(
"1" = c(lower = 0.7, upper = 2),
"2" = c(lower = 0.8, upper = 1.4),
"3" = c(lower = 0.8, upper = 2)
),
control = control_fit_wrc_hcc(
settings = "clocal", pcmp = control_pcmp(ncores = ncpu)))
coef(rfit_cloc_2, gof = TRUE, lc = TRUE, e0 = TRUE)
# ... global optimisation algorithm NLOPT_GN_MLSL
# with sample-specific box-constraints
rfit_uglob_bc <- update(
rfit_uglob,
lower_param = rbind(
"1" = c(alpha = 2.4, n = 1.11, thetar = 0.2, thetas = 0.6),
"2" = c(alpha = 1.2, n = 1.12, thetar = 0.2, thetas = 0.6),
"3" = c(alpha = 1.2, n = 1.13, thetar = 0.2, thetas = 0.6)
),
upper_param = rbind(
"1" = c(alpha = 20.1, n = 2.51, thetar = 0.6, thetas = 0.6),
"2" = c(alpha = 20.2, n = 1.5, thetar = 0.6, thetas = 0.6),
"3" = c(alpha = 1.3, n = 2.53, thetar = 0.6, thetas = 0.6)
)
)
coef(rfit_uglob, gof = TRUE)
coef(rfit_uglob_bc, gof = TRUE)
# }Run the code above in your browser using DataLab