Cette question conduit à un nouveau package R:
wrswoRR par défaut de l'échantillonnage sans remplacement à l'aide de
sample.intsemble exiger quadratique moment de l'exécution, par exemple lors de l'utilisation de poids élaboré à partir d'une distribution uniforme. C'est lent pour un échantillon de grande taille. Quelqu'un sait d'une mise en œuvre plus rapide qui serait utilisable à partir dans R? Deux options sont "rejet de l'échantillonnage avec remplacement" (voir cette question sur les stats.sx) et l'algorithme par Wong et Easton (1980) (avec un Python de mise en œuvre dans un StackOverflow réponse).Merci à Ben Bolker pour faire allusion à la fonction C qui est appelé en interne lors de l'
sample.intest appelée avecreplace=Fet non-uniforme des pondérations:ProbSampleNoReplace. En effet, le code affiche les deux imbriquésforboucles (ligne 420 ff derandom.c).Voici le code pour analyser le temps d'exécution de façon empirique:
library(plyr) sample.int.test <- function(n, p) { sample.int(2 * n, n, replace=F, prob=p); NULL } times <- ldply( 1:7, function(i) { n <- 1024 * (2 ** i) p <- runif(2 * n) data.frame( n=n, user=system.time(sample.int.test(n, p), gcFirst=T)['user.self']) }, .progress='text' ) times library(ggplot2) ggplot(times, aes(x=n, y=user/n)) + geom_point() + scale_x_log10() + ylab('Time per unit (s)') # Output: n user 1 2048 0.008 2 4096 0.028 3 8192 0.100 4 16384 0.408 5 32768 1.645 6 65536 6.604 7 131072 26.558
EDIT: Merci à Arun pour préciser que non pondérée d'échantillonnage ne semble pas avoir cette pénalité de performances.
Réponses
Trop de publicités?
Mise à jour:
Un Rcpp de la mise en œuvre de Efraimidis & Spirakis algorithme (merci à @Hemmo, @Dinrem, @krlmlr et @rtlgrmpf):
library(inline)
library(Rcpp)
src <-
'
int num = as<int>(size), x = as<int>(n);
Rcpp::NumericVector vx = Rcpp::clone<Rcpp::NumericVector>(x);
Rcpp::NumericVector pr = Rcpp::clone<Rcpp::NumericVector>(prob);
Rcpp::NumericVector rnd = rexp(x) / pr;
for(int i= 0; i<vx.size(); ++i) vx[i] = i;
std::partial_sort(vx.begin(), vx.begin() + num, vx.end(), Comp(rnd));
vx = vx[seq(0, num - 1)] + 1;
return vx;
'
incl <-
'
struct Comp{
Comp(const Rcpp::NumericVector& v ) : _v(v) {}
bool operator ()(int a, int b) { return _v[a] < _v[b]; }
const Rcpp::NumericVector& _v;
};
'
funFast <- cxxfunction(signature(n = "Numeric", size = "integer", prob = "numeric"),
src, plugin = "Rcpp", include = incl)
# See the bottom of the answer for comparison
p <- c(995/1000, rep(1/1000, 5))
n <- 100000
system.time(print(table(replicate(funFast(6, 3, p), n = n)) / n))
1 2 3 4 5 6
1.00000 0.39996 0.39969 0.39973 0.40180 0.39882
user system elapsed
3.93 0.00 3.96
# In case of:
# Rcpp::IntegerVector vx = Rcpp::clone<Rcpp::IntegerVector>(x);
# i.e. instead of NumericVector
1 2 3 4 5 6
1.00000 0.40150 0.39888 0.39925 0.40057 0.39980
user system elapsed
1.93 0.00 2.03
Ancienne version:
Nous allons tenter quelques approches possibles:
Simple rejet de l'échantillonnage avec remplacement. C'est beaucoup plus simple de la fonction qu' sample.int.rej offert par @krlmlr, c'est à dire la taille de l'échantillon est toujours égale à n. Comme nous allons le voir, il est toujours très rapide en supposant une distribution uniforme pour les poids, mais extrêmement lent dans une autre situation.
fastSampleReject <- function(all, n, w){
out <- numeric(0)
while(length(out) < n)
out <- unique(c(out, sample(all, n, replace = TRUE, prob = w)))
out[1:n]
}
L'algorithme par Wong et Easton (1980). Voici une implémentation de cette version de Python. Il est stable et j'ai peut-être raté quelque chose, mais il est beaucoup plus lent par rapport à d'autres fonctions.
fastSample1980 <- function(all, n, w){
tws <- w
for(i in (length(tws) - 1):0)
tws[1 + i] <- sum(tws[1 + i], tws[1 + 2 * i + 1],
tws[1 + 2 * i + 2], na.rm = TRUE)
out <- numeric(n)
for(i in 1:n){
gas <- tws[1] * runif(1)
k <- 0
while(gas > w[1 + k]){
gas <- gas - w[1 + k]
k <- 2 * k + 1
if(gas > tws[1 + k]){
gas <- gas - tws[1 + k]
k <- k + 1
}
}
wgh <- w[1 + k]
out[i] <- all[1 + k]
w[1 + k] <- 0
while(1 + k >= 1){
tws[1 + k] <- tws[1 + k] - wgh
k <- floor((k - 1) / 2)
}
}
out
}
Rcpp mise en œuvre de l'algorithme par Wong et Easton. Éventuellement, il peut être optimisé encore plus puisque c'est ma première utilisables Rcpp de la fonction, mais de toute façon, il fonctionne bien.
library(inline)
library(Rcpp)
src <-
'
Rcpp::NumericVector weights = Rcpp::clone<Rcpp::NumericVector>(w);
Rcpp::NumericVector tws = Rcpp::clone<Rcpp::NumericVector>(w);
Rcpp::NumericVector x = Rcpp::NumericVector(all);
int k, num = as<int>(n);
Rcpp::NumericVector out(num);
double gas, wgh;
if((weights.size() - 1) % 2 == 0){
tws[((weights.size()-1)/2)] += tws[weights.size()-1] + tws[weights.size()-2];
}
else
{
tws[floor((weights.size() - 1)/2)] += tws[weights.size() - 1];
}
for (int i = (floor((weights.size() - 1)/2) - 1); i >= 0; i--){
tws[i] += (tws[2 * i + 1]) + (tws[2 * i + 2]);
}
for(int i = 0; i < num; i++){
gas = as<double>(runif(1)) * tws[0];
k = 0;
while(gas > weights[k]){
gas -= weights[k];
k = 2 * k + 1;
if(gas > tws[k]){
gas -= tws[k];
k += 1;
}
}
wgh = weights[k];
out[i] = x[k];
weights[k] = 0;
while(k > 0){
tws[k] -= wgh;
k = floor((k - 1) / 2);
}
tws[0] -= wgh;
}
return out;
'
fun <- cxxfunction(signature(all = "numeric", n = "integer", w = "numeric"),
src, plugin = "Rcpp")
Maintenant quelques résultats:
times1 <- ldply(
1:6,
function(i) {
n <- 1024 * (2 ** i)
p <- runif(2 * n) # Uniform distribution
p <- p/sum(p)
data.frame(
n=n,
user=c(system.time(sample.int.test(n, p), gcFirst=T)['user.self'],
system.time(weighted_Random_Sample(1:(2*n), p, n), gcFirst=T)['user.self'],
system.time(fun(1:(2*n), n, p), gcFirst=T)['user.self'],
system.time(sample.int.rej(2*n, n, p), gcFirst=T)['user.self'],
system.time(fastSampleReject(1:(2*n), n, p), gcFirst=T)['user.self'],
system.time(fastSample1980(1:(2*n), n, p), gcFirst=T)['user.self']),
id=c("Base", "Reservoir", "Rcpp", "Rejection", "Rejection simple", "1980"))
},
.progress='text'
)
times2 <- ldply(
1:6,
function(i) {
n <- 1024 * (2 ** i)
p <- runif(2 * n - 1)
p <- p/sum(p)
p <- c(0.999, 0.001 * p) # Special case
data.frame(
n=n,
user=c(system.time(sample.int.test(n, p), gcFirst=T)['user.self'],
system.time(weighted_Random_Sample(1:(2*n), p, n), gcFirst=T)['user.self'],
system.time(fun(1:(2*n), n, p), gcFirst=T)['user.self'],
system.time(sample.int.rej(2*n, n, p), gcFirst=T)['user.self'],
system.time(fastSampleReject(1:(2*n), n, p), gcFirst=T)['user.self'],
system.time(fastSample1980(1:(2*n), n, p), gcFirst=T)['user.self']),
id=c("Base", "Reservoir", "Rcpp", "Rejection", "Rejection simple", "1980"))
},
.progress='text'
)
arrange(times1, id)
n user id
1 2048 0.53 1980
2 4096 0.94 1980
3 8192 2.00 1980
4 16384 4.32 1980
5 32768 9.10 1980
6 65536 21.32 1980
7 2048 0.02 Base
8 4096 0.05 Base
9 8192 0.18 Base
10 16384 0.75 Base
11 32768 2.99 Base
12 65536 12.23 Base
13 2048 0.00 Rcpp
14 4096 0.01 Rcpp
15 8192 0.03 Rcpp
16 16384 0.07 Rcpp
17 32768 0.14 Rcpp
18 65536 0.31 Rcpp
19 2048 0.00 Rejection
20 4096 0.00 Rejection
21 8192 0.00 Rejection
22 16384 0.02 Rejection
23 32768 0.02 Rejection
24 65536 0.03 Rejection
25 2048 0.00 Rejection simple
26 4096 0.01 Rejection simple
27 8192 0.00 Rejection simple
28 16384 0.01 Rejection simple
29 32768 0.00 Rejection simple
30 65536 0.05 Rejection simple
31 2048 0.00 Reservoir
32 4096 0.00 Reservoir
33 8192 0.00 Reservoir
34 16384 0.02 Reservoir
35 32768 0.03 Reservoir
36 65536 0.05 Reservoir
arrange(times2, id)
n user id
1 2048 0.43 1980
2 4096 0.93 1980
3 8192 2.00 1980
4 16384 4.36 1980
5 32768 9.08 1980
6 65536 19.34 1980
7 2048 0.01 Base
8 4096 0.04 Base
9 8192 0.18 Base
10 16384 0.75 Base
11 32768 3.11 Base
12 65536 12.04 Base
13 2048 0.01 Rcpp
14 4096 0.02 Rcpp
15 8192 0.03 Rcpp
16 16384 0.08 Rcpp
17 32768 0.15 Rcpp
18 65536 0.33 Rcpp
19 2048 0.00 Rejection
20 4096 0.00 Rejection
21 8192 0.02 Rejection
22 16384 0.02 Rejection
23 32768 0.05 Rejection
24 65536 0.08 Rejection
25 2048 1.43 Rejection simple
26 4096 2.87 Rejection simple
27 8192 6.17 Rejection simple
28 16384 13.68 Rejection simple
29 32768 29.74 Rejection simple
30 65536 73.32 Rejection simple
31 2048 0.00 Reservoir
32 4096 0.00 Reservoir
33 8192 0.02 Reservoir
34 16384 0.02 Reservoir
35 32768 0.02 Reservoir
36 65536 0.04 Reservoir
Évidemment, nous pouvons rejeter fonction 1980 parce qu'il est plus lent que d' Base dans les deux cas. Rejection simple dans le pétrin trop quand il y a une probabilité unique 0.999 dans le second cas.
Donc il reste Rejection, Rcpp, Reservoir. La dernière étape est de vérifier si les valeurs sont correctes. Pour être sûr à leur sujet, nous utiliserons sample comme un indice de référence (également pour éliminer la confusion sur les probabilités qui n'ont pas à coïncider avec p en raison de l'échantillonnage sans remplacement).
p <- c(995/1000, rep(1/1000, 5))
n <- 100000
system.time(print(table(replicate(sample(1:6, 3, repl = FALSE, prob = p), n = n))/n))
1 2 3 4 5 6
1.00000 0.39992 0.39886 0.40088 0.39711 0.40323 # Benchmark
user system elapsed
1.90 0.00 2.03
system.time(print(table(replicate(sample.int.rej(2*3, 3, p), n = n))/n))
1 2 3 4 5 6
1.00000 0.40007 0.40099 0.39962 0.40153 0.39779
user system elapsed
76.02 0.03 77.49 # Slow
system.time(print(table(replicate(weighted_Random_Sample(1:6, p, 3), n = n))/n))
1 2 3 4 5 6
1.00000 0.49535 0.41484 0.36432 0.36338 0.36211 # Incorrect
user system elapsed
3.64 0.01 3.67
system.time(print(table(replicate(fun(1:6, 3, p), n = n))/n))
1 2 3 4 5 6
1.00000 0.39876 0.40031 0.40219 0.40039 0.39835
user system elapsed
4.41 0.02 4.47
Avis un peu les choses ici. Pour certaines raisons, weighted_Random_Sample renvoie des valeurs incorrectes (je n'ai pas regardé en cela du tout, mais il fonctionne correct en supposant une distribution uniforme). sample.int.rej est très lent de l'échantillonnage répété.
En conclusion, il semble que l' Rcpp est le choix optimal dans le cas d'un échantillonnage répété, tout sample.int.rej est un peu plus rapide, autrement, et aussi plus facile à utiliser.
J'ai décidé de creuser dans certains commentaires et trouvé le Efraimidis & Spirakis livre fascinant (merci à @Hemmo pour trouver la référence). L'idée générale du livre est: créer une clé par la génération aléatoire uniforme de nombre et de l'élévation à la puissance de l'un sur le poids de chaque élément. Ensuite, il vous suffit de prendre les plus hautes valeurs de la clé de votre échantillon. Cela fonctionne avec brio!
weighted_Random_Sample <- function(
.data,
.weights,
.n
){
key <- runif(length(.data)) ^ (1 / .weights)
return(.data[order(key, decreasing=TRUE)][1:.n])
}
Si vous réglez".n' étant la longueur de '.de données " (ce qui devrait toujours être la longueur de la '.poids'), c'est en fait une pondéré réservoir de permutation, mais la méthode fonctionne bien pour l'échantillonnage et la permutation.
Mise à jour: je devrais probablement mentionner que la fonction ci-dessus s'attend à ce que les coefficients de pondération être supérieure à zéro. Sinon, key <- runif(length(.data)) ^ (1 / .weights) ne sera pas commandé correctement.
Juste pour le plaisir, j'ai aussi utilisé le scénario de test dans l'OP, pour comparer les deux fonctions.
set.seed(1)
times_WRS <- ldply(
1:7,
function(i) {
n <- 1024 * (2 ** i)
p <- runif(2 * n)
n_Set <- 1:(2 * n)
data.frame(
n=n,
user=system.time(weighted_Random_Sample(n_Set, p, n), gcFirst=T)['user.self'])
},
.progress='text'
)
sample.int.test <- function(n, p) {
sample.int(2 * n, n, replace=F, prob=p); NULL }
times_sample.int <- ldply(
1:7,
function(i) {
n <- 1024 * (2 ** i)
p <- runif(2 * n)
data.frame(
n=n,
user=system.time(sample.int.test(n, p), gcFirst=T)['user.self'])
},
.progress='text'
)
times_WRS$group <- "WRS"
times_sample.int$group <- "sample.int"
library(ggplot2)
ggplot(rbind(times_WRS, times_sample.int) , aes(x=n, y=user/n, col=group)) + geom_point() + scale_x_log10() + ylab('Time per unit (s)')
Et voici les temps:
times_WRS
# n user
# 1 2048 0.00
# 2 4096 0.01
# 3 8192 0.00
# 4 16384 0.01
# 5 32768 0.03
# 6 65536 0.06
# 7 131072 0.16
times_sample.int
# n user
# 1 2048 0.02
# 2 4096 0.05
# 3 8192 0.14
# 4 16384 0.58
# 5 32768 2.33
# 6 65536 9.23
# 7 131072 37.79
Laissez-moi jeter dans ma propre mise en œuvre d'une rapide approche basée sur le rejet de l'échantillonnage avec remplacement. L'idée est la suivante:
Générer un échantillon avec remplacement qui est "un peu" plus grande que la taille demandée
Jeter les valeurs en double
Si pas assez de valeurs ont été tirées, appelez la même procédure récursive avec ajusté
n,sizeetprobparamètresReconfigurer le retour de l'index de l'index originaux
Quelle est l'ampleur de l'échantillon devons-nous en tirer? En supposant une distribution uniforme, le résultat est le nombre d'essais pour voir x unique des valeurs de N valeurs totales. C'est une différence de deux nombres harmoniques (H_n et H_{n - taille}). Les premiers nombres harmoniques sont totalisées, sinon une approximation en utilisant le logarithme naturel est utilisé. (Ce n'est qu'une approximation de la figure, pas besoin d'être trop précis ici.) Maintenant, pour une distribution non uniforme, le nombre d'éléments attendus pour être tracée ne peut être plus grande, de sorte que nous ne serons pas en tirer trop de samples. En outre, le nombre d'échantillons est limité par deux fois la taille de la population -- je suppose que c'est plus rapide d'avoir quelques appels récursifs de l'échantillonnage jusqu'à O(n ln n) des éléments.
Le code est disponible dans le package R wrswoR dans la sample.int.rej de routine en sample_int_rej.R. Installer avec:
library(devtools)
install_github('wrswoR', 'muelleki')
Il semble que cela fonctionne "assez rapide", mais pas formel d'exécution des tests ont été réalisés encore. Aussi, il est testé sous Ubuntu uniquement. J'apprécie vos commentaires.
