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
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]
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]
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 = "")
o.pr[p.ok] <- qs
names(o.pr) <- rep("", length(o.pr))
names(o.pr)[p.ok] <- names(qs)
o.pr
\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.}
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.
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)
This archive was generated by hypermail 2.1.8 : Wed 03 Nov 2004 - 22:45:03 EST