quantile() tweak: using *relative* fuzz before rounding, also for type = 1,2,3; Rd: more references

git-svn-id: https://svn.r-project.org/R/trunk@87746 00db46b3-68df-0310-9c12-caf00c1e9a41

Author: maechler

doc/NEWS.Rd

@@ -593,6 +593,10 @@
       when it tars the <pkg> directory, fixing both \PR{18432} and
       \PR{18434}, for \command{vim} and \command{GNU Global}
       \command{emacs} users, thanks to \I{Dirk Eddelbuettel}'s patch.
+
+      \item \code{quantile()}'s default method gets an option \code{fuzz}
+      to protect against floating point inaccuracy before rounding, thus
+      fixing \PR{15811} and, en passant, \PR{17683}.
     }
   }
 }

src/library/stats/R/quantile.R

@@ -1,7 +1,7 @@
 #  File src/library/stats/R/quantile.R
 #  Part of the R package, https://www.R-project.org
 #
-#  Copyright (C) 1995-2022 The R Core Team
+#  Copyright (C) 1995-2025 The R Core Team
 #
 #  This program is free software; you can redistribute it and/or modify
 #  it under the terms of the GNU General Public License as published by
@@ -23,8 +23,11 @@ quantile.POSIXt <- function(x, ...)
 
 quantile.default <-
     function(x, probs = seq(0, 1, 0.25), na.rm = FALSE, names = TRUE,
-             type = 7, digits = 7, ...)
+             type = 7L, digits = 7,
+             fuzz = if(type == 7L) 0 else 4 * .Machine$double.eps, ...)
 {
+    if(length(type) != 1L || is.na(type) || !is.numeric(type) || !any(type == 1:9))
+        stop("'type' must be an integer in 1..9")
     if(is.factor(x)) {
 	if(is.ordered(x)) {
 	   if(!any(type == c(1L, 3L)))
@@ -41,16 +44,19 @@ quantile.default <-
 	x <- x[!is.na(x)]
     else if (anyNA(x))
 	stop("missing values and NaN's not allowed if 'na.rm' is FALSE")
-    eps <- 100*.Machine$double.eps
+    eps <- 100*.Machine$double.eps # allow for slight overshoot
     if (any((p.ok <- !is.na(probs)) & (probs < -eps | probs > 1+eps)))
 	stop("'probs' outside [0,1]")
+    probs <- pmax(0, pmin(1, probs))
+    stopifnot(is.finite(fuzz), fuzz >= 0, length(fuzz) == 1L)
+    ## need to watch for rounding errors --> use *relative* fuzz
+    floorF <- function(np) floor(np * (1+fuzz))
     n <- length(x)
-    probs <- pmax(0, pmin(1, probs)) # allow for slight overshoot
     np <- length(probs)
     {
-        if(type == 7) { # be completely back-compatible
+        if(type == 7) { # be completely back-compatible (=> default fuzz := 0)
             index <- 1 + max(n - 1, 0) * probs
-            lo <- floor(index)
+            lo <- floorF(index)
             hi <- ceiling(index)
             x <- sort(x, partial = if(n == 0) numeric() else unique(c(lo, hi)[p.ok]))
             qs <- x[lo]
@@ -60,28 +66,26 @@ quantile.default <-
 ##	    qs[i] <- ifelse(h == 0, qs[i], (1 - h) * qs[i] + h * x[hi[i]])
 	    qs[i] <- (1 - h) * qs[i] + h * x[hi[i]]
         } else {
-            if (type <= 3) {
+            if (type <= 3L) {
                 ## Types 1, 2 and 3 are discontinuous sample qs.
-                nppm <- if (type == 3) n * probs - .5 # n * probs + m; m = -0.5
-                        else n * probs          # m = 0
-                j <- floor(nppm)
+                nppm <- if (type == 3L) n * probs - .5 # n * probs + m; m = -0.5
+                        else            n * probs      #                m = 0
+                j <- floorF(nppm)
 		h <- switch(type,
 			    !p.ok | (nppm > j),	# type 1
 			    ((nppm > j) + 1)/2, # type 2
 			    !p.ok | (nppm != j) | ((j %% 2L) == 1L)) # type 3
             } else {
                 ## Types 4 through 9 are continuous sample qs.
-                switch(type - 3,
+                switch(type - 3L,
                    {a <- 0; b <- 1},    # type 4
                        a <- b <- 0.5,   # type 5
                        a <- b <- 0,     # type 6
                        a <- b <- 1,     # type 7 (unused here)
                        a <- b <- 1 / 3, # type 8
                        a <- b <- 3 / 8) # type 9
-                ## need to watch for rounding errors here
-                fuzz <- 4 * .Machine$double.eps
-                nppm <- a + probs * (n + 1 - a - b) # n*probs + m
-                j <- floor(nppm + fuzz) # m = a + probs*(1 - a - b)
+                nppm <- a + probs * (n + 1 - a - b) # n*probs + m; m:= a + probs*(1 - a - b)
+                j <- floorF(nppm)
                 h <- nppm - j
                 if(any(sml <- abs(h) < fuzz, na.rm = TRUE)) h[sml] <- 0
             }
@@ -112,7 +116,7 @@ quantile.default <-
 
 ##' Formatting() percentages the same way as quantile(*, names=TRUE).
 ##' Should be exported
-##' (and format.pval() moved to stats; both documented on same page)
+##' NB format.pval() in 'base', needed by print.summary.table() ...
 format_perc <- function(x, digits = max(2L, getOption("digits")),
 			probability = TRUE, use.fC = length(x) < 100, ...)
 {

src/library/stats/man/quantile.Rd

@@ -1,6 +1,6 @@
 % File src/library/stats/man/quantile.Rd
 % Part of the R package, https://www.R-project.org
-% Copyright 1995-2020 R Core Team
+% Copyright 1995-2025 R Core Team
 % Distributed under GPL 2 or later
 
 \name{quantile}
@@ -17,7 +17,9 @@
 quantile(x, \dots)
 
 \method{quantile}{default}(x, probs = seq(0, 1, 0.25), na.rm = FALSE,
-         names = TRUE, type = 7, digits = 7, \dots)
+         names = TRUE, type = 7, digits = 7,
+         fuzz = if(type == 7L) 0 else 4 * .Machine$double.eps,
+         \dots)
 }
 \arguments{
   \item{x}{numeric vector whose sample quantiles are wanted, or an
@@ -36,6 +38,9 @@ quantile(x, \dots)
   \item{digits}{used only when \code{names} is true: the precision to use
     when formatting the percentages.  In \R versions up to 4.0.x, this had
     been set to \code{max(2, getOption("digits"))}, internally.}
+  \item{fuzz}{small non-negative number to protect against rounding errors
+    when \code{j <- \link{floor}(np + m)}, (\eqn{np} \dQuote{a version of}
+    \code{n * probs}, see the formula below), is computed.}
   \item{\dots}{further arguments passed to or from other methods.}
 }
 \details{
@@ -57,6 +62,11 @@ quantile(x, \dots)
   There is a method for the date-time classes (see
   \code{"\link{POSIXt}"}).  Types 1 and 3 can be used for class
   \code{"\link{Date}"} and for ordered factors.
+
+  \code{fuzz := 4 * .Machine$double.eps} has been hard coded and used for
+  \code{type = 4,5,6, 8,9} since \code{type}s were introduced, % r30628, 2004-08-12
+  and is used, since \R 4.5.0, also for the other types (but with default
+  \code{0} for the default \code{type = 7} for back compatibility reasons).
 }
 \section{Types}{
   \code{quantile} returns estimates of underlying distribution quantiles
@@ -139,14 +149,16 @@ quantile(x, \dots)
     }
   }
   Further details are provided in \bibcite{Hyndman and Fan (1996)} who
-  recommended type 8. 
+  recommended type 8.
   The default method is type 7, as used by S and by \R < 2.0.0.
-  \I{Makkonen} argues for type 6, also as already proposed by Weibull in 1939.
+  \bibcite{Makkonen and Pajari (2014)} argue for type 6, also as already
+  proposed by Weibull in 1939.
   The Wikipedia page contains further information about availability of
   these 9 types in software.
 }
 \author{
-  of the version used in \R >= 2.0.0, Ivan Frohne and Rob J Hyndman.
+  of the version used in \R >= 2.0.0, Ivan Frohne and Rob J Hyndman;
+  tweaks, notably use of \code{fuzz}, by the R Core Team.
 }
 \references{
   Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988)
@@ -157,13 +169,16 @@ quantile(x, \dots)
   packages, \emph{American Statistician} \bold{50}, 361--365.
   \doi{10.2307/2684934}.
 
+  Eric Langford (2006)  Quartiles in Elementary Statistics, \emph{Journal
+  of Statistics Education} \bold{14} 3; \doi{10.1080/10691898.2006.11910589}
+
+  Makkonen, L. and Pajari, M. (2014)  Defining Sample Quantiles by the True
+  Rank Probability, \emph{Journal of Probability and Statistics; Hindawi Publ.Corp.}
+  \doi{10.1155/2014/326579}
+
   Wicklin, R. (2017) Sample quantiles: A comparison of 9 definitions; SAS Blog.
   \url{https://blogs.sas.com/content/iml/2017/05/24/definitions-sample-quantiles.html}
 
-  %% Makkonen, L. and Pajari, M. (2014)  Defining Sample Quantiles by the True
-  %% Rank Probability. \emph{Journal of Probability and Statistics; Hindawi Publ.Corp.}
-  %% \doi{10.1155/2014/326579}
-  %%
   Wikipedia: \url{https://en.wikipedia.org/wiki/Quantile#Estimating_quantiles_from_a_sample}
 }
 \seealso{

tests/reg-tests-1e.R

@@ -1846,6 +1846,14 @@ tools::parseLatex("\\Sexpr{ 1 + {1} }")
 assertErrV(tools::parseLatex("\\begin{foo} abc \\end{bar}"))
 
 
+## quantile.default() needing fuzz PR#15811
+x <- 1:1390
+stopifnot(identical(print(quantile(x, 0.7, type=2)),
+                    c(`70%` = 973.5)))
+## was 973  in R <= 4.4.x
+
+
+
 ## keep at end
 rbind(last =  proc.time() - .pt,
       total = proc.time())