R: Extract the feasible parameter sets meeting some criteria.

defineFeasibleSet {hydromad}

R Documentation

Extract the feasible parameter sets meeting some criteria.

Description

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.

Usage

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),
                  ...)

Arguments

`x`	in the `hydromad` method this is a `hydromad` model object which has been run through either `fitBySampling` or `fitByDream` to generate a large number of random parameter sets with associated objective function scores. In the `default` method, the first argument (`x`) is a matrix of parameter values corresponding to the objective function values `objseq`.
`...`	extra arguments to the `hydromad` method are passed on to the `default` method. Extra arguments to the `default` method will result in an error.
`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 `NA` uses the empirical autocorrelation. See `window.mcmc`.
`model`	the `hydromad` model object to be used to run simulations. Does not apply to the `hydromad` method.
`objseq`	objective function values corresponding to the parameter sets `x`. This does not apply to the `hydromad` method, where it is assumed that the objective function values have already been calculated. Note that `objseq` can be omitted if not using the `target.coverage` or `threshold` arguments (i.e. if just defining an error criterion with `frac.within, within.rel, within.abs`), so the original simulation run is not necessary.
`frac.within, within.abs, within.rel`	model simulations are only retained in the feasible set if some fraction `frac.within` of the simulated values are within a fraction `within.rel` of the observed values OR within an absolute difference of `within.abs` (typically mm/day).
`groups, FUN`	`groups` is an optional grouping variable, of the same length as the observed data in `model`, used to aggregate the observed and fitted time series. The function `FUN` is applied to each group. In this case, the error criteria `frac.within, within.rel, within.abs` are evaluated on the aggregated values, not the raw time series. Also `target.coverage` applies to the aggregated values. Typically `groups` would be generated by `cut.Date` (for regular time periods) or `eventseq` (for events). Note that the `feasible.bounds` (e.g. for plotting) in this case will be an aggregated time series, it will not correspond to the original time index. Use `update` or `predict` to generate bounds on the original time index.
`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 `glue.quantiles` values, but to the overall maximum and minimum. The simulated values can be within the tolerance limit of observed values given by `within.abs` to count as coverage.
`threshold`	value of the objective function (the objective function used to generate `objseq`) used to define the feasible set: all parameter sets above this threshold value will be kept. Also this `threshold` value is subtracted from the objective function values to calculate weights when `glue.quantiles` is given. If left as `-Inf`, it will be set to the minimum objective function value in the final feasible set.
`glue.quantiles`	if specified, these GLUE quantiles of the ensemble simulations will be calculated and stored. They can be extracted with `fitted` and shown in `xyplot`. The quantiles are calcualted using the `wtd.quantile` function, and weighted according to `objseq` - `threshold`. If `glue.quantiles` is left as `c(0,1)` then the overall minimum and maximum simulated values are taken on each time step, which is faster.

Value

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.

Author(s)

Felix Andrews felix@nfrac.org

Examples

data(Queanbeyan)
ts74 <- window(Queanbeyan, start = "1974-01-01", end = "1976-12-01")
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 = "1980-01-01", end = "1982-01-01"),
                 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))

[Package hydromad version 0.9-18 Index]