Somme rapide des valeurs d'un vecteur au-dessus de seuils donnés

Question

J'ai un vecteur de valeurs de seuil, thresholds et un autre vecteur, x . J'aimerais créer un nouveau vecteur, par exemple vec_sum de la même longueur que thresholds qui stocke, pour chaque élément de thresholds la somme des valeurs de x plus grand que cet élément.

Quel est le moyen le plus rapide de le faire ? La façon naïve dont je le fais est

vec_sum <- rep(NA,length(thresholds))
for(i in seq_along(thresholds))
{
   vec_sum[i] <- sum(x[x>thresholds[i]])
}

Au cas où ça aiderait, les seuils sont déjà triés.

Answer 1

Voici une autre solution utilisant cumsum :

f1 <- function(v, th){
    v2 <- v[order(v)]
    v2s <- rev(cumsum(rev(v2)))
    return(v2s[findInterval(th, v2) + 1])
}

Voici quelques tests et une comparaison avec l'autre réponse (ainsi que les données de l'exemple) de Ronak :

f2 <- function(x, thresholds){
    if (all(x < thresholds[1])) return(rep(0, length(thresholds)))
    if (all(x > thresholds[length(thresholds)])) return(rep(sum(x), length(thresholds)))
    return(rev(cumsum(rev(tapply(x, 
        findInterval(x, thresholds, left.open = TRUE), sum)[-1]))))
}

test_th <- c(3, 5, 10)
test_x <- c(2, 3, 1, 19, 4, 6, 5, 15, 7:14, 16:18, 20)

vec_sum <- rep(NA,length(test_th))
for(i in seq_along(test_th)) {
    vec_sum[i] <- sum(test_x[test_x>test_th[i]])
}

all(dplyr::near(f1(test_x, test_th), vec_sum))
# [1] TRUE
all(dplyr::near(f2(test_x, test_th), vec_sum))
# [1] TRUE

set.seed(123)
test_x <- rnorm(10000)
test_th <- sort(rnorm(100)) ## f2 requires sorted threshold values

vec_sum <- rep(NA,length(test_th))
for(i in seq_along(test_th)) {
    vec_sum[i] <- sum(test_x[test_x>test_th[i]])
}
all(dplyr::near(f1(test_x, test_th), vec_sum))
# [1] TRUE
all(dplyr::near(f2(test_x, test_th), vec_sum))
# [1] FALSE
# Warning message:
# In x - y : longer object length is not a multiple of shorter object length

library(microbenchmark)
microbenchmark(
    a = f1(test_x, test_th),
    b = f2(test_x, test_th)
)
# Unit: microseconds
#  expr      min       lq      mean   median       uq       max neval
#     a  587.116  682.864  900.3572  694.713  703.726 10647.206   100
#     b 1157.213 1203.063 1260.0663 1223.600 1258.552  2143.069   100

Answer 2

Je ne suis pas sûr que ce soit plus rapide, mais nous pouvons utiliser findInterval pour couper x por thresholds . Nous prenons sum de chaque groupe en utilisant tapply et prendre cumsum à l'envers.

as.integer(rev(cumsum(rev(tapply(x, 
          findInterval(x, thresholds, left.open = TRUE), sum)[-1]))))

Testé sur

thresholds <- c(3, 5, 10)
x <- c(2, 3, 1, 19, 4, 6, 5, 15, 7:14, 16:18, 20) #1:20 in random order
vec_sum <- rep(NA,length(thresholds))

for(i in seq_along(thresholds)) {
  vec_sum[i] <- sum(x[x>thresholds[i]])
}
vec_sum
#[1] 204 195 155

Utilisation de la solution proposée

as.integer(rev(cumsum(rev(tapply(x, 
          findInterval(x, thresholds, left.open = TRUE), sum)[-1]))))
#[1] 204 195 155

---

Expliquer la réponse. findInterval renvoie des groupes où chaque valeur de x appartient à

findInterval(x, thresholds, left.open = TRUE)
#[1] 0 0 0 3 1 2 1 3 2 2 2 2 3 3 3 3 3 3 3 3

Nous utilisons tapply pour obtenir sum de chaque groupe

tapply(x, findInterval(x, thresholds, left.open = TRUE), sum)
#  0   1   2   3 
#  6   9  40 155 

du groupe 0 doivent être exclues car elles sont plus petites que toutes les valeurs de threshold (d'où -1 ). Le groupe 2 doit également contenir la somme du groupe 1 et le groupe 3 doit contenir la somme des groupes 1 et 2. Ainsi, nous rev erse la séquence et prend cumsum

cumsum(rev(tapply(x, findInterval(x, thresholds, left.open = TRUE), sum)[-1]))

#  3   2   1 
#155 195 204 

Pour le remettre en état et le faire correspondre à l'original. threshold nous rev Encore une fois

rev(cumsum(rev(tapply(x, findInterval(x, thresholds, left.open = TRUE), sum)[-1])))
#  1   2   3 
#204 195 155 

---

Cas limites :

S'il y a toutes les valeurs en dessous du seuil ou toutes les valeurs au-dessus du seuil, nous devrons peut-être faire une vérification supplémentaire et retourner ce qui suit.

if (all(x < thresholds[1]))   rep(0, length(thresholds))
if (all(x > thresholds[length(thresholds)])) rep(sum(x), length(thresholds))