Domanda

Ogni volta che eseguo simulazioni di monte carlo su larga scala in S-Plus, finisco sempre per far crescere la barba mentre aspetto che si completi.

Quali sono i trucchi migliori per eseguire simulazioni Monte Carlo in R? Qualche buon esempio di esecuzione di processi in modo distribuito?

È stato utile?

Soluzione

  • L'uso di più core / macchine dovrebbe essere semplice se stai usando repliche indipendenti parallele , ma fai attenzione alle carenze comuni dei generatori di numeri casuali (ad es. se usi l'ora corrente come seed , generando molti processi con un RNG per ciascuno potrebbe produrre numeri casuali correlati, il che porta a risultati non validi - vedi ad esempio questo documento )

  • Potresti voler utilizzare riduzione della varianza per ridurre il numero di repliche richieste , cioè per ridurre le dimensioni del campione richiesto. Tecniche di riduzione della varianza più avanzate sono disponibili in molti libri di testo, ad es. in questo .

Altri suggerimenti

Preallocazione dei tuoi vettori!

> nsims <- 10000
> n <- 100
> 
> system.time({
     res <- NULL
     for (i in 1:nsims) {
         res <- c(res,mean(rnorm(n)))
     }
 })
   user  system elapsed 
  0.761   0.015   0.783 
> 
> system.time({
     res <- rep(NA, nsims)
     for (i in 1:nsims) {
         res[i] <- mean(rnorm(n))
     }
 })
   user  system elapsed 
  0.485   0.001   0.488 
> 

Il campionamento latino Hypercube si applica facilmente e ha una grande influenza sui risultati. Fondamentalmente si prende un campione di ipercubo latino da una distribuzione uniforme (ad es. Usando randomLHS () nel pacchetto lhs) e lo si trasforma nella distribuzione desiderata usando ad es. Qnorm (uniformsample).

So che questo thread è davvero vecchio, ma se qualcuno si imbatte in esso e sta cercando un metodo ancora più veloce, penso che i seguenti lavori:

library(data.table)
library(microbenchmark)

nsims <- 10000
n <- 100

# Answer from @Eduardo_Leoni:
preallocate<-function(nsims, n) {
  res <- rep(NA, nsims)
  for (i in 1:nsims) {
    res[i] <- mean(rnorm(n))
  }
  return(res)
}

# Answer using data.table:
datatable<-function(nsims,n) {
  dt <- data.table(i=1:nsims)[,list(res=mean(rnorm(1:n))),by=i]
  return(dt)
}

# Timing benchmark:
microbenchmark(preallocate(nsims,n), datatable(nsims,n), times=100)
#Unit: milliseconds
#                  expr      min       lq   median       uq      max neval
# preallocate(nsims, n) 428.4022 432.3249 434.2910 436.4806 489.2061   100
#   datatable(nsims, n) 238.9006 242.3517 244.1229 246.5998 303.6133   100
Autorizzato sotto: CC-BY-SA insieme a attribuzione
Non affiliato a StackOverflow
scroll top