Assigning quantiles in R where quantiles are not unique

Question

Let x be a vector of numeric, non-negative data (mostly < 10) and qx <- quantile(x, probs = pq), and where length(pq) is typically > length(x) * (3/4). I am in need of a vector of indices of qx, call it q_i, where x[i] falls in the quantile qx[q_i[i]].

The catch, as the title indicates, is that there may be non-unique values present in qx, e.g. multiple 0-valued quantiles if x is zero-inflated, and potentially other duplicate values. I would like to handle these cases by either (a) recycling the sequence of indices of these equivalent quantiles, or (b) randomly assigning the indices of equivalent quantiles. I think I would prefer option (a), but a solution for either would be useful.

Here is an edit to provide the rules for determining q_i[i] for a particular x[i]: Consider that qx has one or more sequences of duplicate values, i.e. for some j there is (are) sequence(s) qx[j:n] where qx[j] == qx[j + 1] == ... == qx[j + n] < qx[j + n + 1]. Let k = c(j, j + 1,..., j + n). Then q_i[i] <- k[r] where qx[j] <= x[i] <= qx[j + n + 1] if j == 1 or qx[j] < x[i] <= qx[j + n + 1] if j > 1, and where r <- m %% (n + 1) such that x[i] is the m-th occurrence in x where the inequality has been satisfied.

NOTE: based on this rule, I realized I omitted a 4 in my original q_i - this has been changed.

NOTE: @hodgenovice brought up a good point regarding special cases where data values that are strictly smaller than two quantiles may be grouped into the "bin" between two such quantiles. I am not particularly concerned with the special case because, if for example there were no duplicate quantiles but we had the same quantile values, those special cases would correctly be binned together.

I'm thinking there is an efficient way to do this - I have essentially done this using a for loop but I am looking for a vectorized approach.

I started trying to work with cut() which of course doesn't allow non-unique breaks. I found this question here which kind of helped, in that I discovered the .bincode() function, which does allow non-unique breaks. However, it has no rule for "distributing" the indices - it would only use the index of the first of each duplicated quantile value.

Some example code for this problem:

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

Here are the results I would be looking for, if .bincode() was magic and I could talk it into doing what I need:

Under scenario (a) above: (I edited this too, as I was originally missing a value of 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

Under scenario (b), it could, with low probability, look the same as directly above. Or something like:

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

Note here that the full vectors of "equivalent" qx sequences that get recycled are essentially sampled without replacement.

Thanks!

Answer 1

Okay, I've got some code which continues from yours to get to the final q_i under scenario a (recycling). I wish it was a bit prettier, but hope it helps anyway.

Note:
- This assumes length(x) > length(qx) > length(x)/2.
- In the explanation below the code, q_i refers to the value at the end of the question, before any recycling or replacement of values has taken place.

## 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

What's the basic idea?
For each value of q_i, we can fairly easily calculate the number we should recycle up to (if we should recycle at all). We can typically recycle up to one less than the next largest value of q_i. We can then use rep to create a recycled vector to replace what's in q_i e.g. to replace the four 10s with 10 11 10 11.

What else is there to consider?
This basic idea assumes that for each value of q_i, the corresponding value(s) of x can either be all recycled or all not. This is usually the case, but you can also have some value of q_i where all bar one can be recycled i.e. one k such that x[k] < qx[q_i[k]+1], but one or more j where q_i[j] = q_i[k] and also x[j] = qx[q_i[j]+1].

Such 'special' cases (although not present in the question data) should be identified, and care must be taken that this value is not also recycled.

---

Special case in more detail

1. We can make some simple changes to the question data to create this case (see code below). Note that x[5] > x[12], but q_i[5] = q_i[12] = 4. Now, under the 'basic idea' described above, all values where q_i = 4 would be recycled, so we would have q_i_final[12] = 5. This is a problem because we want x[12] to be between qx[q_i_final[12]] and qx[q_i_final[12]+1], which isn't the case since it's strictly less than both. It turns out we are fine to recycle all the values where q_i = 4, except for x[12].

New Code:

# 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

Answer 2

This code is based off of @hodgenovice's answer, but does not consider the special case.

It has an additional condition that correctly recycles the values for the first sequence of duplicated quantiles. This was an error on my part in the question, I originally omitted a q_i of 4 from my desired answer, but it should be one of the indices recycled for data values assigned a q_i of 1 by .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))]