library(airGR) library(coda) library(FME) library(ggmcmc) set.seed(123) load(system.file("vignettesData/vignetteParamMCMC.rda", package = "airGR"))

In this vignette, we give an example of parameter estimation within a Bayesian MCMC approach.

We use the **GR4J** model and we assume that the R global environment contains data and functions from the **airGR** Get Started vignette.

example("Calibration_Michel")

In order for the `RunModel_*()`

functions to run faster during the parameter estimation process, it is recommended that the outputs contain only the simulated flows (see the `Outputs_Sim`

argument in the `CreateRunOptions()`

help page).

RunOptions <- airGR::CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel, IndPeriod_Run = Ind_Run, Outputs_Sim = "Qsim")

Please refer to the 2.1 Plugging in new calibration vignette for explanations on how to plug in a parameter estimation algorithm to **airGR**.

Please note that this vignette is only for illustration purposes and does not provide any guidance about which parameter inference strategy is recommended for the family of the GR models.

We show how to use the DRAM algorithm for SLS Bayesian inference, with the `modMCMC()`

function of the FME package.
First, we need to define a function that returns twice the opposite of the log-likelihood for a given parameter set.

Nota: in the `LogLikeGR4J()`

function, the computation of the log-likelihood is simplified in order to ensure a good computing performance. It corresponds to a translation of the two following lines.

Likelihood <- sum((ObsY - ModY)^2, na.rm = TRUE)^(-sum(!is.na(ObsY)) / 2) LogLike <- -2 * log(Likelihood)

In our simplified setting of Gaussian likelihood with measurement error integrated out, the log of the sum of squared error is related to the log-likelihood.

Note that we do not use here any discharge transformation, although applying Box-Cox transformation is quite popular in hydrological modelling.

LogLikeGR4J <- function(ParamOptim) { ## Transformation to real space RawParamOptim <- airGR::TransfoParam_GR4J(ParamIn = ParamOptim, Direction = "TR") ## Simulation given a parameter set OutputsModel <- airGR::RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions, Param = RawParamOptim) ## Computation of the log-likelihood: N * log(SS) ObsY <- InputsCrit$Obs ModY <- OutputsModel$Qsim LogLike <- sum(!is.na(ObsY)) * log(sum((ObsY - ModY)^2, na.rm = TRUE)) }

We start by using the PORT optimization routine to estimate the best-fit parameters.

optPORT <- stats::nlminb(start = c(4.1, 3.9, -0.9, -8.7), objective = LogLikeGR4J, lower = rep(-9.9, times = 4), upper = rep(9.9, times = 4), control = list(trace = 1)) iniParPORT <- optPORT$par

We run 3 chains with different initial values to assess the convergence of the Markov chains. The number of iterations is fixed to 2000 with a burning length of 0.

Nota: in this example, there are relatively few iterations (2000), in order to limit the running time of this vignette. In addition, the burning length has been set to zero in order to show the convergence process but, in a true exercise, it is better to define more iterations (5000) and to burn the first iterations.

With the DRAM algorithm, the covariance of the proposal is updated every 100 runs and delayed rejection is applied.

iniParPORT <- data.frame(Chain1 = iniParPORT, Chain2 = iniParPORT, Chain3 = iniParPORT, row.names = paste0("X", 1:4)) iniParPORT <- sweep(iniParPORT, MARGIN = 2, STATS = c(1, 0.9, 1.1), FUN = "*") iniParPORT[iniParPORT < -9.9] <- -9.9 iniParPORT[iniParPORT > +9.9] <- +9.9 mcmcDRAM <- apply(iniParPORT, MARGIN = 2, FUN = function(iIniParPORT) { FME::modMCMC(f = LogLikeGR4J, p = iIniParPORT, lower = rep(-9.9, times = 4), ## lower bounds for GR4J upper = rep(+9.9, times = 4), ## upper bounds for GR4J niter = 2000, jump = 0.01, outputlength = 2000, burninlength = 0, updatecov = 100, ## adaptative Metropolis (AM) ntrydr = 2) ## delayed rejection (RD) })

There are several diagnostics that can be used to check the convergence of the chains. The R package coda provides several diagnostic tests.

Among others, the Gelman and Rubin's convergence can be used. A value close to 1 suggests acceptable convergence.
The result will be better with more iterations than 2000. As we kept the iterations during the convergence process, we have to set the `autoburnin`

argument to `TRUE`

in order to consider only the second half of the series.

Note that we rescale model parameter sets of the GR4J model from the transformed space to the real space.

multDRAM <- coda::as.mcmc.list(lapply(mcmcDRAM, FUN = function(x) { coda::as.mcmc(airGR::TransfoParam_GR4J(as.matrix(x$pars), Direction = "TR")) })) gelRub <- coda::gelman.diag(multDRAM, autoburnin = TRUE)$mpsrf

gelRub

In addition, graphical tools can be used, with for example the ggmcmc package.

First, the evolution of the Markov chains can be seen with a traceplot:

parDRAM <- ggmcmc::ggs(multDRAM) ## to convert object for using by all ggs_* graphical functions ggmcmc::ggs_traceplot(parDRAM)

The posterior density for each parameter can then be visualised:

burnParDRAM <- parDRAM[parDRAM$Iteration > 1000, ] # to keep only the second half of the series ggmcmc::ggs_density(burnParDRAM)

Finally, a paired plot can be useful to look at the correlation between parameters. Highly correlated parameters can be seen as an indication for a structural deficit of the model used.

ggmcmc::ggs_pairs(burnParDRAM, lower = list(continuous = "density"))

We only presented one MCMC algorithm and one parameter inference setting. Nonetheless, other approaches can be explored within the same framework.

One can for example use different data transformations to deal with the limitations of the Gaussian error model.

**Any scripts or data that you put into this service are public.**

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.