Quantile mapping using distribution derived transformations

Share:

Description

fitQmapDIST fits a theoretical distribution to observed and to modelled time series and returns these parameters as well as a transfer function derived from the distribution. doQmapDIST uses the transfer function to transform the distribution of the modelled data to match the distribution of the observations.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
fitQmapDIST(obs, mod, ...)
## Default S3 method:
fitQmapDIST(obs,mod,distr="berngamma",start.fun,
qstep=NULL,mlepar,...)
## S3 method for class 'matrix'
fitQmapDIST(obs, mod, ...)
## S3 method for class 'data.frame'
fitQmapDIST(obs, mod, ...)

doQmapDIST(x,fobj,...)
## Default S3 method:
doQmapDIST(x,fobj,...)
## S3 method for class 'matrix'
doQmapDIST(x,fobj,...)
## S3 method for class 'data.frame'
doQmapDIST(x,fobj,...)

Arguments

obs

numeric vector, column matrix or data.frame with observed time series.

mod

numeric vector, column matrix or data.frame with modelled time series, corresponding to obs.

distr

A character string "name" naming a distribution for which the corresponding density function (dname), the corresponding distribution function (pname) and the quantile function (qname) must be defined (see for example GammaDist, berngamma or bernweibull.

start.fun

function estimating starting values for parameter optimisation. Default starting values are provided for berngamma, bernweibull, bernlnorm, bernexp and the distributions mentioned in the documentation of mledist.

qstep

NULL or a numeric value between 0 and 1. If !is.null(qstep) than mod and obs are aggregated to quantiles before model identification as:

quantile(x,probs=seq(0,1,by=qstep). This effectively reduces the sample-size and can be used to speedup computations - but may render estimates less reliable.

mlepar

a named list. Names correspond to parameters passed to mledist note that start may be overwritten by start.fun See examples.

x

numeric vector or a column matrix of modelled time series

fobj

output from fitQmapDIST

...

Further arguments passed to methods

Details

Quantile mapping using distribution derived transformations to adjust the distribution of a modelled variable (P_m) such that it matches the distribution of an observed variable (P_o). The distribution derived transfer function is defined as

P_o=F^{-1}_o(F_m(P_m))

where F is a CDF and F^{-1} is the corresponding quantile function (inverse CDF). The subscripts o and m indicate parameters of the distribution that correspond to observed and modelled data respectively.

Value

fitQmapDIST returns an object of class fitQmapDIST containing following elements:

tfun

The function used to transform the distribution of modelled values such that the distribution of observations. The function is build internally based on the distribution function ("pname") and quantile function ("qname") corresponding to distr.

par

A matrix. The (named) columns correspond to the parameters of the distribution specified in distr estimated for the observed (suffix .o) and the modelled (suffix .m) data. The rows correspond to each pair of time series in obs and mod.

doQmapDIST returns a numeric vector, matrix or data.frame depending on the format of x.

Author(s)

Lukas Gudmundsson

References

Piani, C.; Haerter, J. & Coppola, E. Statistical bias correction for daily precipitation in regional climate models over Europe. Theoretical and Applied Climatology, 2010, 99, 187-192, doi:10.1007/s00704-009-0134-9.

Li, H.; Sheffield, J. & Wood, E. F. Bias correction of monthly precipitation and temperature fields from Intergovernmental Panel on Climate Change AR4 models using equidistant quantile matching. J. Geophys. Res., 2010, 115, D10101, doi:10.1029/2009JD012882.

Ines, A. V. & Hansen, J. W. Bias correction of daily GCM rainfall for crop simulation studies. Agricultural and Forest Meteorology, 2006, 138, 44 - 53, doi: 10.1016/j.agrformet.2006.03.009.

For a general assessment of the methods see:

Gudmundsson, L.; Bremnes, J. B.; Haugen, J. E. & Engen-Skaugen, T. Technical Note: Downscaling RCM precipitation to the station scale using statistical transformations - a comparison of methods. Hydrology and Earth System Sciences, 2012, 16, 3383-3390, doi:10.5194/hess-16-3383-2012.

See Also

doQmap, startberngamma, berngamma, startbernweibull, bernweibull, startbernlnorm, bernlnorm, startbernexp, bernexp, mledist, fitdist

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
data(obsprecip)
data(modprecip)

qm.fit <- fitQmapDIST(obsprecip[,1],modprecip[,1],
                      distr="berngamma",
                      qstep=0.001)
qm <- doQmapDIST(modprecip[,1],qm.fit)


qm.lnorm.fit <- fitQmapDIST(obsprecip[,1],modprecip[,1],
                      distr="bernlnorm",
                      qstep=0.001)
qm.lnorm <- doQmapDIST(modprecip[,1],qm.lnorm.fit)


qm.weibu.fit <- fitQmapDIST(obsprecip[,1],modprecip[,1],
                      distr="bernweibull",
                      qstep=0.001)
qm.weibu <- doQmapDIST(modprecip[,1],qm.weibu.fit)

qm.exp.fit <- fitQmapDIST(sqrt(obsprecip[,1]),sqrt(modprecip[,1]),
                      distr="bernexp",
                      qstep=0.001)
qm.exp <- doQmapDIST(sqrt(modprecip[,1]),qm.exp.fit)^2


## utility function. 
## plots are easier to investigate if
## precipitation data are sqrt transformed
sqrtquant <- function(x,qstep=0.01){
  qq <- quantile(x,prob=seq(0,1,by=qstep))
  sqrt(qq)
}

plot(sqrtquant(modprecip[,1]),
     sqrtquant(obsprecip[,1]))
lines(sqrtquant(modprecip[,1]),
      sqrtquant(qm),col="red")
lines(sqrtquant(modprecip[,1]),
      sqrtquant(qm.lnorm),col="blue")
lines(sqrtquant(modprecip[,1]),
      sqrtquant(qm.weibu),col="green")
lines(sqrtquant(modprecip[,1]),
      sqrtquant(qm.exp),col="orange")
legend("topleft",
       legend=c("berngamma","bernlnorm","bernweibull","bernexp"),
       lty=1,
       col=c("red","blue","green","orange"))

## effect of qstep on speed of fitting process:
system.time(
qm.a.fit <- fitQmapDIST(obsprecip[,2],modprecip[,2],
                       distr="berngamma",
                       start.fun=startberngamma,
                       qstep=0.001)
)

system.time(
qm.b.fit <- fitQmapDIST(obsprecip[,2],modprecip[,2],
                       distr="berngamma",
                       start.fun=startberngamma,
                       qstep=0.01)
)

qm.a <- doQmapDIST(modprecip[,2],qm.a.fit)
qm.b <- doQmapDIST(modprecip[,2],qm.b.fit)

plot(sqrtquant(modprecip[,2]),
     sqrtquant(obsprecip[,2]))
lines(sqrtquant(modprecip[,2]),
     sqrtquant(qm.a),col="red")
lines(sqrtquant(modprecip[,2]),
     sqrtquant(qm.b),col="blue")
legend("topleft",
       legend=c("qstep=0.001","qstep=0.01"),
       col=c("red","blue"),
       lty=1)


## method for matrix
## the sqrt() transformation renders the
## fitting procedure more stable
qm2.fit <- fitQmapDIST(sqrt(obsprecip),sqrt(modprecip),
                       distr="berngamma",
                       qstep=0.001)
qm2 <- doQmapDIST(sqrt(modprecip),qm2.fit)^2

if(!any(is.na(qm2.fit$par))){
  op <- par(mfrow=c(1,3))
  for(i in 1:3){
    plot(sqrtquant(modprecip[,i]),
         sqrtquant(obsprecip[,i]))
    lines(sqrtquant(modprecip[,i]),
          sqrtquant(qm2[,i]),col="red")
  }
 par(op)
}