Grandi anelli appesi in R?

Question

Supponiamo che io voglio eseguire una simulazione utilizzando il seguente function:

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

Per molto grande N, calcolo sembra bloccarsi. Ci sono modi migliori per farlo?

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

Answer 1

Per loop in R sono notoriamente lento, ma qui c'è un altro problema. È molto più veloce preallocare i risultati vettore, res, piuttosto aggiungere a res a ogni iterazione.

Di seguito è possibile confrontare la velocità della versione precedente con una versione che inizia semplicemente con un vettore, res, di lunghezza N e modifica l'elemento ith durante il 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 

Inoltre, come Sharpie points out, siamo in grado di rendere questo leggermente più veloce utilizzando le funzioni di ricerca, come apply (o dei suoi parenti, sapply e lapply).

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

L'efficienza di spire può essere aumentata enormemente in R attraverso l'uso delle funzioni che elaborano applicare essenzialmente interi vettori di dati contemporaneamente anziché scorrendo loro. Per il ciclo mostrato sopra, ci sono due operazioni fondamentali che occorrono durante ogni iterazione:

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

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

In questo caso la funzione appropriata sarebbe sapply(). sapply() opera su un elenco di oggetti, come ad esempio il vettore generato dal l'istruzione ciclica 1:N e restituisce un vettore di risultati:

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

noti che il valore di indice i è disponibile durante la chiamata di funzione e successivamente assume i valori tra 1 e N, tuttavia non è necessario in questo caso.

Ottenere l'abitudine di riconoscere dove apply può essere utilizzato più for è un grande valore skill- molte librerie R per il calcolo parallelo forniscono plug-and-play parallelizzazione attraverso apply funzioni. L'utilizzo di apply può spesso consentire l'accesso a significativi aumenti delle prestazioni su sistemi multicore con zero refactoring del codice.

Answer 3

A volte il ciclo non è necessario. Dal momento che RNorm dà campione iid (in teoria), si otterrà lo stesso risultato (campionamento X-Y dove X e Y sono N (0,1)) facendo:

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

Answer 4

Ampliando il mio commento alla risposta di chris_dubois, ecco alcune informazioni di temporizzazione:

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

Contrasto questo con FN3 da quella stessa risposta:

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

La prima cosa da notare è che il mio giro la parte superiore è più lenta della macchina di chris_dubois.:)

Il secondo, e più importante, punto è che l'approccio vettoriale, qui applicabile, è un ordine di grandezza più veloce. (Anche sottolineato da Richie Cotton in un commento alla stessa risposta).

Questo mi porta al punto finale: si tratta di un mito che apply ed i suoi amici sono molto più veloce di for loop sono sulla stesso ordine nella maggior parte delle misure che ho visto in R.. Perché sono solo for loop dietro le quinte. Vedi anche questo post:

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

Secondo il professor Brian Ripley, "apply() è solo un wrapper per un ciclo." L'unico vantaggio dell'utilizzo di apply() è che rende il tuo codice più ordinato!

Esattamente

. Dovresti usare apply se è più espressivo, specialmente se stai programmando in uno stile funzionale. Non perché è più veloce.

Answer 5

Forse la sostituzione più efficiente per la vostra funzione sarebbe semplicemente:

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

che è due volte più veloce come prendere differenza di variates normali iid. Più in generale, tuttavia, se il tuo obiettivo è quello di eseguire semplici simulazioni, la preallocazione di vettori/array e le chiamate alle funzioni native accelerano notevolmente il processo.

Se si desidera eseguire simulazioni di monte-carlo per stime statistiche (ad es. MCMC), R ha un numero di pacchetti nativi. Per la simulazione stocastica generale, non sono a conoscenza dei pacchetti R ma si potrebbe voler provare Simpy (http://simpy.sourceforge.net/), che è eccellente.