fitQmapDIST: Quantile mapping using distribution derived transformations

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

fitQmapDIST(obs, mod, ...)
# S3 method for default
fitQmapDIST(obs,mod,distr="berngamma",start.fun,
qstep=NULL,mlepar,...)
# S3 method for matrix
fitQmapDIST(obs, mod, ...)
# S3 method for data.frame
fitQmapDIST(obs, mod, ...)

doQmapDIST(x,fobj,...)
# S3 method for default
doQmapDIST(x,fobj,...)
# S3 method for matrix
doQmapDIST(x,fobj,...)
# S3 method for data.frame
doQmapDIST(x,fobj,...)

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.

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

Author

Lukas Gudmundsson

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_{o}^{- 1} (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.

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>.

Examples

Run this code

data(obsprecip)
data(modprecip)

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


# adjusting settings of the maximum likelyhood estimator 
# here changing the convergence criterion in the optimizer 'optim'
# Note: the selected value might be to large for true applications.
# Please check for your context.
qm.lnorm.fit <- fitQmapDIST(obsprecip[,1],modprecip[,1],
                      distr="bernlnorm",
                      qstep=0.001,
                      mlepar = list(control = list(
                            reltol = 1e-2
                      )))
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)
}

Run the code above in your browser using DataLab

The Learning Leader's Guide to AI Literacy