# 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))