# generatore random con valori che seguono la distribuzione del # del prodotto dei valori delle facce di due dadi regolari # matrice delle 36 possibilità mp <- outer(1:6, 1:6, "*") # tabella del numero di occorrenze dei prodotti tmp <- table(as.vector(mp)) # tabella delle probabilità (casi favorevoli su casi possibili) round( pmp <- tmp / sum(tmp[]), 3) plot(pmp, xlab='prod', ylab='P(prod)', main='Lancio di due dadi') # estraiamo i possibili valori dei prodotti # nella tabella erano dei "nomi", che vanno opportunamente # convertiti (eseguire un passaggio alla volta per capire meglio) p <- as.integer(as.vector(dimnames(tm1))[[1]]) # pdf fp <- as.vector(pmp[]) # cumulativa Fx <- cumsum(fx) rFx <- runif(1) # si ricerca l'indice mediante ricerca binaria for (i in 1:length(Fx)) { if (rFx < Fx[i]) break } # il numero random sarà quello in corrispondenza di quell'indice xr <- p[i] # ci costruiamo una funzione ad hoc (lenta, ma che funziona) find.index <- function(rFx, Fx) { for (i in 1:length(Fx)) { if (rFx < Fx[i]) break } return(i) } # => si lascia come ESERCIZIO fare una function con una ricerca # più inteligente # finalmente possiamo generare tanti numeri casuali nr <- 10000 xr <- rep(0, nr) for (i in 1:nr) { xr[i] <- p[ find.index(runif(1), Fx) ] } # plottiamo il risultato plot(table(xr)) #----------------------------------------------------------------------- # OVVIAMENTE questo è un esempio didattico nel casi si abbia # la tabella dei valori con le rispettive probabilità (tmp dell'esempio) # conoscendo il fenomeno aleatorio sottostante, la simulazione è MOLTO # più SEMPLICE e FACILE nr <- 10000 xr <- rep(0, nr) for (i in 1:nr) { d1 <- sample(1:6)[1] d2 <- sample(1:6)[1] xr[i] <- d1*d2 } # plottiamo il risultato plot(table(xr))