在分位数不唯一的 R 中分配分位数

Assigning quantiles in R where quantiles are not unique

x 为数值向量,non-negative 数据(主要是 < 10)和 qx <- quantile(x, probs = pq),其中 length(pq) 通常是 > length(x) * (3/4) . 我需要 qx 的索引向量,称之为 q_i,其中 x[i] 落在分位数 qx[q_i[i]].

如标题所示,问题是 qx 中可能存在 non-unique 个值,例如如果 x 是 zero-inflated,则多个 0 值分位数,并且可能还有其他重复值。我想通过(a)回收这些等效分位数的索引序列,或(b)随机分配等效分位数的索引来处理这些情况。我想我更喜欢选项 (a),但是任何一个的解决方案都是有用的。

这里是为特定 x[i] 提供确定 q_i[i] 的规则的编辑: 考虑 qx 有一个或多个重复值序列,即对于某些 j 有(是)序列 qx[j:n] 其中 qx[j] == qx[j + 1] == ... == qx[j + n] < qx[j + n + 1]。让k = c(j, j + 1,..., j + n)。然后 q_i[i] <- k[r] 其中 qx[j] <= x[i] <= qx[j + n + 1] 如果 j == 1qx[j] < x[i] <= qx[j + n + 1] 如果 j > 1,并且 r <- m %% (n + 1) 这样 x[i] 就是 mx 中第 1 次出现,其中不等式已得到满足。

注意:根据这条规则,我意识到我在原来的 q_i 中省略了一个 4 - 这已经改变了。

注意:@hodgenovice 就特殊情况提出了一个很好的观点,在这种情况下,严格小于两个分位数的数据值可能会被分组到两个这样的分位数之间的 "bin" 中。我并不特别关心特殊情况,因为例如,如果没有重复的分位数但我们有相同的分位数值,那么这些特殊情况将正确地合并在一起。

我认为有一种有效的方法可以做到这一点 - 我基本上是使用 for 循环完成的,但我正在寻找一种矢量化方法。

我开始尝试使用 cut() ,它当然不允许 non-unique 中断。我发现 this question here 有帮助,因为我发现了 .bincode() 函数,它确实允许 non-unique 中断。但是,它没有 "distributing" 索引的规则 - 它只会使用每个重复分位数的第一个索引。

这个问题的一些示例代码:

x <- c(5.8,  0.0, 16.1,  5.8,  3.5, 13.8,  6.9,  5.8, 11.5,  9.2, 11.5,
       3.5,  0.0,  8.1,  0.0,  4.6,  5.8,  3.5,  0.0, 10.3,  0.0,  0.0,
       3.5, 6.9, 3.5)
pq <- seq(0, 1, length.out = 20)
qx <- quantile(x, pq)

# quantiles for reference, rounded for readability
round(as.numeric(qx), 2)
[1]  0.00  0.00  0.00  0.00  0.18  3.50  3.50  3.50  3.62  5.04  5.80 5.80  5.97
[14] 6.90  7.72  9.14 10.55 11.50 13.19 16.10

q_i <- .bincode(x, qx, include.lowest = TRUE)
q_i
[1] 10  1 19 10  5 19 13 10 17 16 17  5  1 15  1  9 10  5  1 16  1  1  5 13 5

这是我想要的结果,如果 .bincode() 很神奇,我可以说服它做我需要的事情:

在上述情况 (a) 下: (我也编辑了这个,因为我最初缺少 4 的值)

q_i
[1] 10 1 19 11 5 19 13 10 17 16 17 6 2 15 3 9 11 7 4 16 1 2 5 13 6

在场景 (b) 下,它看起来与上面的一样的可能性很小。或者类似的东西:

q_i
[1] 10 1 19 10 6 19 13 11 17 16 17 5 3 15 2 9 11 6 2 16 1 4 5 13 7

请注意,"equivalent" qx 序列的完整向量被回收基本上是在没有替换的情况下采样的。

谢谢!

好的,我有一些代码从你的代码继续到方案 a(回收)下的最终 q_i。我希望它更漂亮一点,但无论如何希望它能有所帮助。

注:
- 这假设 length(x) > length(qx) > length(x)/2.
- 在代码下方的解释中,q_i 指的是问题末尾的值,在任何回收或替换值发生之前。

## Start off with the code provided in the question...
#  1. For each distinct q_i, calculate the number of occurrances, and how far we can recycle it
df <- data.frame(lower=sort(unique(q_i)), freq=as.integer(table(q_i)))
df$upper <- c(df$lower[-1] - df$lower[-nrow(df)], 1) + df$lower - 1
df$upper <- df$upper - as.numeric(df$upper > df$lower & qx[df$upper] < qx[df$upper + 1])

