Commit 91d6c8ec authored by Delaigue Olivier's avatar Delaigue Olivier
Browse files

v1.1.2.18 CLEAN: syntaxe revision of vignette param_mcmc

parent 5f79785c
Package: airGR
Type: Package
Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
Version: 1.1.2.17
Version: 1.1.2.18
Date: 2018-10-23
Authors@R: c(
person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
......
......@@ -13,7 +13,7 @@ output:
### 1.1.2.17 Release Notes (2018-10-23)
### 1.1.2.18 Release Notes (2018-10-23)
......
......@@ -18,7 +18,7 @@ library(ggmcmc)
library(dplyr)
# source("airGR.R")
set.seed(123)
load(system.file("vignettes_data/Vignette_Param.rda", package = "airGR"))
load(system.file("vignettesData/vignetteParamMCMC.rda", package = "airGR"))
```
......@@ -31,10 +31,10 @@ In this vignette, we give an example of parameter estimation within a Bayesian M
We use the **GR4J** model and we assume that the R global environment contains data and functions from the **airGR** [Get Started](V01_get_started.html) vignette.
```{r, warning=FALSE, fig.keep='none', results='hide', fig.height=10, fig.width=10, eval=TRUE, echo=FALSE, message=FALSE}
example("Calibration_Michel", echo = FALSE, ask = FALSE)
```
```{r, echo=TRUE, eval=FALSE}
<!-- ```{r, warning=FALSE, fig.keep='none', results='hide', fig.height=10, fig.width=10, eval=TRUE, echo=FALSE, message=FALSE} -->
<!-- example("Calibration_Michel", echo = FALSE, ask = FALSE) -->
<!-- ``` -->
```{r, echo=TRUE, eval=FALSE, eval=FALSE}
example("Calibration_Michel")
```
......@@ -48,29 +48,29 @@ Please note that this vignette is only for illustration purposes and does not pr
We show how to use the DRAM algorithm for SLS Bayesian inference, with the `modMCMC()` function of the [FME](https://cran.r-project.org/package=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 `RunAirGR4J()` 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.
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.
```{r, echo=TRUE, eval=FALSE}
Likelihood <- sum((ObsY - ModY)^2, na.rm = TRUE)^(-sum(!is.na(ObsY)) / 2)
LL <- -2 * log(Likelihood)
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.
```{r, results='hide'}
RunAirGR4J <- function(Param_Optim) {
```{r, results='hide', eval=FALSE}
LogLikeGR4J <- function(ParamOptim) {
## Transformation to real space
Param_Optim_Vre <- airGR::TransfoParam_GR4J(ParamIn = Param_Optim,
RawParamOptim <- airGR::TransfoParam_GR4J(ParamIn = ParamOptim,
Direction = "TR")
## Simulation given a parameter set
OutputsModel <- airGR::RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions,
Param = Param_Optim_Vre)
Param = RawParamOptim)
## Computation of the log-likelihood: N * log(SS)
ObsY <- InputsCrit$Qobs
ModY <- OutputsModel$Qsim
LL <- sum(!is.na(ObsY)) * log(sum((ObsY - ModY)^2, na.rm = TRUE))
LogLike <- sum(!is.na(ObsY)) * log(sum((ObsY - ModY)^2, na.rm = TRUE))
}
```
......@@ -80,10 +80,10 @@ RunAirGR4J <- function(Param_Optim) {
We start by using the PORT optimization routine to estimate the best-fit parameters.
```{r, results='hide', eval=FALSE}
optPORT <- stats::nlminb(start = c(4.1, 3.9, -0.9, -8.7),
objective = RunAirGR4J,
objective = LogLikeGR4J,
lower = rep(-9.9, times = 4), upper = rep(9.9, times = 4),
control = list(trace = 1))
IniParam <- optPORT$par
iniParPORT <- optPORT$par
```
## Running 3 chains for convergence assessment
......@@ -96,15 +96,15 @@ Nota: in this example, there are relatively few iterations (2000), in order to l
With the DRAM algorithm, the covariance of the proposal is updated every 100 runs and delayed rejection is applied.
```{r, results='hide', eval=FALSE}
ListIniParam <- data.frame(Chain1 = IniParam, Chain2 = IniParam, Chain3 = IniParam,
iniParPORT <- data.frame(Chain1 = iniParPORT, Chain2 = iniParPORT, Chain3 = iniParPORT,
row.names = paste0("X", 1:4))
ListIniParam <- sweep(ListIniParam, MARGIN = 2, STATS = c(1, 0.9, 1.1), FUN = "*")
ListIniParam[ListIniParam < -9.9] <- -9.9
ListIniParam[ListIniParam > +9.9] <- +9.9
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(ListIniParam, 2, function(iListIniParam) {
FME::modMCMC(f = RunAirGR4J,
p = iListIniParam,
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,
......@@ -126,12 +126,14 @@ The result will be better with more iterations than 2000. As we kept the iterati
Note that we rescale model parameter sets of the GR4J model from the transformed space to the real space.
```{r, results='hide'}
MultDRAM <- coda::as.mcmc.list(lapply(mcmcDRAM, function(x) {
```{r, results='hide', eval=FALSE}
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
gelRub <- coda::gelman.diag(multDRAM, autoburnin = TRUE)$mpsrf
```
```{r}
gelRub
```
In addition, graphical tools can be used, with for example the [ggmcmc](https://cran.r-project.org/package=ggmcmc) package.
......@@ -141,20 +143,20 @@ In addition, graphical tools can be used, with for example the [ggmcmc](https://
First, the evolution of the Markov chains can be seen with a traceplot:
```{r, fig.width=6, fig.height=9, warning=FALSE}
ParamDRAM <- ggmcmc::ggs(MultDRAM) ## to convert objet for using by all ggs_* graphical functions
ggmcmc::ggs_traceplot(ParamDRAM)
parDRAM <- ggmcmc::ggs(multDRAM) ## to convert objet for using by all ggs_* graphical functions
ggmcmc::ggs_traceplot(parDRAM)
```
The posterior density for each parameter can then be visualised:
```{r, fig.width=6, fig.height=9, warning=FALSE}
ParamDRAM_burn <- dplyr::filter(ParamDRAM, Iteration > 1000) # to keep only the second half of the series
ggmcmc::ggs_density(ParamDRAM_burn)
burnParDRAM <- dplyr::filter(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.
```{r, fig.width=6, fig.height=6}
ggmcmc::ggs_pairs(ParamDRAM_burn, lower = list(continuous = "density"))
```{r, fig.width=6, fig.height=6, results='hide'}
ggmcmc::ggs_pairs(burnParDRAM, lower = list(continuous = "density"))
```
## Exploring further possibilities
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment