Buongiorno a tutti,

sono uno studente di economia e devo svolgere alcuni esperimenti, utilizzando il software R, per un progetto di statistica economica.

Uno di questi prevede di lanciare un tetraedo regolare 59 volte e di calcolare il punteggio totale; questo esperimento andrebbe ripetuto 100000 volte in modo da stimare con il metodo Monte Carlo la probabilità che il totale sia superiore a 164.

Io ho impostato il problema con il ciclo for in modo da avere come output un vettore y che contenga tutti i risultati delle 100000 prove

> tetraedo<-c(1:4)

> tetraedo

[1] 1 2 3 4

> x<-sample(tetraedo,size=59,replace=TRUE)

> x

 [1] 3 3 1 2 3 3 2 2 1 4 2 4 2 3 4 3 3 4 4 1 2 2 3 1 3 3 1 1 4 4 3 2 3 3 2 4 2 3

[39] 4 3 3 1 1 3 4 1 2 1 2 4 3 3 1 4 3 1 4 4 4

> sum(x)

[1] 156

> for(i in 1:100000)

+ {

+ x<-sample(tetraedo,1)

+ y[i]<-x

+ }

> sum(y)

[1] 249819

> 164<sum(y)

[1] TRUE

Secondo voi può essere corretto risolvere l'esperiemento in questo modo?

Visualizzazioni: 87

Risposte a questa discussione

per risolvere il tuo problema ho creato questo script composto da due funzioni,

- la prima funzione lancioTetraedro() è annidata in getResult() ed esegue il lancio per 59 volte e restituisce la somma del risultato

- la funzione getResult() conta il numero di occorrenze in cui la variabile esperimento è maggiore di 164 e le divide per 100000 (con sapply ho eseguito 100000 volte il lancio del tetratedro memorizzando il risultato sulle celle di un vettore lungo appunto 100000)

getResult <- function(y=rep(NA,100000)){

      lancioTetraedro <- function(y=y){
                                    tetraedro<-c(1:4)
                                    x<-sample(tetraedro,size=59,replace=TRUE)
                                    out<- sum(x)
                                    return(out)
       }

     esperimento <- sapply(y,lancioTetraedro)
     result <- length(esperimento[esperimento>164])/length(esperimento)
     return(result)
}

Se modifichi per esempio il valore di default di y che passi alla funzione può modificare il numero di ripetizioni a tuo piacimento (es. getResult(y=rep(NA,10^7)  hai un test fatto con 10 milioni di lanci)

Non sono un esperto, probabilmente in R ci sono già funzioni predefinite che lo fanno meglio.

Ottima soluzione!

Solo per il gusto di imparare si può efficientare un po' (la soluzione è ottima per il problema proposto, che credo non necessiti per forza tempi rapidi)

per esempio: 

  1. tetraedro<-c(1:4) può essere definito una sola volta all'esterno di lancioTetraedro (è vero che si perde l'integrità funzionale, ma nn siamo puristi). 
  2. y=rep(NA, 100000) sarebbe meglio sostituirlo con y = vector(length = 100000)
  3. Possiamo costruire out come vettore di logical invece che che numeric, inserendo il confronto (out<- sum(x) > 164) all'interno della funzione
  4. Su un vettore di logical possiamo usare sum() invece che length() (è più il subset esperimento[...] che rallenta)
  5. sample! Debuggando (con il package profvis) mi sono accorto che la maggior parte del tempo è usata dalla funzione sample. Ho quindi scoperto che esiste una funzione sample.int (che ovviamente funziona solo si voglia estrarre numeri interi), che è più rapida:

benchmark(sample(1:4, 59, T), sample.int(4, 59, T), replications = 1000000, columns = c('test', 'elapsed', 'relative'))
                         test   elapsed   relative
1 sample(1:4, 59, T)       7.19        1.546
2 sample.int(4, 59, T)     4.65        1.000

Sommando il tutto, questa è la versione più efficiente che son riuscito ad ottenere, ed è quasi due volte e mezza più veloce dell'originale:

sampleint_func <- function(n = 4, lanci = 59, rep = 100000, soglia = 164){

        esperimento <- sapply(1:rep, function(i) {

            risultatoLanci <- sample.int(n, lanci, T)
            sum(risultatoLanci) > soglia
       })

       sum(esperimento) / rep
}

Di seguito codice e risultati del benchmark

original_func <- function(y=rep(NA,100000)){
    set.seed(1234)

    lancioTetraedro <- function(y=y){
        tetraedro<-c(1:4)
        x<-sample(tetraedro,size=59,replace=TRUE)
        out<- sum(x)
        return(out)
    }

    esperimento <- sapply(y,lancioTetraedro)
    result <- length(esperimento[esperimento>164])/length(esperimento)
    return(result)
}


eff_func <- function(y=rep(NA,100000)){
    set.seed(1234)
    tetraedro<-c(1:4)
    lancioTetraedro <- function(y=y){

        x<-sample(tetraedro,size=59,replace=TRUE)
        out<- sum(x) > 164
        return(out)
    }

    esperimento <- sapply(y,lancioTetraedro)
    result <- length(esperimento[esperimento])/length(esperimento)
    return(result)
}


sapply_func <- function(){
    set.seed(1234)
    tetraedo<-c(1:4)

    n <- 100000

    esperimento <- sapply(1:n, function(i) {
        lanci <- sample(tetraedo, 59, T)
        sum(lanci) > 164
    })

    sum(esperimento) / n
}


for_cycle_func <- function() {
    set.seed(1234)
    tetraedo<-c(1:4)

    n <- 100000

    y <- vector(length = n)

    for(i in 1:n) {
        lanci <- sample(tetraedo, 59, T)
        y[i] <- sum(lanci) > 164
    }

    sum(y) / n
}


sample_int_func <- function(n = 4, lanci = 59, rep = 100000, soglia = 164){

    esperimento <- sapply(1:rep, function(i) {

        risultatoLanci <- sample.int(n, lanci, T)
        sum(risultatoLanci) > soglia
    })

    sum(esperimento) / rep
}


tests = list(
    original = quote(original_func()),
    eff = quote(eff_func()),
    sapply = quote(sapply_func()),
    for_cycle = quote(for_cycle_func()),
    sample_int = quote(sample_int_func())
)


library(rbenchmark)
do.call(benchmark,
    c(tests, list(replications=10,
      columns=c('test', 'elapsed', 'relative'),
      order='elapsed')))
#>                 test    elapsed    relative
#> 5     sample_int         3.76     1.000
#> 3           sapply         5.72     1.521
#> 4       for_cycle         5.94     1.580
#> 2                eff          7.12     1.894
#> 1          original         8.73     2.322

Grazie Giacomo!

Da un problema apparentemente semplice si possono sempre imparare nuove cose.

Mi ero accorto anch'io che a livello di benchmark la mia soluzione era poco efficiente..

R lo uso da poco, in passato programmavo un po' in java.

A distanza di anni, non facendolo per mestiere, sto perdendo pezzi.

RSS

Social

 

Gruppi

© 2017   Creato da Duccio Schiavon.   Tecnologia

Badge  |  Segnala un problema  |  Politica sulla privacy  |  Termini del servizio