[Rd] Suggestion for quantile.default()

[Rob Hyndman]

I'm not sure who is responsible for quantile(), but I assume they read 
this list. Ivan Frohne and I have produced a revision of the 
quantile.default() function which enables the computation of alternative 
sample quantile definitions. The code and .Rd file are attached.

This enables the user to produce quantiles that are equivalent to those 
in various statistics package. There is a type argument that allows one 
to choose between the various sample quantile methods. type=7 gives 
identical results to the current R function quantile.default(). In our 
revised function, type=8 is the default following the recommendation of 
Hyndman and Fan (American Statistician, 1996).

We suggest the attached function replaces the current quantile.default() 
function in R as it provides additional functionality without increasing 
computation time or affecting ease-of-use. If backwards-compatibility is 
important, we are happy to set type=7 as the default. However, we prefer 
type=8. For moderate to large sample sizes, the difference is negligble.

The code and associated documentation is as close as possible to what 
already exists.

Regards,
Rob Hyndman


__________________________________________________
Professor Rob J Hyndman
Department of Econometrics & Business Statistics
Monash University, VIC 3800, Australia.
http://www-personal.buseco.monash.edu.au/~hyndman/
-------------- next part --------------
quantile.default <- function(x, probs = seq(0, 1, 0.25), na.rm = FALSE,
                           names = TRUE, type = 8, ...)                                                                             
{
    if (na.rm) 
        x <- x[!is.na(x)]
    else if (any(is.na(x)))
        stop("Missing values and NaN's not allowed if `na.rm' is FALSE")
    if (any((p.ok <- !is.na(probs)) & (probs < 0 | probs > 1))) 
        stop("probs outside [0,1]")
    if (na.p <- any(!p.ok)) {
        o.pr <- probs
        probs <- probs[p.ok]
    }
    type <- as.integer(type)   
    if (is.na(type) || (type < 1 | type > 9))
        stop("type outside range [1,9]")
    np <- length(probs)
    n <- length(x)
    if (n > 0 && np > 0) {                              
        if (type <= 3) {
            ## Types 1, 2 and 3 are discontinuous sample qs.
            if (type == 3)
                nppm <- n * probs - .5          # n * probs + m; m = -0.5
            else 
                nppm <- n * probs               # m = 0
            j <- floor(nppm) 
            switch(type,
                        h <- ifelse(nppm > j, 1, 0),        # type 1
                        h <- ifelse(nppm > j, 1, 0.5),      # type 2
                        h <- ifelse((nppm == j) &&
                             ((j %% 2) == 0), 0, 1))        # type 3
        }
        else {
            ## Types 4 through 9 are continuous sample qs.
            switch(type - 3,
                            {a <- 0; b <- 1},               # type 4
                            a <- b <- 0.5,                  # type 5
                            a <- b <- 0,                    # type 6
                            a <- b <- 1,                    # type 7
                            a <- b <- 1 / 3,                # type 8
                            a <- b <- 3 / 8)                # type 9
            nppm <- a + probs * (n + 1 - a - b) # n*probs + m
            j <- floor(nppm)                    # m = a + probs*(1 - a - b)
            h <- nppm - j 
        }
        x <- sort(x, partial = unique(c(1,j[j>0 & j<=n],(j+1)[j>0 & j<n],n)))
        x <- c(x[1], x[1], x, x[n], x[n])
        qs <- (1 - h) * x[j + 2] + h * x[j + 3]
    }
    else {
        qs <- rep(as.numeric(NA), np)
    }
    if (names && np > 0) {
        dig <- max(2, getOption("digits"))
        names(qs) <- paste(if (np < 100) 
            formatC(100 * probs, format = "fg", wid = 1, digits = dig)
        else format(100 * probs, trim = TRUE, digits = dig), 
            "%", sep = "")
    }
    if (na.p) {
        o.pr[p.ok] <- qs
        names(o.pr) <- rep("", length(o.pr))
        names(o.pr)[p.ok] <- names(qs)
        o.pr
    }
    else qs
}
-------------- next part --------------
\name{quantile}
\title{Sample Quantiles}
\alias{quantile}
\alias{quantile.default}
\description{
    The generic function \code{quantile} produces sample quantiles corresponding to the
    given probabilities. The smallest observation corresponds to a probability of
    0 and the largest to a probability of 1.
}
\usage{
quantile(x, \dots)

\method{quantile}{default}(x, probs = seq(0, 1, 0.25), na.rm = FALSE,
         names = TRUE, type = 8, ...)
}
\arguments{
    \item{x}{numeric vector whose sample quantiles are wanted.}
    \item{probs}{numeric vector of probabilities with values in \eqn{[0,1]}{[0,1]}.}
    \item{na.rm}{logical; if \code{TRUE} any \code{NA} or \code{NaN}
        is removed from \code{x} before the quantiles are computed.  If
        \code{FALSE} the presence of \code{NA} or \code{NaN} in \code{x}
        aborts the function.}
    \item{names}{logical; if \code{TRUE}, the result has a \code{\link[base]{names}} attribute.
    Set to FALSE for efficiency when \code{probs} is long.}
    \item{type}{an integer between 1 and 9 selecting one of the
         nine quantile algorithms detailed below to be used.}
    \item{\dots}{further arguments passed to or from other methods}
}
\value{
A vector of the same length as \code{probs} is returned;
if \code{names = TRUE}, it has a \code{\link[base]{names}} attribute.

\code{\link[base]{NA}} and \code{\link[base]{NaN}} values in \code{probs} are propagated to the result.
}
\details{
\code{quantile} returns estimates of underlying distribution quantiles
based on one or two order statistics from the supplied elements in \code{x} at
probabilities in \code{probs}.  One of the nine quantile algorithms discussed in
Hyndman and Fan (1996), selected by \code{type}, is employed.

Sample quantiles of type \eqn{i} are defined by
\deqn{Q_{i}(p) = (1 - \gamma)x_{j} + \gamma x_{j+1}}{Q[i](p) = (1 - gamma) x[j] + gamma x[j+1],} where
\eqn{1 \le i \le 9}{1 <= i <= 9}, \eqn{  }
\eqn{\frac{j - m}{n} \le p < \frac{j - m + 1}{n}}{(j - m) / n <= p < (j - m + 1) / n}, \eqn{  }
\eqn{x_{j}}{x[j]} is the \eqn{j}th order statistic,
\eqn{n} is the sample size, and \eqn{m} is a constant determined by
the sample quantile type.
Here \eqn{\gamma}{gamma} depends on
\eqn{j = \lfloor np + m\rfloor}{j = floor(np + m),} and
\eqn{g = np + m - j.}

For the continuous sample quantile types (4 through 9), the sample quantiles can be obtained
by linear interpolation between the \eqn{k}th order statistic and \eqn{p(k)}:
\deqn{p(k) = \frac{k - \alpha} {n - \alpha - \beta + 1}}{p(k) = (k - alpha) / (n - alpha - beta + 1),}
where \eqn{\alpha}{alpha} and \eqn{\beta}{beta} are constants determined by the type.
Further,
\eqn{m = \alpha + p \left( 1 - \alpha - \beta \right)}{m = alpha + p(1 - alpha - beta),}
and \eqn{\gamma = g = np + m - j}{gamma = g = np + m - j.}


            \strong{Discontinuous sample quantile types 1, 2, and 3}

\describe{
\item{Type 1}{Inverse of empirical distribution function.}
\item{Type 2}{Similar to type 1 but with averaging at discontinuities.}
\item{Type 3}{SAS definition: nearest even order statistic.}
}


            \strong{Continuous sample quantile types 4 through 9}

\describe{
\item{Type 4}{\eqn{p(k) = \frac{k}{n}}{p(k) = k / n}.
    That is, linear interpolation of the empirical cdf.}

\item{Type 5}{\eqn{p(k) = \frac{k - 0.5}{n}}{p(k) = (k - 0.5) / n}.
    That is a piecewise linear function where the knots are the values
    midway through the steps of the empirical cdf. This is popular amongst hydrologists.}

\item{Type 6}{\eqn{p(k) = \frac{k}{n + 1}}{p(k) = k / (n + 1)}.
    Thus \eqn{p(k) = \mbox{E}[F(x_{k})]}{p(k) = E[F(x[k])]}. This is used by Minitab and by SPSS.}

\item{Type 7}{\eqn{p(k) = \frac{k - 1}{n - 1}}{p(k) = (k - 1) / (n - 1)}.
    In this case, \eqn{p(k) = \mbox{mode}[F(x_{k})]}{p(k) = mode[F(x[k])]}. This is used by S-Plus.}

\item{Type 8}{\eqn{p(k) = \frac{k - \frac{1}{3}}{n + \frac{1}{3}}}{p(k) = (k - 1/3) / (n + 1/3)}.
    Then \eqn{p(k) \approx \mbox{median}[F(x_{k})]}{p(k) =~ median[F(x[k])]}.
    The resulting quantile estimates are approximately median-unbiased regardless of the
    distribution of \code{x}.}

\item{Type 9}{\eqn{p(k) = \frac{k - \frac{3}{8}}{n + \frac{1}{4}}}{p(k) = (k - 3/8) / (n + 1/4)}.
    The resulting quantile estimates are approximately unbiased if \code{x} is normally distributed.}
}

Hyndman and Fan (1996) recommend type 8 which is the default method for this function.
}
\references{
Hyndman, R. J. and Fan, Y. (1996) Sample quantiles in
statistical packages, \emph{American Statistician}, \bold{50}, 361-365.
}
\seealso{
\code{\link[stats]{ecdf}} for empirical distributions of which quantile is
the \dQuote{inverse}; \code{\link[graphics]{boxplot.stats}} and \code{\link[stats]{fivenum}}
for computing \dQuote{versions} of quartiles, etc.
}
\examples{
quantile(x <- rnorm(1001))# Extremes & Quartiles by default

### Compare different methods
p <- c(0.1,0.5,1,2,5,10,50)/100
quantile(x,  p, type=1)
quantile(x,  p, type=2)
quantile(x,  p, type=3)
quantile(x,  p, type=4)
quantile(x,  p, type=5)
quantile(x,  p, type=6)
quantile(x,  p, type=7)
quantile(x,  p, type=8)
quantile(x,  p, type=9)
}
\keyword{univar}
\author{Ivan Frohne and Rob J Hyndman}