解释 R 中的分位数()函数

我一整天都被分位函数搞糊涂了。

我对分位数的工作原理有一个直观的概念，并且在统计学上有一个理学硕士学位，但是天哪，天哪，它的文档让我很困惑。

来自文件:

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

到目前为止我还是这么想的。对于 我分位数类型，它是基于某个神秘常数 Γ的 x [ j ]和 x [ j + 1]之间的插值

其中1 < = i < = 9，(j-m)/n < = p < (j-m + 1)/n，x [ j ]是第 j 阶 统计量，n 是样本量，m 是样本量 是由样品决定的常数 分位数类型，这里伽玛取决于 G = np + m-j 的分数部分。

那么，如何计算 j? m？

对于连续样本分位数 类别(4至9)、样本 分位数可以通过线性 K 阶之间的插值 统计和 p (k) :

P (k) = (k-alpha)/(n-alpha-beta + 1) , 其中 α 和 β 是常数 此外，m = alpha + p (1 - α-β)和 γ = g。

现在我真的迷路了。 p 以前是个常数，现在明显是个函数。

所以对于类型7的分位数，默认的..。

类型7

P (k) = (k-1)/(n-1)在这种情况下，p (k) = mode [ F (x [ k ])]。

有人想帮我吗？特别是 p 是一个函数和一个常数的表示法，我是什么鬼东西，现在要计算某个特定 P的 j，这让我感到困惑。

我希望基于这里的答案，我们可以提交一些经过修订的文件，更好地解释这里发生的情况。

来源代码 或类型: Quantile.default

34951 次浏览
小开

当你给它一个向量并且没有已知的 CDF 时，有很多计算分位数的方法。

考虑一下这样的问题: 当你的观察结果不完全落在分位数上时，你该怎么办。

“类型”只是决定如何做到这一点。因此，这些方法说，“使用 k 阶统计量和 p (k)之间的一个线性插值。”。

P (k)是多少？一个人说“我喜欢用 k/n”。另一个家伙说，“我喜欢用(k-1)/(n-1)”等等。每个方法都有不同的属性，这些属性更适合于这个或那个问题。

Alpha 和 beta 只是参数化函数 p 的方法。在一种情况下，他们是1和1。在另一种情况下，他们是3/8和 -1/4。我不认为 P 在文档中是一个常量。它们只是不总是显式地显示依赖关系。

看看当你放入向量1:5和1:6时，不同的类型会发生什么。

(还要注意，即使你的观察恰好落在分位数上，某些类型仍然会使用线性插值)。

小开
最佳答案

你很困惑是可以理解的。那些文件太可怕了。我不得不回到报纸的基础上(Hyndman，R.J. ; Fan，Y。(1996年11月)。“统计包中的抽样分位数”。美国统计学家50(4) : 361-365.Doi: 10.2307/2684934)了解。让我们从第一个问题开始。

其中1 < = i < = 9，(j-m)/n < = p < (j-m + 1)/n，x [ j ]是 jth 阶统计量，n 是样本量，m 是由样本分位数类型确定的常数。这里 γ 取决于 g = np + m-j 的分数部分。

第一部分直接来自论文，但是文档作者忽略了 j = int(pn+m)。这意味着 Q[i](p)只依赖于最接近 p分数的两阶统计量通过(排序)观测的方式。(对于像我这样不熟悉这个术语的人来说，一系列观察结果的“顺序统计”就是排序后的序列。)

还有，最后一句话也不对，应该改成

这里 γ 取决于 np + m 的分数部分，g = np + m-j

至于 m，这很简单。m取决于选择了9种算法中的哪一种。因此，就像 Q[i]是分位函数一样，m也应该被认为是 m[i]。对于算法1和算法2，m是0，对于算法3，m是 -1/2，对于其他算法，这在下一部分。

就连续样本分位数类别(4至9)而言，样本分位数可用 kth 阶统计量与 p (k)之间的线性插值计算:

P (k) = (k-alpha)/(n-alpha-beta + 1) ，其中 α 和 β 是由类型决定的常数。此外，m = alpha + p (1-alpha-beta)和 γ = g。

这真的很让人困惑。文档调用的 p(k)与之前的 p不同。p(k)就是 p6。在论文中，作者把它写成 pp7，这有所帮助。特别是在 m的表达中，p是原来的 p，而 m = alpha + p * (1 - alpha - beta)是原来的 p。从概念上讲，对于算法4-9，点(pp7，p1)是插值得到的解决方案(p，p3)。每种算法只是在 pp7的算法中有所不同。

至于最后一点，R 只是说明 S 使用什么。

原文给出了6个“样本分位数的期望性质”函数的列表，并指出了满足所有1的 # 8的偏好。第五条满足了他们所有的要求，但是他们不喜欢其他的理由(它更多的是现象学的而不是从原则中衍生出来的)。第二点是像我这样的非统计学极客会考虑的分位数，也是维基百科所描述的。

顺便说一句，对于 Dreeves 回答，Mathematica 做事情的方式明显不同。我想我明白地图的意思了。虽然 Mathematica 的算法更容易理解，(a)用无意义的参数打自己的脚更容易，(b)它不能做 R 的算法 # 2。(这里是 数学世界的分位数页面，它指出 Mathematica 不能执行 # 2，但是用四个参数对所有其他算法进行了更简单的推广。)

小开

我相信 R 帮助文档在@RobHyndman 的评论中注明了修改后是清晰的，但是我发现它有点让人喘不过气来。我张贴这个答案，以防它帮助有人快速通过选项和他们的假设。

