Difference between revisions of "GLUE analysis"

From SourceWiki
Jump to navigation Jump to search
Line 16: Line 16:
  
 
===Procedure===
 
===Procedure===
Sample a parameter set from a prior parameter distribution. This example uses the uniform distribution, which can be sampled between 0 and 1 with ''runif()''. We then scale to the appropriate range and offset.
+
Sample a parameter set from a prior parameter distribution. This example uses the uniform distribution, which can be sampled between 0 and 1 with ''runif()''. We then scale to the appropriate range and offset. The parameters vch and psi are not used in this example but need to be initialised.
  
 
   qs0  <- runif(1)*0.00012+0.00012
 
   qs0  <- runif(1)*0.00012+0.00012
Line 47: Line 47:
 
Decide whether the parameter set is behavioural or not and retain the parameter set if behavioural
 
Decide whether the parameter set is behavioural or not and retain the parameter set if behavioural
  
NOTE: this decision is again subjective. For a more scientifically sound determination of the behavioural limit, see Beven (2006). Here we will use an efficiency of 0.6 as a threshold.
+
NOTE: this decision is again subjective. For a more scientifically sound determination of the behavioural limit, see Beven (2006). Here we will use an efficiency of 0.6 as a threshold. The efficiency, parameter set and simulated discharge of a behavioural run are stored in resp. the objects ''total.eff'', ''total.param.set'' and ''total.qsim''
The efficiency, parameter set and simulated discharge of a behavioural run are stored in resp. the objects total.eff, total.param.set and total.qsim
 
  
 
   if(eff > 0.6) {
 
   if(eff > 0.6) {
Line 57: Line 56:
 
The above procedure should be repeated until enough behavioural runs are obtained (e.g. using a while()-loop)
 
The above procedure should be repeated until enough behavioural runs are obtained (e.g. using a while()-loop)
  
In the next section it is assumed that each column of behavioural.parameters contains a behavioural parameter set and that the corresponding efficiency is found at the same location in the total.eff vector
+
If the above code is used, each column of the matrix ''behavioural.parameters'' contains a behavioural parameter set. The corresponding performance is found at the same location in the vector ''total.eff''
  
 
Rerun the model for the prediction period, using each of the behavioural parameter sets. The simulated discharges are stored in the columns of a matrix called predicted.qsim
 
Rerun the model for the prediction period, using each of the behavioural parameter sets. The simulated discharges are stored in the columns of a matrix called predicted.qsim

Revision as of 13:03, 26 June 2008


Introduction

This is an example of R code to perform a Generalised Likelihood Uncertainty Estimation (GLUE) on a hydrological model. In the example we use topmodel (implemented as topmodel()) on a catchment in the Ecuadorian Andes (Huagrahuma). Other models can be used as long as they are implemented as an R function.

Libraries needed

Library Hmisc (CRAN) is needed for the wgt.quantile() function. See ?wgt.quantile for details.

  > library(Hmisc)
  > library(topmodel)
  > data(Huagrahuma)

Procedure

Sample a parameter set from a prior parameter distribution. This example uses the uniform distribution, which can be sampled between 0 and 1 with runif(). We then scale to the appropriate range and offset. The parameters vch and psi are not used in this example but need to be initialised.

  qs0   <- runif(1)*0.00012+0.00012
  lnTe  <- runif(1)*5-2
  m     <- runif(1)*0.1
  Sr0   <- runif(1)*0.2
  Srmax <- runif(1)*0.1
  td    <- runif(1)*3
  vch   <- 1000
  vr    <- 100+runif(1)*2500
  k0    <- runif(1)*10
  psi   <- 1
  dtheta<- runif(1)*5
  dt    <- 0.25
  parameters<-c(qs0,lnTe,m,Sr0,Srmax,td,vch,vr,k0,psi,dtheta,dt)


(TODO: explain the parameters and the chosen ranges) Run the model for the calibration period with the generated parameter set to obtain the simulated discharge

  Qsim <- topmodel(parameters,topidx,delay,rain,ET0)

Calculate the likelihood of this parameter set using the simulated and observed discharge. The choice of the likelihood function is up to the user, but the Nash - Sutcliffe efficiency is given here as an example:

  eff <- 1 - sum((Qobs - Qsim)^2) / sum((Qobs-mean(Qobs))^2)

NOTE: calculation of the Nash-Sutcliffe efficiency is also implimented in NSeff() of the topmodel-package.

Decide whether the parameter set is behavioural or not and retain the parameter set if behavioural

NOTE: this decision is again subjective. For a more scientifically sound determination of the behavioural limit, see Beven (2006). Here we will use an efficiency of 0.6 as a threshold. The efficiency, parameter set and simulated discharge of a behavioural run are stored in resp. the objects total.eff, total.param.set and total.qsim

  if(eff > 0.6) {
     total.eff <- c(total.eff,eff)
     behavioural.parameters <- cbind(behavioural.parameters,parameters)
  }

The above procedure should be repeated until enough behavioural runs are obtained (e.g. using a while()-loop)

If the above code is used, each column of the matrix behavioural.parameters contains a behavioural parameter set. The corresponding performance is found at the same location in the vector total.eff

Rerun the model for the prediction period, using each of the behavioural parameter sets. The simulated discharges are stored in the columns of a matrix called predicted.qsim

  predicted.qsim <- model(behavioural.parameters[,1], rain, ...)
  for(i in 2:dim(param.set)[2]) {
     qsim <- model(behavioural.parameters[,i], rain, ...)
     predicted.qsim <- cbind(predicted.obs,qobs)
  }

Normalise the efficiencies so that they sum up to 1:

  eff <- eff - 0.6
  eff <- eff/sum(eff)

Define a quantile for the prediction bounds. Here we take the 0.05 and 0.95 quantiles resulting in 90% prediction limits.

  lower <- 0.05
  upper <- 0.95

Create the objects in which we will store the prediction limits:

  Ulimit <- 0
  Llimit <- 0

Now we calculate the quantiles for each timestep (this can also be done with apply() function)

  for(i in 1:dim(predicted.qsim)[1]) {
     Llimit[i] <- wtd.quantile(predicted.qsim[i,],weights = eff, probs = lower, normwt=T)
     Ulimit[i] <- wtd.quantile(predicted.qsim[i,],weights = eff, probs = upper, normwt=T)
  }

If everything went well, the final prediciton limits are in Llimit and Ulimit.

Final notes

  • If topmodel is used, some loops can be avoided because topmodel() can work on entire parameter set matrices. It can also return the Nash-Sutcliffe efficiency directly (see the topmodel page)
  • The procedure can be very memory intensive because all simulated discharges for all parameter sets are stored in memory (the matrix predicted.qsim). If the model can give output per timestep, the above procedure can be repeated for each timestep separately to reduce memory usage

References

  • Beven, K., and Binley, A. The future of distributed models: Model calibration and uncertainty prediction. Hydrological Processes 6 (1992), 279-298.
  • Beven, K. A manifesto for the equifinality thesis. Journal of Hydrology 320 (2006), 18-36.