HydroPortailStats (version 1.1.0)

GetEstimate_HBay: Bayesian estimation using historical data

Description

Bayesian estimation of a GEV or Gumbel distribution based on a mixed sample containing point (i.e. perfectly known) or interval (i.e. known to be within bounds) data. Systematic errors induced by rating curve errors can also be accounted for. Returns MCMC samples from the posterior distribution.

Usage

GetEstimate_HBay(
  y,
  dist,
  prior,
  SystErrorIndex = rep(0, NROW(y)),
  SystErrorPrior = list(),
  par0 = GetEstimate_ROUGH(0.5 * (y[, 1] + y[, 2])[is.finite(y[, 1] + y[, 2])], dist)$par,
  SystError0 = rep(1, length(SystErrorPrior)),
  mult = 0.1,
  eps = 0.1,
  batch.length = 100,
  batch.n = 100,
  moverate.min = 0.1,
  moverate.max = 0.5,
  mult.down = 0.9,
  mult.up = 1.1
)

Value

A list with the following components:

x

numeric matrix nsim * (length(par0)+length(SystError0)), MCMC simulations

fx

numeric vector, corresponding values f(x)

Arguments

y

numeric 2-column matrix, data. The first column gives the lower bound, the second column gives the upper bound. Where y[i,1]==y[i,2], the value is assumed perfectly known (up to systematic errors, see below). Where y[i,1]<y[i,2], the value is assumed to be in the interval [y[i,1];y[i,2]] -Inf and +Inf are allowed for data being only right- or left-censored (i.e. values known to be smaller than or larger than some threshold).

dist

character, distribution name. Only distributions 'GEV' and 'Gumbel' are supported.

prior

list of lists, prior distributions. For each parameter to be estimated, the prior is a list of the form pr=list(dist=..., par=...). See example below.

SystErrorIndex

integer vector, length NROW(y). Index of systematic errors. Rows where SystErrorIndex==k are all affected by the same multiplicative error gamma_k, typically induced by the kth rating curve. SystErrorIndex==0 means no systematic error. Should only contain integer values between 0 and N_{systematic errors}.

SystErrorPrior

list of lists, prior distribution for each systematic error. For instance for a systematic error in the range +/- 20%, you may use a Uniform between 0.8 and 1.2, or a triangular distribution with the same bounds and peak at 1.

par0

numeric vector, initial parameter guess.

SystError0

numeric vector, initial guess for systematic errors. Typically a vector of 1.

mult

numeric, initial jump standard deviations are set to mult * abs(par0)

eps

numeric, where par0 is zero, initial jump standard deviations are set to eps (to avoid jumps of size zero)

batch.length

integer, MCMC parameter: length of each non-adaptive batch

batch.n

integer, MCMC parameter: number of batches (= adaptation period). Total number of simulations is nsim=batch.n*batch.length

moverate.min

numeric in (0;1), MCMC parameter: lower bound for the desired move rate interval

moverate.max

numeric in (0;1), MCMC parameter: upper bound for the desired move rate interval

mult.down

numeric in (0;1), MCMC parameter: multiplication factor used to decrease jump size when move rate is too low.

mult.up

numeric (>1, avoid 1/mult.down), MCMC parameter: multiplication factor used to increase jump size when move rate is too high.

Examples

set.seed(98765)
n=50;n_censored=30
y0=Generate('GEV',c(100,50,-0.2),n)
y=cbind(y0,y0)
# Mimics censoring between 0 and 300 for first n_censored years
y[1:n_censored,1][y0[1:n_censored]<300]=0
y[1:n_censored,2][y0[1:n_censored]<300]=300
plot(y[,1]);points(y[,2])
# Systematic errors
SystErrorIndex=c(rep(1,n_censored),rep(2,n-n_censored))
SystErrorPrior=list(list(dist="Triangle",par=c(1,0.7,1.3)),
                    list(dist="Triangle",par=c(1,0.95,1.05)))
# Priors on GEV parameters
prior=list(list(dist="FlatPrior",par=NULL),
           list(dist="FlatPrior",par=NULL),
           list(dist="Normal",par=c(0,0.25)))
# Go!
mcmc=GetEstimate_HBay(y=y,dist='GEV',prior=prior,
                      SystErrorIndex=SystErrorIndex,
                      SystErrorPrior=SystErrorPrior,
                      # The values below aim at making this example fast to run.
                      # In practice, it is recommended to use the default values
                      # (batch.length=100,batch.n=100) or larger.
                      batch.length=25,batch.n=20) 
graphicalpar=par(mfrow=c(2,3))
for(i in 1:5){hist(mcmc$x[,i])}
par(graphicalpar)