Los bucles grandes se cuelgan en R?

Question

Supongamos que quiero realizar una simulación utilizando la siguiente function:

fn1 <- function(N) { 
    res <- c() 
    for (i in 1:N) { 
    x <- rnorm(2) 
    res <- c(res, x[2]-x[1]) 
    } 
    res 
} 

Por muy grande N, computación parece bloquearse. ¿Hay mejores formas de hacer esto?

(Inspirado por: https://stat.ethz.ch/pipermail/r-help/2008-February/155591.html)

Answer 1

Para bucles en R son notoriamente lentos, pero aquí hay otro problema. Es mucho más rápido preasignar el vector de resultados, res, en lugar de anexar res en cada iteración.

A continuación podemos comparar la velocidad de la versión anterior con una versión que simplemente comienza con un vector, res, de longitud N y cambia el elemento ith durante el ciclo.

fn1 <- function(N) { 
    res <- c() 
    for (i in 1:N) { 
    x <- rnorm(2) 
    res <- c(res,x[2]-x[1]) 
    } 
    res 
} 
fn2 <- function(N) { 
    res <- rep(0,N) 
    for (i in 1:N) { 
    x <- rnorm(2) 
    res[i] <- x[2]-x[1] 
    } 
    res 
} 
> N <- 50000 
> system.time(res1 <- fn1(N)) 
    user system elapsed 
    6.568 0.256 6.826 
> system.time(res2 <- fn2(N)) 
    user system elapsed 
    0.452 0.004 0.496 

También, como Sharpie points out, podemos hacer de este un poco más rápido mediante el uso de funciones como R apply (o sus parientes, y sapplylapply).

fn3 <- function(N) { 
    sapply(1:N, function(i){ x <- rnorm(2); return(x[2] - x[1]) }) 
} 
> system.time(res3 <- fn3(N)) 
    user system elapsed 
    0.397 0.004 0.397 

Answer 2

La eficiencia de bucles se puede aumentar enormemente en R a través del uso de las funciones de aplicación que procesan esencialmente vectores enteros de datos a la vez en lugar de bucle a través de ellos. Para el bucle mostrado anteriormente, hay dos operaciones básicas que ocurren durante cada iteración:

# A vector of two random numbers is generated 
x <- rnorm(2) 

# The difference between those numbers is calculated 
x[2] - x[1] 

En este caso la función apropiada sería sapply(). sapply() opera en una lista de objetos, tales como el vector generado por la sentencia de bucle 1:N y devuelve un vector de los resultados:

sapply(1:N, function(i){ x <- rnorm(2); return(x[2] - x[1]) }) 

Tenga en cuenta que el valor del índice i está disponible durante la llamada de función y sucesivamente toma los valores entre 1 y N, sin embargo, no es necesario en este caso.

Adquirir el hábito de reconocer dónde se puede usar apply en for es una habilidad muy valiosa: muchas bibliotecas R para cómputo paralelo proporcionan paralelización plug-and-play a través de las funciones apply. El uso de apply a menudo permite el acceso a aumentos de rendimiento significativos en sistemas multinúcleo con cero refactorización de código.

Answer 3

A veces, el bucle no es necesario. Desde rnorm da muestra iid (teóricamente), va a lograr el mismo resultado (muestreo X-Y donde X e Y son N (0,1)) haciendo:

res <- rnorm(N)-rnorm(N) 

Answer 4

Ampliando mi comentario a la respuesta de chris_dubois, aquí está alguna información de temporización:

> system.time(res <- rnorm(50000) - rnorm(50000)) 
user system elapsed 
0.06 0.00 0.06 

Contraste esto con na3 de esa misma respuesta:

> system.time(res3 <- fn3(50000)) 
user system elapsed 
1.33 0.01 1.36 

La primera cosa a notar es que mi vuelta la parte superior es más lenta que la máquina de chris_dubois.:)

El segundo punto, y más importante, es que el enfoque vectorial, bastante aplicable aquí, es un orden de magnitud más rápido. (También señalado por Richie Cotton en un comentario a la misma respuesta).

Esto me lleva al último punto: es un mito que apply y sus amigos son mucho más rápido que for bucles en R. Están en el mismo orden en la mayoría de las mediciones que he visto. Porque solo son for bucles detrás de escena. Ver también este post:

http://yusung.blogspot.com/2008/04/speed-issue-in-r-computing-apply-vs.html

Según el profesor Brian Ripley "apply() es sólo un envoltorio para un bucle." ¡La única ventaja de usar apply() es que hace que su código sea más ordenado!

Exactamente. Debe usar apply si es más expresivo, especialmente si está programando en un estilo funcional. No porque sea más rápido.

Answer 5

Tal vez el sustituto más eficiente para su función sería simplemente:

fn <- function(n) rnorm(N,0,sqrt(2)) 

que es dos veces tan rápido como teniendo diferencia de variables normales IID. En términos más generales, si su objetivo es ejecutar simulaciones simples, la preasignación de vector/matriz y las llamadas a funciones nativas aceleran el proceso en gran medida.

Si desea realizar simulaciones de monte-carlo para estimaciones estadísticas (por ejemplo, MCMC), R tiene una serie de paquetes nativos. Para la simulación estocástica general, no conozco los paquetes R, pero es posible que desee probar Simpy (http://simpy.sourceforge.net/), que es excelente.