diff --git a/vignettes/V02.1_param_optim.Rmd b/vignettes/V02.1_param_optim.Rmd
index 2dd121e92d89e6838e3c7db3af73322e74080400..b71d6b8fa6006974ef663e1f38a60a1bab43e43e 100644
--- a/vignettes/V02.1_param_optim.Rmd
+++ b/vignettes/V02.1_param_optim.Rmd
@@ -45,6 +45,13 @@ Please note that the calibration period is defined in the `CreateRunOptions()` f
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).
+
+```{r, results='hide', eval=FALSE}
+RunOptions <- airGR::CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
+ IndPeriod_Run = Ind_Run,
+ Outputs_Sim = "Qsim")
+```
Regarding the different optimization strategies presented here, we refer to each package for in-depth information about the description of the methods used.
@@ -58,7 +65,6 @@ Here we choose to minimize the root mean square error.
The change of the repository from the "real" parameter space to a "transformed" space ensures homogeneity of displacement in the different dimensions of the parameter space during the step-by-step procedure of the calibration algorithm of the model.
-
```{r, warning=FALSE, results='hide', eval=FALSE}
OptimGR4J <- function(ParamOptim) {
## Transformation of the parameter set to real space
@@ -78,6 +84,7 @@ OptimGR4J <- function(ParamOptim) {
In addition, we need to define the lower and upper bounds of the four **GR4J** parameters in the transformed parameter space:
+
```{r, warning=FALSE, results='hide', eval=FALSE}
lowerGR4J <- rep(-9.99, times = 4)
upperGR4J <- rep(+9.99, times = 4)
@@ -86,6 +93,7 @@ upperGR4J <- rep(+9.99, times = 4)
# Local optimization
We start with a local optimization strategy by using the PORT routines (using the `nlminb()` of the `stats` package) and by setting a starting point in the transformed parameter space:
+
```{r, warning=FALSE, results='hide', eval=FALSE}
startGR4J <- c(4.1, 3.9, -0.9, -8.7)
optPORT <- stats::nlminb(start = startGR4J,
@@ -93,13 +101,17 @@ optPORT <- stats::nlminb(start = startGR4J,
lower = lowerGR4J, upper = upperGR4J,
control = list(trace = 1))
```
+
The RMSE value reaches a local minimum value after 35 iterations.
We can also try a multi-start approach to test the consistency of the local optimization.
Here we use the same grid used for the filtering step of the Michel's calibration strategy (`Calibration_Michel()` function).
For each starting point, a local optimization is performed.
+
```{r, warning=FALSE, results='hide', eval=FALSE}
-startGR4J <- expand.grid(data.frame(CalibOptions$StartParamDistrib))
+StarttGR4JDistrib <- TransfoParam_GR4J(ParamIn = CalibOptions$StartParamDistrib,
+ Direction = "RT")
+startGR4J <- expand.grid(data.frame(StarttGR4JDistrib))
optPORT_ <- function(x) {
opt <- stats::nlminb(start = x,
objective = OptimGR4J,
@@ -110,6 +122,7 @@ listOptPORT <- apply(startGR4J, MARGIN = 1, FUN = optPORT_)
```
We can then extract the best parameter sets and the value of the performance criteria:
+
```{r, warning=FALSE, results='hide', eval=FALSE}
parPORT <- t(sapply(listOptPORT, function(x) x$par))
objPORT <- sapply(listOptPORT, function(x) x$objective)
@@ -117,6 +130,7 @@ resPORT <- data.frame(parPORT, RMSE = objPORT)
```
As can be seen below, the optimum performance criterion values (column *objective*) can differ from the global optimum value in many cases, resulting in various parameter sets.
+
```{r, warning=FALSE}
summary(resPORT)
```
@@ -134,6 +148,7 @@ Here we use the following R implementation of some popular strategies:
* [Rmalschains: memetic algorithms](https://cran.r-project.org/package=Rmalschains)
## Differential Evolution
+
```{r, warning=FALSE, results='hide', eval=FALSE}
optDE <- DEoptim::DEoptim(fn = OptimGR4J,
lower = lowerGR4J, upper = upperGR4J,
@@ -142,6 +157,7 @@ optDE <- DEoptim::DEoptim(fn = OptimGR4J,
## Particle Swarm
+
```{r, warning=FALSE, results='hide', message=FALSE, eval=FALSE}
optPSO <- hydroPSO::hydroPSO(fn = OptimGR4J,
lower = lowerGR4J, upper = upperGR4J,
@@ -149,6 +165,7 @@ optPSO <- hydroPSO::hydroPSO(fn = OptimGR4J,
```
## MA-LS-Chains
+
```{r, warning=FALSE, results='hide', eval=FALSE}
optMALS <- Rmalschains::malschains(fn = OptimGR4J,
lower = lowerGR4J, upper = upperGR4J,
@@ -214,6 +231,7 @@ MOptimGR4J <- function(i) {
```
## caRamel
+
caRamel is a multiobjective evolutionary algorithm combining the MEAS algorithm and the NGSA-II algorithm.
```{r, warning=FALSE, results='hide', eval=FALSE}
diff --git a/vignettes/V02.2_param_mcmc.Rmd b/vignettes/V02.2_param_mcmc.Rmd
index 02e9c5520c37833183840c03c4cb226d062b8888..7d6091ceafe3c3433ef3a65db84ffb8a9b428648 100644
--- a/vignettes/V02.2_param_mcmc.Rmd
+++ b/vignettes/V02.2_param_mcmc.Rmd
@@ -38,6 +38,14 @@ We use the **GR4J** model and we assume that the R global environment contains d
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).
+
+```{r, results='hide', eval=FALSE}
+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](V02.1_param_optim.html) vignette for explanations on how to plug in a parameter estimation algorithm to **airGR**.
@@ -49,6 +57,7 @@ We show how to use the DRAM algorithm for SLS Bayesian inference, with the `modM
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.
+
```{r, echo=TRUE, eval=FALSE, purl=FALSE}
Likelihood <- sum((ObsY - ModY)^2, na.rm = TRUE)^(-sum(!is.na(ObsY)) / 2)
LogLike <- -2 * log(Likelihood)
@@ -96,8 +105,10 @@ 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}
-iniParPORT <- data.frame(Chain1 = iniParPORT, Chain2 = iniParPORT, Chain3 = iniParPORT,
- row.names = paste0("X", 1:4))
+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
@@ -111,8 +122,8 @@ mcmcDRAM <- apply(iniParPORT, MARGIN = 2, FUN = function(iIniParPORT) {
jump = 0.01,
outputlength = 2000,
burninlength = 0,
- updatecov = 100, ## Adaptative Metropolis
- ntrydr = 2) ## Delayed Rejection
+ updatecov = 100, ## adaptative Metropolis (AM)
+ ntrydr = 2) ## delayed rejection (RD)
})
```
@@ -143,7 +154,7 @@ 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}
-parDRAM <- ggmcmc::ggs(multDRAM) ## to convert objet for using by all ggs_* graphical functions
+parDRAM <- ggmcmc::ggs(multDRAM) ## to convert object for using by all ggs_* graphical functions
ggmcmc::ggs_traceplot(parDRAM)
```