Duże pętle zawieszone w R?

Question

Załóżmy, że chcemy przeprowadzić symulację przy użyciu następujących function:

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

dla bardzo dużych N, obliczeń wydaje się powiesić. Czy istnieją lepsze sposoby robienia tego?

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

Answer 1

Na pętli w R są notorycznie powolny, ale tu jest inny problem. O wiele szybszy jest wstępny przydział wektora wyników, res, raczej dołączanie do res w każdej iteracji.

Poniżej możemy porównać prędkość powyżej wersji z wersji, która po prostu zaczyna się z wektorem, res, o długości N i-ty element zmienia się podczas pętli.

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 

Ponadto, jak Sharpie points out, możemy zrobić to nieco szybciej za pomocą funkcji R jak apply (lub jego krewnych, sapply i 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

Wydajność pętli może być znacznie zwiększona, w R poprzez zastosowanie funkcji, które mają zastosowanie zasadniczo przetwarzać cały wektorów danych na raz zamiast pętli przez nie. Dla pętli pokazanej powyżej, istnieją dwa podstawowe działania zachodzące w każdej iteracji:

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

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

W tym przypadku odpowiednie byłoby sapply() funkcji. sapply() działa na liście obiektów, takich jak wektora wygenerowanego przez instrukcję pętli 1:N i zwraca wektor wyników:

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

pamiętać, że wartość indeksu i jest dostępny podczas wywołania funkcji i sukcesywnie zaczyna się na wartościach między 1 i N, jednak nie jest to konieczne w tym przypadku.

Pierwsze w zwyczaju uznając gdzie apply mogą być wykorzystane przez for jest bardzo cennym skill- wiele bibliotek R dla obliczeń równoległych zapewnić plug-and-play Zrównoleglanie przez apply funkcji. Korzystanie apply często może zezwolić na dostęp do znacznego wzrostu wydajności w systemach wielordzeniowych z zerowy refaktoryzacji kodu.

Answer 3

Czasami pętla nie jest potrzebne. Od rnorm daje próbkę IID (teoretycznie), można osiągnąć ten sam rezultat (próbkowanie X-Y gdzie X i Y są N (0,1)), wykonując:

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

Answer 4

Rozszerzenie na moim komentarzu do odpowiedzi chris_dubois jest, tutaj niektóre informacje taktowania:

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

Kontrast to z FN3 z tą samą odpowiedź:

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

pierwszą rzeczą, aby zauważyć, że kolanach góra jest wolniejsza od maszyny chris_dubois.:)

Drugim i ważniejszym punktem jest to, że podejście wektorowe, mające tu zastosowanie, jest o rząd wielkości szybsze. (Zwrócił także uwagę Richie Cotton w komentarzu do tej samej odpowiedzi).

To doprowadza mnie do ostatniego punktu: to mit że apply i jego przyjaciele są znacznie szybsze niż for pętli w R. Są w tej samej kolejności, w większości pomiarów widziałem. Ponieważ są tylko for pętle za kulisami. Zobacz także tego posta:

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

Według profesora Briana Ripley "apply() jest tylko nakładką na pętli." Jedyną korzyścią wynikającą z zastosowania apply() jest to, że czyni twój kod bardziej uporządkowanym!

Dokładnie. Powinieneś użyć apply, jeśli jest bardziej ekspresowy , zwłaszcza jeśli programujesz w stylu funkcjonalnym. Nie dlatego, że jest szybszy.

Answer 5

Może najbardziej efektywny zamiennik dla funkcji byłoby po prostu:

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

który jest dwa razy szybciej niż przy różnicy IID normalnych zmiennymi. Mówiąc bardziej ogólnie, jeśli Twoim celem jest uruchomienie prostych symulacji, wstępna alokacja wektorów/tablic i wywoływanie funkcji natywnych znacznie przyspiesza ten proces.

Jeśli chcesz wykonywać symulacje monte-carlo do estymacji statystycznych (np. MCMC), R ma wiele rodzimych pakietów. Do ogólnej symulacji stochastycznej nie jestem świadomy pakietów R, ale możesz spróbować Simpy (http://simpy.sourceforge.net/), która jest doskonała.