defineFeasibleSet {hydromad}  R Documentation 
Extract the feasible (or behavioural) parameter sets meeting some criteria. These could be used as a feasible set or to estimate prediction quantiles according to the GLUE (Generalised Likelihood Uncertainty Estimation) method.
defineFeasibleSet(x, ...) ## S3 method for class 'hydromad' defineFeasibleSet(x, ..., thin = NA) ## Default S3 method: defineFeasibleSet(x, model, objseq = rep(1, NROW(x)), frac.within = 0, within.rel = 0.01, within.abs = 0.1, groups = NULL, FUN = sum, target.coverage = 1, threshold = Inf, glue.quantiles = c(0, 1), ...)
x 
in the 
... 
extra arguments to the 
thin 
interval between samples for results from DREAM.
As it is a Markov Chain Monte Carlo method, the
sequences should be thinned first to remove autocorrelation and
achieve an efficient sample of the parameter distributions. The
default thinning interval 
model 
the 
objseq 
objective function values corresponding to the parameter sets

frac.within, within.abs, within.rel 
model simulations are only retained in the feasible set if some
fraction 
groups, FUN 

target.coverage 
fraction of the observed values to be contained within the overall ranges
of simulated values (minimum and maximum on each time step) from the
feasible set of parameters. Note, this does not refer to the

threshold 
value of the objective function (the objective function used to
generate 
glue.quantiles 
if specified, these GLUE quantiles of the ensemble simulations will be
calculated and stored. They can be extracted with 
a modified version of model
, with added elements feasible.set
,
feasible.scores
, feasible.fitted
, glue.quantiles
and feasible.threshold
.
Can be passed to xyplot
, fitted
, predict
,
update
, coef
, and print
.
Felix Andrews felix@nfrac.org
predict.hydromad
,
update.hydromad
,
fitBySampling
,
fitByDream
data(Queanbeyan) ts74 < window(Queanbeyan, start = "19740101", end = "19761201") mod < hydromad(ts74, routing = "expuh", rfit = list("inverse", order = c(2,1))) mod < update(mod, sma = "cwi", tw = c(0,100), f = c(0, 8), loss = c(0.1, 0.1)) ## Calculate the set of simulations within 15% error (or 1 mm/day) 90% of time. ## In this case we do not need to calculate objective function values ## beforehand. For GLUE quantiles, however, need to give 'objective'. psets < parameterSets(coef(mod), samples = 300) feas < defineFeasibleSet(psets, model = mod, frac.within = 0.9, within.rel = 0.15, within.abs = 1) ## How many of the 300 possible parameter sets were retained? nrow(coef(feas, feasible.set = TRUE)) ## View ranges of parameters in feasible set feas ## Plot simulation bounds xyplot(feas, feasible.bounds = TRUE, cut = 3) ## Generate set of simulations with NSE > 0.5, for GLUE. ## First, need to calculate objective function values: fit < fitBySampling(mod, samples = 300, objective = hmadstat("r.squared")) ## Calculate 5 percent and 95 percent GLUE quantiles (i.e. weighted). fitglu < defineFeasibleSet(fit, threshold = 0.5, glue.quantiles = c(0.05, 0.95)) ## Coverage of the GLUE quantile simulations (within 0.1 mm/day) sim < fitted(fitglu, feasible.bounds = TRUE) head(sim) mean((sim[,1] < observed(fitglu) + 0.1) & (sim[,2] > observed(fitglu)  0.1)) ## Or  keep adding parameter sets until we reach a target coverage: ## Calculate 5 percent and 95 percent GLUE quantiles (i.e. weighted). fitglu < defineFeasibleSet(fit, target.coverage = 0.9, glue.quantiles = c(0.05, 0.95)) ## Coverage of the GLUE quantile simulations (within 0.1 mm/day) ## (not to be confused with the target.coverage to define overall feasible set) sim < fitted(fitglu, feasible.bounds = TRUE) mean((sim[,1] < observed(fitglu) + 0.1) & (sim[,2] > observed(fitglu)  0.1)) ## Plot simulated GLUE quantiles xyplot(fitglu, feasible.bounds = TRUE, cut = 3) xyplot(fitglu, feasible.bounds = TRUE, cut = 3, scales = list(y = list(log = TRUE))) ## Summarise size of the simulation bounds: lower as fraction of upper summary(coredata(sim[,1] / sim[,2])) ## Simulate on a new data period newglu < update(fitglu, newdata = window(Queanbeyan, start = "19800101", end = "19820101"), glue.quantiles = c(0.05, 0.95)) ## The new period is very dry, all model simulations overestimate flow. xyplot(newglu, feasible.bounds = TRUE, cut = 3) ## Coverage of the GLUE quantile simulations (within 0.1 mm/day) sim < fitted(newglu, feasible.bounds = TRUE) mean((sim[,1] < observed(newglu) + 0.1) & (sim[,2] > observed(newglu)  0.1))