Fit a hydromad model using general optimisation algorithms.

Fits a hydromad model using R's optim or nlminb functions. Has multi-start and pre-sampling options.

fitByOptim(
  MODEL,
  objective = hydromad.getOption("objective"),
  method = hydromad.getOption("optim.method"),
  control = list(),
  samples = hydromad.getOption("fit.samples"),
  sampletype = c("latin.hypercube", "random", "all.combinations"),
  initpars = NULL,
  multistart = FALSE,
  vcov = FALSE,
  hessian = vcov
)

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.
objective: objective function to maximise, given as a function(Q, X, ...). See objFunVal.
method, control: optimisation algorithm and settings. See optim. The default is "Nelder-Mead" (a simplex algorithm). An additional method "PORT" can be given, which calls nlminb.
samples: number of parameter sets to test.
sampletype: sampling scheme -- see parameterSets.
initpars: initial parameter set.
multistart: if this is TRUE, then each of the initial parameter samples (based on samples and sampletype) are used as a starting value for a separate optimisation run, and the best result is kept.
vcov, hessian: if vcov = TRUE, the parameter variance-covariance matrix will be estimated as the inverse of the hessian matrix returned by optim(). This may be a very poor estimate! It can be extract using vcov. Does not work with method = "PORT".

Value

the best model from those sampled, according to the given objective function. Also, these extra elements are inserted:

fit.result: the result from optim.
objective: the objective function used.
funevals: total number of evaluations of the model simulation function.
timing: timing vector as returned by system.time.

Details

See optim for a brief description of the available algorithms, including references to literature.

fitByOptim1 handles a single free parameter using optimize; it is called by fitByOptim, with a warning, in that case.

Author

Felix Andrews felix@nfrac.org

Examples


data(Cotter)
x <- Cotter[1:1000]

## IHACRES CWI model with power law unit hydrograph
modx <- hydromad(x, sma = "cwi", routing = "powuh")
modx
#> 
#> Hydromad model with "cwi" SMA and "powuh" routing:
#> Start = 1966-05-01, End = 1969-01-24
#> 
#> SMA Parameters:
#>       lower upper     
#> tw        0   100     
#> f         0     8     
#> scale    NA    NA     
#> l         0     0 (==)
#> p         1     1 (==)
#> t_ref    20    20 (==)
#> Routing Parameters:
#>   lower upper  
#> a  0.01    60  
#> b  0.50     3  
#> c  0.50     2  

set.seed(0)
foo <- fitByOptim(modx, samples = 100)

summary(foo)
#> 
#> Call:
#> hydromad(DATA = x, sma = "cwi", routing = "powuh", a = 7.87367, 
#>     b = 1.41609, c = 2, tw = 100, f = 4.77557, scale = 0.00128754)
#> 
#> Time steps: 900 (0 missing).
#> Runoff ratio (Q/P): (0.7028 / 2.285) = 0.3075
#> rel bias: -2.036e-16
#> r squared: 0.7508
#> r sq sqrt: 0.8364
#> r sq log: 0.8341
#> 
#> For definitions see ?hydromad.stats
#> 

## return value from optim (with 'objseq' added):
str(foo$fit.result)
#> List of 7
#>  $ par        : Named num [1:5] 7.87 1.42 2 100 4.78
#>   ..- attr(*, "names")= chr [1:5] "a" "b" "c" "tw" ...
#>  $ objective  : num -0.855
#>  $ convergence: int 0
#>  $ iterations : int 25
#>  $ evaluations: Named int [1:2] 26 149
#>   ..- attr(*, "names")= chr [1:2] "function" "gradient"
#>  $ message    : chr "relative convergence (4)"
#>  $ objseq     : num [1:276] 0.705 0.773 0.782 -1.302 0.706 ...

## plot objective function value convergence over time
xyplot(optimtrace(foo),
  type = "b",
  xlab = "function evaluations", ylab = "objective fn. value"
)


## repeat optimisation with single random starting points
fooreps <-
  replicate(4,
    fitByOptim(modx,
      samples = 1, sampletype = "random",
      objective = hmadstat("r.squared"),
      method = "Nelder-Mead"
    ),
    simplify = FALSE
  )
#> Warning: optim() reached maximum iterations
#> Warning: optim() reached maximum iterations
#> Warning: optim() reached maximum iterations
#> Warning: optim() reached maximum iterations
names(fooreps) <- paste("rep.", seq_along(fooreps))
## extract and plot the optimisation traces
traces <- lapply(fooreps, optimtrace)
tracesraw <- lapply(fooreps, optimtrace, raw = TRUE)
xyplot(do.call("merge", traces),
  superpose = TRUE,
  sub = 'method = "Nelder-Mead"',
  xlab = "Fn. evaluations", ylab = "Objective value",
  auto.key = list(corner = c(1, 0))
) +
  xyplot(do.call("merge", tracesraw),
    superpose = TRUE,
    type = "p", cex = 0.5
  )


## if you try it again with method = "PORT" you will find that all
## replicates converge to the optimum regardless of starting point.