R/qdks.R

In waternumbers/FloodForT: Flood Forecasting Toolbox

#' A quick and dirty non-parametric quantile estimator based on pooling and linear interpolation
#' 
#' This function runs a estimator of the quantiles of y depending upon x.
#' It pools the values of y into bins based on the quantiles of x then interpolates in 
#' a linear fashion assuming the quantiles occur at the bin centre. Out of range extrapolation
#' depends upon rule - see \code{approx} function in the stats package
#'
#' @param y Values whose quantiles are sought
#' @param x Explanatory variable
#' @param q Quantiles estimated. Default is 0.05,0.1,...,0.95
#' @param xout Value of x to return quantiles at. NA (default) returns values for seq(min(x),max(x),length=100)
#' @param rule See \code{approx} function in stats. default is 2, take closest end point
#' 
#' @keywords FloorForT, ModelBuild
#' @export
#' @examples
#' # Not Run
#' # ModelBuild()
qdks <- function(y,x,q = seq(5,95,5)/100, xout=NA,rule = 2){
  
  # remove missing values
  idx <- is.finite(y) & is.finite(x)
  y <- y[idx]
  x <- x[idx]
  
  if(length(x)==0) return(matrix(NA,length(xout),length(q)))
  
  if(all(is.na(xout))) xout <- seq(min(x,na.rm=TRUE),max(x,na.rm=TRUE),length=100)
  # number of bins
  nbin <- max(1, floor(length(x)/100))
  
  # splits and centres
  if(nbin==1){
    xtholds <- max(x)+1
    xx <- mean(x)
  }else{   
    xtholds <- quantile(x,(1:(nbin-1))/nbin)
    xx <- (c(min(x),xtholds)+c(xtholds,max(x)))/2
  }
  
  # quantiles for bins
  qxx <- matrix(NA,nbin,length(q))
  for(ii in 1:nbin){
      if(nbin==1){
          samp <- y
      }else{
           if(ii == 1){
               samp <- y[x<xtholds[ii]]
           }
           if(ii==nbin){
               samp <- y[x>=xtholds[ii-1]]
           }
           if( ii >1 & ii<nbin){
               samp <- y[x>=xtholds[ii-1] &x<xtholds[ii] ]
           }
       }
      samp <- samp[is.finite(samp)]
      if(length(samp)<3){
          qxx[ii,] <- rep(NA,length(q))      
      }else{
           tmp <- density(samp)
           cdf <- cumsum(tmp$y)/sum(tmp$y)
           qxx[ii,] <- approx(cdf,tmp$x,q)$y
       }  
  }
  # quantiles for chosen values of x
  out <- matrix(NA,length(xout),length(q))
  for(ii in 1:length(q)){
      if(nbin==1){
          out[,ii] <- qxx[,ii]
      }else{
           if( sum(is.finite(qxx[,ii]))  > 1){
               out[,ii] <- approx(xx,qxx[,ii],xout,rule=rule)$y
           }
       }
  }
  
  
  
  return(list(x = xout,y = out,q=q))
}