RBaM (version 1.1.1)

SPD_estimate: Estimation of a BaRatin-SPD model

Description

Run BaM to estimate a BaRatin-SPD model. The following assumptions hold:

• Only parameters b/k's and a's can be variable. Exponents c's are stable.

• Incremental changes affecting parameters b/k's are additive, while incremental changes affecting parameters a's are multiplicative.

• The prior distribution for parameters b/k's and associated incremental changes has to be 'Gaussian'.

• The prior distribution for parameters a's and associated incremental changes has to be 'LogNormal'.

Usage

SPD_estimate(
  workspace,
  controlMatrix,
  pars,
  bVAR,
  aVAR,
  deltaPars,
  periods,
  H,
  Q,
  uQ = 0 * Q,
  nPeriods = lapply(periods, max),
  BaRatinFlavor = "BaRatinBAC",
  remnant = remnantErrorModel(funk = "Linear", par = list(parameter("g1", 1, "LogNormal",
    c(0, 10)), parameter("g2", 0.1, "LogNormal", c(log(0.1), 10)))),
  mcmcOpt = mcmcOptions(),
  mcmcCook = mcmcCooking()
)

Value

A data frame containing the MCMC simulations performed by BaM.

Arguments

workspace

Character, directory where config and result files are stored.

controlMatrix

Integer matrix, control matrix, dimension nControl*nControl.

pars

list of parameter objects, parameters of the model. For VAR parameters, they will be interpreted as the prior for period 1.

bVAR

Logical vector, size nControl. bVAR[i]=TRUE means that the b/k parameter of control i is variable, otherwise it is stable.

aVAR

Logical vector, size nControl. aVaR[i]=TRUE means that the a parameter of control i is variable, otherwise it is stable.

deltaPars

list, prior parameters of incremental changes for each VAR parameter. deltaPars should be a named list, with the names corresponding to the names of the parameters that have been declared variable. Each element of deltaPars is then a numeric vector of size 2. For b/k's, the 2 values are the mean/sd of the Gaussian prior for ADDITIVE incremental changes. For a's, the 2 parameters are the meanlog/sdlog of the LogNormal prior for MULTIPLICATIVE incremental changes.

periods

list, period index for each VAR parameter. periods should be a named list as previously. Each element of the list is an integer vector (starting at 1) with same length as the calibration data. Periods do not need to be the same for all VAR parameters.

H

numeric vector, gauging stages.

Q

numeric vector, gauging discharges.

uQ

numeric vector, gauging discharge uncertainties.

nPeriods

list, number of periods for each VAR parameter. nPeriods should be a named list as deltaPars and periods. In general and by default, nPeriods[[i]] is just the max of periods[[i]], but this is not compulsory: there could be one or several additional periods with no gaugings.

BaRatinFlavor

character, either 'BaRatinBAC' (default) or 'BaRatin' (the original k-a-c parameterization). It is in general easier to specify priors on changes affecting b's than k's, hence the default choice. However, 'BaRatinBAC' requires some numerical resolution and it is hence a bit slower and more prone to failures than 'BaRatin'.

remnant

remnantErrorModel object, by default the structural standard deviation varies as an affine function of simulated discharges, with very wide priors on coefficients g1 and g2.

mcmcOpt

mcmcOptions object, MCMC options passed to BaM.

mcmcCook

mcmcCooking object, MCMC cooking options (burn and slice) passed to BaM.

Examples

# Calibration data
H=MeyrasGaugings$h
Q=MeyrasGaugings$Q
uQ=MeyrasGaugings$uQ
# Control matrix
controlMatrix=rbind(c(1,0,0),c(0,1,0),c(0,1,1))
# Declare variable parameters.
bVAR=c(TRUE, TRUE, FALSE) # b's for first 2 controls (k1 and k2) are VAR
aVAR=c(TRUE, FALSE, FALSE) # a for first control (a1) is VAR
# Define priors.
b1=parameter(name='b1',init=-0.6,prior.dist='Gaussian',prior.par=c(-0.6,0.5))
a1=parameter(name='a1',init=exp(2.65),prior.dist='LogNormal',prior.par=c(2.65,0.35))
c1=parameter(name='c1',init=1.5,prior.dist='Gaussian',prior.par=c(1.5,0.025))
b2=parameter(name='b2',init=0,prior.dist='Gaussian',prior.par=c(-0.6,0.5))
a2=parameter(name='a2',init=exp(3.28),prior.dist='LogNormal',prior.par=c(3.28,0.33))
c2=parameter(name='c2',init=1.67,prior.dist='Gaussian',prior.par=c(1.67,0.025))
b3=parameter(name='b3',init=1.2,prior.dist='Gaussian',prior.par=c(1.2,0.2))
a3=parameter(name='a3',init=exp(3.48),prior.dist='LogNormal',prior.par=c(3.46,0.38))
c3=parameter(name='c3',init=1.67,prior.dist='Gaussian',prior.par=c(1.67,0.025))
pars=list(b1,a1,c1,b2,a2,c2,b3,a3,c3)
# Define properties of VAR parameters.
deltaPars=list(b1=c(0,0.25),a1=c(0,0.2),b2=c(0,0.5))
periods=list(b1=MeyrasGaugings$Period,a1=c(rep(1,49),rep(2,55)),b2=MeyrasGaugings$Period)
# Run BaM and estimate SPD parameters
mcmcOpt=mcmcOptions(nAdapt=20,nCycles=25) # only few iterations so that the example runs fast.
mcmc=SPD_estimate(workspace=tempdir(),controlMatrix=controlMatrix,pars=pars,
                  bVAR=bVAR,aVAR=aVAR,deltaPars=deltaPars,periods=periods,
                  H=H,Q=Q,uQ=uQ,mcmcOpt=mcmcOpt)