simple - muestreo con r studio

Muestreo ponderado más rápido sin reemplazo

Question

Esta pregunta llevó a un nuevo paquete R: wrswoR

El muestreo predeterminado de R sin reemplazo usando sample.int parece requerir un tiempo de ejecución cuadrático, por ejemplo, cuando se usan pesos extraídos de una distribución uniforme. Esto es lento para muestras de gran tamaño. ¿Alguien sabe una implementación más rápida que se pueda utilizar desde R ? Dos opciones son "muestreo de rechazo con reemplazo" (consulte esta pregunta en stats.sx) y el algoritmo de Wong y Easton (1980) (con una implementación de Python en una respuesta de StackOverflow ).

Gracias a Ben Bolker por sugerir la función C que se llama internamente cuando se llama a sample.int con replace=F y pesos no uniformes: ProbSampleNoReplace . De hecho, el código muestra dos anidados for bucles (línea 420 ff de random.c ).

Aquí está el código para analizar empíricamente el tiempo de ejecución:

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

EDITAR : Gracias a Arun por señalar que el muestreo no ponderado no parece tener esta penalización de rendimiento.

Answer 1

Decidí profundizar en algunos de los comentarios y encontré que el artículo de Efraimidis y Spirakis era fascinante (gracias a @Hemmo por encontrar la referencia). La idea general en el documento es esta: cree una clave generando un número aleatorio uniforme y elevándolo a la potencia de uno sobre el peso de cada artículo. Luego, simplemente toma los valores clave más altos como su muestra. ¡Esto funciona brillantemente!

weighted_Random_Sample <- function(
.data,
.weights,
.n
){
key <- runif(length(.data)) ^ (1 / .weights)
return(.data[order(key, decreasing=TRUE)][1:.n])
}

Si configura ''.n'' para que sea la longitud de ''.data'' (que siempre debería ser la longitud de ''.weights''), esto es en realidad una permutación ponderada del reservorio, pero el método funciona bien tanto para el muestreo como para la permutación.

Actualización : probablemente debería mencionar que la función anterior espera que los pesos sean mayores que cero. De lo contrario, la key <- runif(length(.data)) ^ (1 / .weights) no se ordenará correctamente.

Sólo por diversión, también usé el escenario de prueba en el OP para comparar ambas funciones.

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)'')

Y aquí están los tiempos:

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

Answer 2

Permítame incluir mi propia implementación de un enfoque más rápido basado en el muestreo de rechazo con reemplazo . La idea es esta:

• Genere una muestra con reemplazo que sea "algo" más grande que el tamaño solicitado

• Deseche los valores duplicados

• Si no se han dibujado suficientes valores, llame al mismo procedimiento de forma recursiva con los parámetros n , size y prob ajustados

• Reasignar los índices devueltos a los índices originales

¿Qué tamaño de muestra necesitamos dibujar? Suponiendo una distribución uniforme, el resultado es el número esperado de ensayos para ver x valores únicos de N valores totales . Es una diferencia de dos números armónicos (H_n y H_ {n - tamaño}). Los primeros números de armónicos se tabulan, de lo contrario se usa una aproximación utilizando el logaritmo natural. (Esto es solo una cifra aproximada, no es necesario que sea demasiado preciso aquí). Ahora, para una distribución no uniforme, el número esperado de elementos a dibujar solo puede ser mayor, por lo que no extraeremos demasiadas muestras. Además, el número de muestras extraídas está limitado por el doble del tamaño de la población. Supongo que es más rápido tener algunas llamadas recursivas que muestrear hasta O (n ln n) elementos.

El código está disponible en el paquete de R wrswoR en la rutina sample_int_rej.R en sample_int_rej.R . Instalar con:

library(devtools)
install_github(''wrswoR'', ''muelleki'')

Parece funcionar "lo suficientemente rápido", sin embargo, aún no se han realizado pruebas formales de tiempo de ejecución. Además, el paquete se prueba solo en Ubuntu. Aprecio sus comentarios.

Answer 3

Actualizar:

Una implementación Rcpp del algoritmo Efraimidis y Spirakis (gracias a @Hemmo, @Dinrem, @krlmlr y @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 

Versión antigua:

Probemos algunos enfoques posibles:

Muestreo de rechazo simple con reemplazo. Esta es una función mucho más simple que sample.int.rej ofrece @krlmlr, es decir, el tamaño de la muestra siempre es igual a n . Como veremos, sigue siendo muy rápido, asumiendo una distribución uniforme para los pesos, pero extremadamente lento en otra situación.

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]
}

El algoritmo de Wong y Easton (1980) . Aquí hay una implementación de esta versión de Python. Es estable y puede que me falte algo, pero es mucho más lento en comparación con otras funciones.

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
}

Implementación Rcpp del algoritmo por Wong y Easton. Posiblemente se puede optimizar aún más ya que esta es mi primera función Rcpp utilizable, pero de todos modos funciona 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")

Ahora algunos resultados:

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

Obviamente, podemos rechazar la función 1980 porque es más lenta que la Base en ambos casos. Rejection simple mete en problemas también cuando hay una probabilidad única de 0.999 en el segundo caso.

Así que queda el Rejection , Rcpp , Reservoir . El último paso es verificar si los valores en sí son correctos. Para estar seguros acerca de ellos, utilizaremos la sample como punto de referencia (también para eliminar la confusión sobre las probabilidades que no tienen que coincidir con p debido al muestreo sin reemplazo).

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 

Note algunas cosas aquí. Por alguna razón, weighted_Random_Sample devuelve valores incorrectos (no lo he buscado en absoluto, pero funciona correctamente suponiendo una distribución uniforme). sample.int.rej es muy lento en el muestreo repetido.

En conclusión, parece que Rcpp es la opción óptima en caso de un muestreo repetido, mientras que sample.int.rej es un poco más rápido y también más fácil de usar.