Hydro3_HBay: Bayesian estimation using historical data

In HydroPortailStats: 'HydroPortail' Statistical Functions

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 an Hydro3 object

Usage

Hydro3_HBay(
  y,
  dist,
  prior = GetDefaultPrior(GetParNumber(dist)),
  SystErrorIndex = rep(0, NROW(y)),
  SystErrorPrior = list(),
  options = options_def,
  mcmcoptions = mcmcoptions_def,
  do.KS = TRUE,
  do.MK = TRUE,
  do.Pettitt = TRUE
)

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.
options	list, options, see details below.
mcmcoptions	list, MCMC options, see details below.
do.KS, do.MK, do.Pettitt	logical, perform KS/MK/Pettitt tests?

Details

The argument 'options' allows controlling various properties of the analysis and results. It is a list with the following components:

• FreqFormula, character, formula for computing nonexceedance frequency, see ?GetEmpFreq.

• pgrid, numeric vector, probabilities defining the x values where pdf f(x) and cdf F(x) are computed. These x values are quantiles from the estimated distribution with probabilities pgrid.

• Tgrid, numeric vector, return periods where quantile function q(T) is computed.

• IClevel, numeric, level of uncertainty interval.

• p2T, numeric, conversion factor between nonexceedance probability p and return period T. p=1-1/(p2T*T). Here p2T=1 in general since GEV/Gumbel are applied to annual maxima in general.

• invertT, logical, when invertT=TRUE, LARGE return periods correspond to SMALL data values. This is typically used for low-flow statistics. Unused here.

• splitZeros, logical, when splitZeros=TRUE zero and negative values are removed from the data y before estimating the distribution,and are used to estimate the probability of zeros p0. This is typically used for low-flow statistics to estimate the probability of zero streamflow. Unused here.

• lang, chanracter, language ('fr' or 'en').

• nsim, integer, number of replicated parameters representing uncertainty. Unused here (derives from mcmc options)

The argument 'mcmcoptions' is a list controlling MCMC properties:

• mult, numeric, see ?Metropolis_OAAT_adaptive

• eps, numeric, see ?Metropolis_OAAT_adaptive

• batch.length, integer, see ?Metropolis_OAAT_adaptive

• batch.n, integer, see ?Metropolis_OAAT_adaptive

• moverate.min, numeric, see ?Metropolis_OAAT_adaptive

• moverate.max, numeric, see ?Metropolis_OAAT_adaptive

• mult.down, numeric, see ?Metropolis_OAAT_adaptive

• mult.up, numeric, see ?Metropolis_OAAT_adaptive

• burn, numeric, burn-in factor, e.g. if burn=0.2 the first 20 percents of MCMC samples are discarded

• slim, integer, sliming factor, e.g. if slim=5 only one MCMC sample every 5 is kept (after burn-in)

Value

A list with the following components:

dist	character, estimated distribution.
ok	logical, did estimation succeed?
err	integer, error code (0 if ok).
message	error message.
empirical	data frame, sorted data and empirical estimates (nonexceedance frequency, return period and reduced variate). NOTE: interval data are replaced by a value randomly sampled from a GEV constrained in this interval. See ?GenerateWithinBounds.
pcdf	data frame, estimated pdf and cdf
quantile	data frame, estimated quantiles and uncertainty intervals
par	data frame, estimated parameters and uncertainty intervals
KS	list, result of the Kolmogorov-Smirnov test, see ?KS. NOTE: interval data are replaced by a value randomly sampled from a GEV constrained in this interval. See ?HBay_simGEV.
MK	list, result of the Mann-Kendall test, see ?MK. Same note as KS test.
Pettitt	list, result of the Pettitt test, see ?Pettitt. Same note as KS test.
u	list, parameter uncertainty in the form of a covariance matrix ($cov) and simulated parameter replicates ($sim). Also contains error-handling flags $ok, $err and $message.

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)))
# Handle MCMC options
# 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.
mcmcoptions=mcmcoptions_def
mcmcoptions$batch.length=25
mcmcoptions$batch.n=20
# Go!
H3=Hydro3_HBay(y=y,dist='GEV',prior=prior,
               SystErrorIndex=SystErrorIndex,
               SystErrorPrior=SystErrorPrior,
               mcmcoptions=mcmcoptions) 
Hydro3_Plot(H3)

HydroPortailStats documentation built on Sept. 12, 2024, 9:36 a.m.