GetEstimate_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 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
)

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.

Value

A list with the following components:

x	numeric matrix nsim * (length(par0)+length(SystError0)), MCMC simulations
fx	numeric vector, corresponding values f(x)

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)

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