为了掌握 quantile(x, probs=probs)，我想查看一下源代码。这也是比我在 R 中预期的更棘手，所以我实际上只是从一个 Github 回购抓取它看起来足够近运行。我对默认的(类型7)行为感兴趣，所以我注释了一些，但是没有对每个选项做同样的事情。

您可以在代码中看到“ type 7”方法是如何一步一步地进行插值的，我还添加了几行代码来打印一些重要的值。

quantile.default <-function(x, probs = seq(0, 1, 0.25), na.rm = FALSE, names = TRUE
, type = 7, ...){
if(is.factor(x)) { #worry about non-numeric data
if(!is.ordered(x) || ! type %in% c(1L, 3L))
stop("factors are not allowed")
lx <- levels(x)
} else lx <- NULL
if (na.rm){
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 #this is to deal with rounding things sensibly
if (any((p.ok <- !is.na(probs)) & (probs < -eps | probs > 1+eps)))
stop("'probs' outside [0,1]")


#####################################
# here is where terms really used in default type==7 situation get defined


n <- length(x) #how many observations are in sample?


if(na.p <- any(!p.ok)) { # set aside NA & NaN
o.pr <- probs
probs <- probs[p.ok]
probs <- pmax(0, pmin(1, probs)) # allow for slight overshoot
}


np <- length(probs) #how many quantiles are you computing?


if (n > 0 && np > 0) { #have positive observations and # quantiles to compute
if(type == 7) { # be completely back-compatible


index <- 1 + (n - 1) * probs #this gives the order statistic of the quantiles
lo <- floor(index)  #this is the observed order statistic just below each quantile
hi <- ceiling(index) #above
x <- sort(x, partial = unique(c(lo, hi))) #the partial thing is to reduce time to sort,
#and it only guarantees that sorting is "right" at these order statistics, important for large vectors
#ties are not broken and tied elements just stay in their original order
qs <- x[lo] #the values associated with the "floor" order statistics
i <- which(index > lo) #which of the order statistics for the quantiles do not land on an order statistic for an observed value


#this is the difference between the order statistic and the available ranks, i think
h <- (index - lo)[i] # > 0  by construction
##      qs[i] <- qs[i] + .minus(x[hi[i]], x[lo[i]]) * (index[i] - lo[i])
##      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]] # This is the interpolation step: assemble the estimated quantile by removing h*low and adding back in h*high.
# h is the arithmetic difference between the desired order statistic amd the available ranks
#interpolation only occurs if the desired order statistic is not observed, e.g. .5 quantile is the actual observed median if n is odd.
# This means having a more extreme 99th observation doesn't matter when computing the .75 quantile




###################################
# print all of these things


cat("floor pos=", c(lo))
cat("\nceiling pos=", c(hi))
cat("\nfloor values= ", c(x[lo]))
cat( "\nwhich floors not targets? ", c(i))
cat("\ninterpolate between ", c(x[lo[i]]), ";", c(x[hi[i]]))
cat( "\nadjustment values= ", c(h))
cat("\nquantile estimates:")


}else if (type <= 3){## 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)
h <- switch(type,
(nppm > j),     # type 1
((nppm > j) + 1)/2, # type 2
(nppm != j) | ((j %% 2L) == 1L)) # 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 (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)
h <- nppm - j


if(any(sml <- abs(h) < fuzz)) h[sml] <- 0


x <- sort(x, partial =
unique(c(1, j[j>0L & j<=n], (j+1)[j>0L & j<n], n))
)
x <- c(x[1L], x[1L], x, x[n], x[n])
## h can be zero or one (types 1 to 3), and infinities matter
####        qs <- (1 - h) * x[j + 2] + h * x[j + 3]
## also h*x might be invalid ... e.g. Dates and ordered factors
qs <- x[j+2L]
qs[h == 1] <- x[j+3L][h == 1]
other <- (0 < h) & (h < 1)
if(any(other)) qs[other] <- ((1-h)*x[j+2L] + h*x[j+3L])[other]


}
} else {
qs <- rep(NA_real_, np)}


if(is.character(lx)){
qs <- factor(qs, levels = seq_along(lx), labels = lx, ordered = TRUE)}
if(names && np > 0L) {
names(qs) <- format_perc(probs)
}
if(na.p) { # do this more elegantly (?!)
o.pr[p.ok] <- qs
names(o.pr) <- rep("", length(o.pr)) # suppress <NA> names
names(o.pr)[p.ok] <- names(qs)
o.pr
} else qs
}


####################


# fake data
x<-c(1,2,2,2,3,3,3,4,4,4,4,4,5,5,5,5,5,5,5,5,5,6,6,7,99)
y<-c(1,2,2,2,3,3,3,4,4,4,4,4,5,5,5,5,5,5,5,5,5,6,6,7,9)
z<-c(1,2,2,2,3,3,3,4,4,4,4,4,5,5,5,5,5,5,5,5,5,6,6,7)


#quantiles "of interest"
probs<-c(0.5, 0.75, 0.95, 0.975)


# a tiny bit of illustrative behavior
quantile.default(x,probs=probs, names=F)
quantile.default(y,probs=probs, names=F) #only difference is .975 quantile since that is driven by highest 2 observations
quantile.default(z,probs=probs, names=F) # This shifts everything b/c now none of the quantiles fall on an observation (and of course the distribution changed...)... but
#.75 quantile is stil 5.0 b/c the observations just above and below the order statistic for that quantile are still 5. However, it got there for a different reason.


#how does rescaling affect quantile estimates?
sqrt(quantile.default(x^2, probs=probs, names=F))
exp(quantile.default(log(x), probs=probs, names=F))

提交答案