#  2. Identify when there's a (single) number we can't recycle, and identify which position it's in
#     e.g. is it the third time q_i == 10?
df$special_case <- rep(NA, nrow(df))
df$special_case[df$lower < df$upper] <- sapply(df$lower[df$lower < df$upper], function(low) {
                                        bin <- x[q_i==low]
                                        if(length(unique(bin)) > 1) {
                                          return(match(min(bin), bin))} 
                                        else return(NA)})

# 3. For each row of df, get a vector of (possibly recycled) numbers
recycled <- apply(df, 1, function(x) {
  out <- rep(x["lower"]:x["upper"], length.out=x["freq"])

  # This part modifies the vector created to handle the 'special case'
  if(!is.na(x["special_case"])) {
    out[x["special_case"]] <- x["lower"]
    if(x["special_case"] < x["freq"]) {
      out[(x["special_case"]+1):x["freq"]] <- out[x["special_case"]:(x["freq"]-1)]
    }
  }
  return(out)
})

# 3b. Make this follow the same order as q_i
q_i_final <- unlist(recycled)[order(order(q_i))]

q_i_final
[1] 10  1 19 11  5 19 13 10 17 16 17  6  2 15  3  9 11  7  1 16  2  3  5 13  6

基本思想是什么?
对于 q_i 的每个值,我们可以很容易地计算出我们应该回收的数量(如果我们应该回收的话)。我们通常最多可以回收比 q_i 的下一个最大值少一个。然后我们可以使用 rep 创建一个回收向量来替换 q_i 中的内容,例如将四个 10 替换为 10 11 10 11.

还有什么要考虑的?
这个基本思想假设对于 q_i 的每个值,x 的相应值可以全部回收或全部不回收。这是 通常 的情况,但你也可以有一些 q_i 的值,其中 all bar one 可以被回收,即一个 k 这样的x[k] < qx[q_i[k]+1],但一个或多个 j 其中 q_i[j] = q_i[k] 以及 x[j] = qx[q_i[j]+1].

应识别此类 'special' 案例(尽管问题数据中不存在),并且必须注意不要回收此值。


特例更详细

  1. 我们可以对问题数据做一些简单的更改来创建这个案例(见下面的代码)。请注意 x[5] > x[12],但 q_i[5] = q_i[12] = 4。现在,在上述 'basic idea' 下,所有 q_i = 4 的值都将被回收,因此我们将有 q_i_final[12] = 5。这是一个问题,因为我们希望 x[12] 介于 qx[q_i_final[12]]qx[q_i_final[12]+1] 之间,但事实并非如此,因为它严格小于两者。事实证明我们可以回收 q_i = 4 的所有值,除了 x[12].

新代码:

# Code copied from question, changes as follows:
# x[12] changed from 3.5 to 3.4
# x[13] and x[21] changed from 0.0 to 10.0
x <- c(5.8,  0.0, 16.1,  5.8,  3.5, 13.8,  6.9,  5.8, 11.5,  9.2, 11.5,
       3.4,  10.0,  8.1,  0.0,  4.6,  5.8,  3.5,  0.0, 10.3,  10.0,  0.0,
       3.5, 6.9, 3.5)
pq <- seq(0, 1, length.out = 20)
qx <- quantile(x, pq)
q_i <- .bincode(x, qx, include.lowest = T, right=T)

q_i
[1]  8  1 19  8  4 19 12  8 17 14 17  4 15 13  1  8  8  4  1 16 15  1  4 12  4

此代码基于@hodgenovice 的回答,但未考虑特殊情况。

它有一个附加条件,可以正确回收第一个重复分位数序列的值。这是我在问题中的一个错误,我最初从我想要的答案中省略了 4q_i,但它应该是为分配了 q_i 的数据值回收的索引之一1 来自 .bincode().

df <- data.frame(lower=sort(unique(q_i)), freq=as.integer(table(q_i)))
df$upper <- c(df$lower[-1] - df$lower[-nrow(df)], 1) + df$lower - 1
# want to omit this adjustment if the first quantile is also the first
#   duplicate, to follow rule described in question edit
ub <- df$lower != 1
df$upper[ub] <- df$upper[ub] - as.numeric(df$upper[ub] > df$lower[ub] & 
                  qx[df$upper[ub]] < qx[df$upper[ub] + 1])

recycled <- apply(df, 1, function(x) {
  out <- rep(x["lower"]:x["upper"], length.out=x["freq"])

  return(out)
})

q_i_final <- unlist(recycled)[order(order(q_i))]