tags:

views:

3647

answers:

3

I've been mystified by the R quantile function all day.

I have an intuitive notion of how quantiles work, and an M.S. in stats, but boy oh boy, the documentation for it is confusing to me.

From the docs:

Q[i](p) = (1 - gamma) x[j] + gamma x[j+1],

I'm with it so far. For a type i quantile, it's an interpolation between x[j] and x [j+1], based on some mysterious constant gamma

where 1 <= i <= 9, (j-m)/n <= p < (j-m+1)/ n, x[j] is the jth order statistic, n is the sample size, and m is a constant determined by the sample quantile type. Here gamma depends on the fractional part of g = np+m-j.

So, how calculate j? m?

For the continuous sample quantile types (4 through 9), the sample quantiles can be obtained by linear interpolation between the kth order statistic and p(k):

p(k) = (k - alpha) / (n - alpha - beta + 1), where α and β are constants determined by the type. Further, m = alpha + p(1 - alpha - beta), and gamma = g.

Now I'm really lost. p, which was a constant before, is now apparently a function.

So for Type 7 quantiles, the default...

Type 7

p(k) = (k - 1) / (n - 1). In this case, p(k) = mode[F(x[k])]. This is used by S.

Anyone want to help me out? In particular I'm confused by the notation of p being a function and a constant, what the heck m is, and now to calculate j for some particular p.

I hope that based on the answers here, we can submit some revised documentation that better explains what is going on here.

quantile.R source code or type: quantile.default

+1  A: 

There are various ways of computing quantiles when you give it a vector, and don't have a known CDF.

Consider the question of what to do when your observations don't fall on quantiles exactly.

The "types" are just determining how to do that. So, the methods say, "use a linear interpolation between the k-th order statistic and p(k)".

So, what's p(k)? One guy says, "well, I like to use k/n". Another guy says, "I like to use (k-1)/(n-1)" etc. Each of these methods have different properties that are better suited for one problem or another.

The \alpha's and \beta's are just ways to parameterize the functions p. In one case, they're 1 and 1. In another case, they're 3/8 and -1/4. I don't think the p's are ever a constant in the documentation. They just don't always show the dependency explicitly.

See what happens with the different types when you put in vectors like 1:5 and 1:6.

(also note that even if your observations fall exactly on the quantiles, certain types will still use linear interpolation).

Baltimark
+1  A: 

Mathematica's quantile function is likely similar and its documentation is possibly more elucidating:

http://reference.wolfram.com/mathematica/ref/Quantile.html

Here is general information about quantile functions:

http://en.wikipedia.org/wiki/Quantile_function

dreeves
+5  A: 

You're understandably confused. That documentation is terrible. I had to go back to the paper its based on (Hyndman, R.J.; Fan, Y. (November 1996). "Sample Quantiles in Statistical Packages". American Statistician 50 (4): 361–365. doi:10.2307/2684934) to get an understanding. Let's start with the first problem.

where 1 <= i <= 9, (j-m)/n <= p < (j-m+1)/ n, x[j] is the jth order statistic, n is the sample size, and m is a constant determined by the sample quantile type. Here gamma depends on the fractional part of g = np+m-j.

The first part comes straight from the paper, but what the documentation writers omitted was that j = int(pn+m). This means Q[i](p) only depends on the two order statistics closest to being p fraction of the way through the (sorted) observations. (For those, like me, who are unfamiliar with the term, the "order statistics" of a series of observations is the sorted series.)

Also, that last sentence is just wrong. It should read

Here gamma depends on the fractional part of np+m, g = np+m-j

As for m that's straightforward. m depends on which of the 9 algorithms was chosen. So just like Q[i] is the quantile function, m should be considered m[i]. For algorithms 1 and 2, m is 0, for 3, m is -1/2, and for the others, that's in the next part.

For the continuous sample quantile types (4 through 9), the sample quantiles can be obtained by linear interpolation between the kth order statistic and p(k):

p(k) = (k - alpha) / (n - alpha - beta + 1), where α and β are constants determined by the type. Further, m = alpha + p(1 - alpha - beta), and gamma = g.

This is really confusing. What the documentation calls p(k) is not the same as the p from before. p(k) is the plotting position. In the paper, the authors write it as pk, which helps. Especially since in the expression for m, the p is the original p, and the m = alpha + p * (1 - alpha - beta). Conceptually, for algorithms 4-9, the points (pk, x[k]) are interpolated to get the solution (p, Q[i](p)). Each algorithm only differs in the algorithm for the pk.

As for the last bit, R is just stating what S uses.

The original paper gives a list of 6 "desirable properties for a sample quantile" function, and states a preference for #8 which satisfies all by 1. #5 satisfies all of them, but they don't like it on other grounds (it's more phenomenological than derived from principles). #2 is what non-stat geeks like myself would consider the quantiles and is what's described in wikipedia.

BTW, in response to dreeves answer, Mathematica does things significantly differently. I think I understand the mapping. While Mathematica's is easier to understand, (a) it's easier to shoot yourself in the foot with nonsensical parameters, and (b) it can't do R's algorithm #2. (Here's Mathworld's Quantile page, which states Mathematica can't do #2, but gives a simpler generalization of all the other algorithms in terms of four parameters.)

AFoglia
Thank you for actually answering my question :) That was a serious amount of detective work.
Gregg Lind
No problem. I'm trying to write a quantile function for Python/Numpy for our group, which lead me to this question. When I eventually found the answer, I figured I'd share.
AFoglia
I wrote the quantile() function and the associated help file and submitted it to the R core team in August 2004 (replacing the previous versions). I've just checked and all of these errors were caused by my help file being changed after I submitted it. (I am responsible for the use of p and p[k] though.) I'd never noticed it as I assumed my file would be left untouched. I'll see if I can get the help file fixed for R 2.10.0.
Rob Hyndman
@AFoglia. I've put a proposed new help file at http://www.robjhyndman.com/quantile.html. Comments before I submit to Rcore?
Rob Hyndman
The new one is much better. I have a minor suggestion to add the definitions of gamma for methods one through three, although that might not be necessary for the statistical-knowledgable R audience. Otherwise, it looks great.
AFoglia
Just to complete this discussion, the new help file is now part of base Rv2.10.0.
Rob Hyndman