From 91d6c8ec47c28812ed56d3b052044fa9deb5a2db Mon Sep 17 00:00:00 2001
From: Delaigue Olivier <olivier.delaigue@irstea.priv>
Date: Tue, 23 Oct 2018 15:06:53 +0200
Subject: [PATCH] v1.1.2.18 CLEAN: syntaxe revision of vignette param_mcmc
---
DESCRIPTION | 2 +-
NEWS.rmd | 2 +-
vignettes/V02.2_param_mcmc.Rmd | 64 ++++++++++++++++++----------------
3 files changed, 35 insertions(+), 33 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index df670acc..9d86cda5 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
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")),
diff --git a/NEWS.rmd b/NEWS.rmd
index 57069879..94922af1 100644
--- a/NEWS.rmd
+++ b/NEWS.rmd
@@ -13,7 +13,7 @@ output:
-### 1.1.2.17 Release Notes (2018-10-23)
+### 1.1.2.18 Release Notes (2018-10-23)
diff --git a/vignettes/V02.2_param_mcmc.Rmd b/vignettes/V02.2_param_mcmc.Rmd
index f618b78b..4555cf0c 100644
--- a/vignettes/V02.2_param_mcmc.Rmd
+++ b/vignettes/V02.2_param_mcmc.Rmd
@@ -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
--
GitLab