# R/fitQmapRQUANT.R In qmap: Statistical Transformations for Post-Processing Climate Model Output

```fitQmapRQUANT <- function(obs,mod,...)
UseMethod("fitQmapRQUANT")

fitQmapRQUANT.default <- function(obs,mod,wet.day=TRUE,qstep=0.01,
nlls = 10,nboot = 10,...){
### fits an nonparametric transfere function for quantile
### mapping using local linear least square regression
###
### Based on code from John Bjornar Bremnes
ys <- na.omit(obs)
xs <- na.omit(mod)
if(length(xs)!=length(ys)){
hn <- min(length(xs),length(ys))
## quantile algorithm 'type=8' appeares to
## be reccomended. See help(quantile) and
## Hyndman & Fan (1996) mentioned therein
ys <- quantile(ys,seq(0,1,length.out=hn),type=8)
xs <- quantile(xs,seq(0,1,length.out=hn),type=8)
} else {
xs <- sort(xs)
ys <- sort(ys)
}
if(is.numeric(wet.day)){
q0 <- ys>=wet.day
ys <- ys[q0]
xs <- xs[q0]
} else if(is.logical(wet.day)){
if(wet.day){
q0 <- ys>0
ys <- ys[q0]
xs <- xs[q0]
wet.day <- xs[1]
names(wet.day) <- NULL
} else {
wet.day <- NULL
}
} else {
stop("'wet.day' should be 'numeric' or 'logical'")
}
nn <- length(ys)
## define predictor values at which a fit is wanted
newx <- quantile(xs, probs=seq(0,1,by=qstep),type=8)
fit   <- array(NA, dim=c(length(newx),2,nboot))
nlls2 <- min(nlls, nn)
for (j in 1:nboot) {
if(nboot==1){
xss <- sort.default(xs)
yss <- sort.default(ys)
} else {
xss <- sort.default(sample(xs, size=nn, replace=FALSE))
yss <- sort.default(sample(ys, size=nn, replace=FALSE))
}
for (i in 1:length(newx)) {
xc     <- xss - newx[i]
mdist  <- sort.default(abs(xc))[nlls2]
k      <- abs(xc) <= mdist
xc     <- cbind(1, xc[k])
a      <- t.default(xc) %*% xc
if ( abs(d <- a[1,1]*a[2,2] - a[1,2]*a[2,1]) < 1e-10 )  # singular? <=> det(x^tx)=0?
fit[i,1,j] <- mean(yss[k])
else {
b         <- t.default(xc) %*% yss[k]
fit[i,,j] <- c(b[1]*a[2,2]-b[2]*a[1,2], b[2]*a[1,1]-b[1]*a[2,1]) / d
}
}
}
if (!is.null(wet.day))
if (any(k <- fit[,1,] < 0))
fit[,1,][which(k)] <- 0
fitted <- if(nboot>1){
rowMeans(fit[,1,],na.rm=TRUE)
} else fit[,1,]
slope <- if(nboot>1){
rowMeans(fit[,2,][,,drop=FALSE],na.rm=TRUE)
} else fit[,2,]
slope <- slope[range(which(!is.na(slope)))]
slope <- matrix(slope,ncol=1,
dimnames=list(c("lower","upper"),NULL))
ppar <- list(modq=matrix(newx,ncol=1),
fitq=matrix(fitted,ncol=1),
slope=slope)
op <- list(par=ppar,
wet.day=wet.day)
class(op) <- "fitQmapRQUANT"
return(op)
}

fitQmapRQUANT.matrix <- function(obs,mod,...){
if(ncol(mod)!=ncol(obs))
stop("'mod' and 'obs' need the same number of columns")
NN <- ncol(mod)
hind <- 1:NN
names(hind) <- colnames(mod)
xx <- lapply(hind,function(i){
tr <- try(fitQmapRQUANT.default(obs=obs[,i],mod=mod[,i],...),
silent=TRUE)
if(any(class(tr)=="try-error")){
warning("model identification for ",names(hind)[i],
" failed. NA's produced")
NULL
} else{
tr
}
})
xx.NULL <- sapply(xx,is.null)
## tfun <- xx[!xx.NULL][[1]]\$tfun
modq <- lapply(xx,function(x)x\$par\$modq)
fitq <- lapply(xx,function(x)x\$par\$fitq)
slope <- lapply(xx,function(x)x\$par\$slope)
nq <- nrow(xx[!xx.NULL][[1]]\$par\$modq)
nq <- matrix(NA,nrow=nq)
modq[xx.NULL] <- list(nq)
fitq[xx.NULL] <- list(nq)
slope[xx.NULL] <- list(matrix(NA,ncol=1,nrow=2,
dimnames=list(c("lower","upper"),NULL)))
modq <- do.call(cbind,modq)
fitq <- do.call(cbind,fitq)
slope <- do.call(cbind,slope)
colnames(modq) <- names(xx)
colnames(fitq) <- names(xx)
colnames(slope) <- names(xx)
wday <- lapply(xx,function(x)x\$wet.day)
wday[xx.NULL] <- NA
wday <- do.call(c,wday)
xx <- list(par=list(modq=modq,fitq=fitq,slope=slope),
wet.day=wday)
class(xx) <- c("fitQmapRQUANT")
return(xx)
}

fitQmapRQUANT.data.frame <- function(obs,mod,...){
obs <- as.matrix(obs)
mod <- as.matrix(mod)
fitQmapRQUANT.matrix(obs,mod,...)
}
```

## Try the qmap package in your browser

Any scripts or data that you put into this service are public.

qmap documentation built on May 1, 2019, 7:31 p.m.