GLUE analysis - Revision history

Wbuytaert at 18:13, 12 August 2009

2009-08-12T18:13:13Z

Wbuytaert at 21:47, 23 May 2009

2009-05-23T21:47:27Z

Wbuytaert at 01:09, 12 January 2009

2009-01-12T01:09:10Z

Wbuytaert at 01:04, 12 January 2009

2009-01-12T01:04:05Z

Wbuytaert at 10:06, 2 August 2008

2008-08-02T10:06:29Z

Wbuytaert at 10:03, 2 August 2008

2008-08-02T10:03:38Z

Wbuytaert at 07:54, 7 July 2008

2008-07-07T07:54:21Z

Wbuytaert at 07:53, 7 July 2008

2008-07-07T07:53:13Z

Wbuytaert at 13:06, 26 June 2008

2008-06-26T13:06:19Z

Wbuytaert at 13:03, 26 June 2008

2008-06-26T13:03:53Z

← Older revision		Revision as of 18:13, 12 August 2009
Line 2:		Line 2:
	[[Category:Hydrology in R]]		[[Category:Hydrology in R]]

−		+	Since I am now at Imperial College London, this page has moved and can now be found on [http://wiki.r-project.org/rwiki/doku.php?id=guides:tutorials:hydrological_data_analysis the wiki page of the R-project].
−	~~===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.
−
−	~~More examples~~ and ~~scripts for doing uncertainty analysis in R~~ can be found on [http://~~www.glue-uncertainty.org/ Keith Beven's uncertainty pages].~~
−
−	~~===Libraries ===~~
−
−	~~Library Hmisc ([http://cran~~.r-project.org/~~web/packages/Hmisc~~/~~index~~.~~html CRAN]) is needed for the ''wgt.quantile()'' function. See~~ ?~~wgt.quantile for details.~~
−
−	~~library(Hmisc)~~
−	~~library(topmodel)~~
−	~~data(huagrahuma)~~
−
−	~~===Procedure===~~
−	Sample parameter sets from a prior parameter distribution. This example uses the uniform distribution, which can be sampled with "runif()". The parameters vch and psi are not used in this example but need to be initialised. For an explanation of the parameter values see the page on [[Running Topmodel]]. We take initially 1000 parameter sets. Depending on the performance of the model, many of them will be rejected as non-behavioural. If you end up with too few parameter behavioural parameter sets, you will need to increase your sample size.
−
−	~~qs0 <- runif(1000, min = 0.0001, max = 0.00025)~~
−	~~lnTe <- runif(1000, min = -2, max = 3)~~
−	~~m <- runif(1000, min = 0, max = 0.1)~~
−	~~Sr0 <- runif(1000, min = 0, max = 0.2)~~
−	~~Srmax <- runif(1000, min = 0, max = 0.1)~~
−	~~td <- runif(1000, min = 0, max = 3)~~
−	~~vch <- 1000~~
−	~~vr <- runif(1000, min = 100, max = 2500)~~
−	~~k0 <- runif(1000, min = 0, max~~ = ~~10)~~
−	~~CD <- runif(1000, min = 0, max = 5)~~
−	~~dt <- 0.25~~
−
−	~~parameters <- cbind(qs0,lnTe,m,Sr0,Srmax,td,vch,vr,k0,CD,dt)~~
−
−	~~Now each row in the matrix ''parameters'' contains a sampled parameter set. Run the model for the calibration period with the each set to obtain the simulated discharge:~~
−
−	~~Qsim <- topmodel(parameters[1,],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 <- NSeff(Qobs,Qsim)~~
−
−	~~Decide whether the parameter set is behavioural or not and retain the parameter set if behavioural~~
−
−	~~NOTE~~: ~~in GLUE, this decision is 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)~~
−
−	~~Calculate the quantiles for each timestep. Here we take the 0.05 and 0.95 quantiles, resulting in 90% prediction limits.~~
−
−	~~limits <- apply(Qsim,1,"wtd.quantile", weights = Eff, probs = c(0.05,0.95), normwt=T)~~
−
−	~~===Final notes===~~
−
−	* If topmodel is used, some loops can be avoided. topmodel() can work on entire parameter set matrices, which makes things faster. 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.

@@ Line 32: / Line 32: @@
     dt    <- 0.25
-    parameters <- cbind(qs0,lnTe,m,Sr0,Srmax,td,vch,vr,k0,psi,dtheta,dt)
+    parameters <- cbind(qs0,lnTe,m,Sr0,Srmax,td,vch,vr,k0,CD,dt)
 Now each row in the matrix ''parameters'' contains a sampled parameter set. Run the model for the calibration period with the each set to obtain the simulated discharge:

@@ Line 18: / Line 18: @@
 ===Procedure===
-Sample a parameter set from a prior parameter distribution. This example uses the uniform distribution, which can be sampled with "runif()". The parameters vch and psi are not used in this example but need to be initialised. For an explanation of the parameter values see the page on [[Running Topmodel]].
+Sample parameter sets from a prior parameter distribution. This example uses the uniform distribution, which can be sampled with "runif()". The parameters vch and psi are not used in this example but need to be initialised. For an explanation of the parameter values see the page on [[Running Topmodel]]. We take initially 1000 parameter sets. Depending on the performance of the model, many of them will be rejected as non-behavioural. If you end up with too few parameter behavioural parameter sets, you will need to increase your sample size.
-    qs0   <- runif(1, min = 0.0001, max = 0.00025)
+    qs0   <- runif(1000, min = 0.0001, max = 0.00025)
-    lnTe  <- runif(1, min = -2, max = 3)
+    lnTe  <- runif(1000, min = -2, max = 3)
-    m     <- runif(1, min = 0, max = 0.1)
+    m     <- runif(1000, min = 0, max = 0.1)
-    Sr0   <- runif(1, min = 0, max = 0.2)
+    Sr0   <- runif(1000, min = 0, max = 0.2)
-    Srmax <- runif(1, min = 0, max = 0.1)
+    Srmax <- runif(1000, min = 0, max = 0.1)
-    td    <- runif(1, min = 0, max = 3)
+    td    <- runif(1000, min = 0, max = 3)
     vch   <- 1000
-    vr    <- runif(1, min = 100, max = 2500)
+    vr    <- runif(1000, min = 100, max = 2500)
-    k0    <- runif(1, min = 0, max = 10)
+    k0    <- runif(1000, min = 0, max = 10)
-    CD    <- runif(1, min = 0, max = 5)
+    CD    <- runif(1000, min = 0, max = 5)
     dt    <- 0.25
-    parameters <- c(qs0,lnTe,m,Sr0,Srmax,td,vch,vr,k0,psi,dtheta,dt)
+    parameters <- cbind(qs0,lnTe,m,Sr0,Srmax,td,vch,vr,k0,psi,dtheta,dt)
-Run the model for the calibration period with the generated parameter set to obtain the simulated discharge
+Now each row in the matrix ''parameters'' contains a sampled parameter set. Run the model for the calibration period with the each set to obtain the simulated discharge:
-    Qsim <- topmodel(parameters,topidx,delay,rain,ET0)
+    Qsim <- topmodel(parameters[1,],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:
@@ Line 44: / Line 44: @@
 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''
+NOTE: in GLUE, this decision is 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) {

@@ Line 15: / Line 15: @@
     library(Hmisc)
     library(topmodel)
-    data(Huagrahuma)
+    data(huagrahuma)
 ===Procedure===
@@ Line 29: / Line 29: @@
     vr    <- runif(1, min = 100, max = 2500)
     k0    <- runif(1, min = 0, max = 10)
-    psi   <- 1
+    CD    <- runif(1, min = 0, max = 5)
     dt    <- 0.25
-    parameters<-c(qs0,lnTe,m,Sr0,Srmax,td,vch,vr,k0,psi,dtheta,dt)
+    parameters <- c(qs0,lnTe,m,Sr0,Srmax,td,vch,vr,k0,psi,dtheta,dt)
 Run the model for the calibration period with the generated parameter set to obtain the simulated discharge
@@ Line 41: / Line 40: @@
 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)
+    eff <- NSeff(Qobs,Qsim)
 Decide whether the parameter set is behavioural or not and retain the parameter set if behavioural
@@ Line 72: / Line 69: @@
     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.
+Calculate the quantiles for each timestep. Here we take the 0.05 and 0.95 quantiles, resulting in 90% prediction limits.
-The final prediciton limits are stored in Llimit and Ulimit.
+   limits <- apply(Qsim,1,"wtd.quantile", weights = Eff, probs = c(0.05,0.95), normwt=T)
 ===Final notes===
-* If topmodel is used, some loops can be avoided because topmodel() can work on entire parameter set matrices, which makes things faster. It can also return the Nash-Sutcliffe efficiency directly (see the topmodel page)
+* If topmodel is used, some loops can be avoided. topmodel() can work on entire parameter set matrices, which makes things faster. 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

@@ Line 32: / Line 32: @@
     dtheta<- runif(1, min = 0, max = 5)
     dt    <- 0.25
-−
     parameters<-c(qs0,lnTe,m,Sr0,Srmax,td,vch,vr,k0,psi,dtheta,dt)
@@ Line 61: / Line 61: @@
     predicted.qsim <- model(behavioural.parameters[,1], rain, ...)
-−
     for(i in 2:dim(param.set)[2]) {
        qsim <- model(behavioural.parameters[,i], rain, ...)

@@ Line 7: / Line 7: @@
 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.
-More inspiration for doing uncertainty analysis in R can be found on [Keith Beven's uncertainty pages http://www.glue-uncertainty.org/]
+More examples and scripts for doing uncertainty analysis in R can be found on [http://www.glue-uncertainty.org/ Keith Beven's uncertainty pages].
 ===Libraries ===

← Older revision		Revision as of 07:53, 7 July 2008
Line 6:		Line 6:

	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.		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.
		+
		+	More inspiration for doing uncertainty analysis in R can be found on [Keith Beven's uncertainty pages http://www.glue-uncertainty.org/]

	===Libraries ===		===Libraries ===

@@ Line 82: / Line 82: @@
     Llimit <- 0
-Now we calculate the quantiles for each timestep (this can also be done with apply() function)
+Now we calculate the quantiles for each timestep (this can also be done with the function apply())
     for(i in 1:dim(predicted.qsim)[1]) {
@@ Line 89: / Line 89: @@
+    }
-If everything went well, the final prediciton limits are in Llimit and Ulimit.
+The final prediciton limits are stored 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)
+* If topmodel is used, some loops can be avoided because topmodel() can work on entire parameter set matrices, which makes things faster. 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