######################
## Simulazione in R ##
######################

# le distribuzioni piu' comuni sono contenute nella distribuzione base di R

# esempio: distribuzione normale
# dnorm: densita'
# pnorm: funzione di ripartizione
# qnorm: funzione quantile (inversa generalizzata)
# rnorm: simulazione

# Beta: rbeta
# Binomiale: rbinom
# Cauchy: rcauchy
# Chi-quadrato: rchisq
# Esponenziale: rexp
# F: rf
# Gamma: rgamma
# Geometrica: rgeom
# Ipergeometrica: rhyper
# Lognormale: rlnorm
# Logistica: rlogis
# Normale: rnorm
# Poisson: rpois
# Student: rt
# Uniforme: runif
# Weibull: rweibull

# molte altre distribuzioni sono contenute in pacchetti aggiuntivi
# e.g.: distribuzione di pareto nel pacchetto 'actuar'

# runif fornisce il generatore casuale in (0,1)
# il seme viene controllato tramite il comando set.seed

# 1) imposta il seme a 1 e genera 5 uniformi U(0,1)
set.seed(1)
runif(n = 5, min = 0, max = 1)


# 2) genera 5 uniformi U(6,10)
set.seed(1)
runif(n = 5, min = 6, max = 10)

set.seed(1)
4 * runif(n = 5) + 6


# 3) controllare graficamente l'uniformita del generatore casuale
# usare 10000 simulazioni e i comandi plot e hist / density
set.seed(1)
U1 <- runif(n = 10000)
plot(U1)

hist(x = U1, freq = F, breaks = 20)
plot(density(x = U1))


# 4) testare l'indipendenza delle uniformi simulate al punto precedente con il test Box-Pierce (o Ljung-Box)
# H0: i dati vengono da variabili indipendenti

Box.test(x = U1)


# 5) set.seed si puo' usare per impostare il generatore casuale oppure usare RNGkind
# NON CAMBIARE GENERATORE CASUALE SE NON SI HANNO RAGIONI SPECIFICHE PER FARLO!!!!

RNGkind(kind = NULL, normal.kind = NULL)

# kind permette di scegliere il generatore casuale
# in R sono implementati i seguenti:
# "Wichmann-Hill": periodo ~ 7*10^12
# "Marsaglia-Multicarry": periodo ~ 2^60
# "Super-Duper": periodo ~ 4.6*10^18
# "Mersenne-Twister" (default): periodo ~ 2^19937-1 ********
# "Knuth-TAOCP-2002": periodo ~ 2^129.
# "Knuth-TAOCP"
# "L'Ecuyer-CMRG"
# "user-supplied"

# 6) NON USARE METODI COSTRUITI "ALLA MANO" SE METODI DIRETTI SONO GIA' PRESENTI IN R
# => CONTROPRODUCENTE
# => INEFFICIENTE
# => PERICOLOSO

# esempio: esponenziale usando rexp oppure con il metodo dell'inversa => poca differenza

nsim <- 10 ^ 6 # un milione di simulazioni
set.seed(1)

E1 <- function()rexp(n = nsim, rate = 1)
system.time(E1())

E2 <- function()-log(runif(n = nsim))
system.time(E2())


# 7) altro esempio: chi-quadrato con 2n gradi di liberta' (n intero) puo' essere ottenuta come 2(X1 + X2 + ... + Xn) dove Xi sono esponenziali di parametro 1 indipendenti oppure come Y1^2 + ... + Y2n^2 dove le Yi sono normali standard
# usare "apply" e "colSums"
set.seed(0)
nsim <- 10 ^ 7

C1 <- function(df)rchisq(n = nsim, df = df)

C2 <- function(df){
  stopifnot(df %% 2 == 0)
 
  nexp <- df / 2
  U <- matrix(data = runif(n = nsim * nexp), nrow = nexp) # matrice di U(0,1)
  X <- -log(U) # matrice di esponenziali unitarie
  Y <- 2 * apply(X = X, MARGIN = 2, FUN = sum)
}
  
system.time(C1(6))
system.time(C2(6))

# si velocizza se usiamo "colSums"?

C2.1 <- function(df){
  stopifnot(df %% 2 == 0)
  nexp <- df / 2
  2 * colSums(-log(matrix(data = runif(n = nsim * nexp), nrow = nexp)))
}

system.time(C2.1(6))


C3 <- function(df)colSums(matrix(data = rnorm(n = nsim * df), nrow = df) ^ 2)

system.time(C3(6))

C2.2 <- function(df){
  stopifnot(df %% 2 == 0)
  nexp <- df / 2
  
  X <- -log(matrix(data = runif(n = nsim * nexp), nrow = nexp))
  Y <- vector(mode = "numeric", length = nsim)
  for (i in 1 : nsim){
    Y[i] <- 2 * sum(X[, i])
  }
  return(Y)
}

system.time(C2.2(6))
system.time(C1(6))
system.time(C2.1(6))
system.time(C3(6))

?apply
?sapply
?lapply
...

# se si aumentano i gradi di liberta':

nsim <- 10 ^ 6

system.time(C2.2(30))
system.time(C1(30))
system.time(C2.1(30))
system.time(C3(30))

# 8) simulare 500 lanci di moneta equa, T = 0, C = 1 
# con "rbinom", "runif" e "sample"

set.seed(1)
TC1 <- rbinom(n = 500, size = 1, prob = 0.5)

# verificare distribuzione casuale usando "table" e "replicate" (valutazioni ripetute di una data espressione)
TC1.1 <- replicate(n = 10, expr = rbinom(n = 500, size = 1, prob = 0.5))
dim(TC1.1)
apply(X = TC1.1, MARGIN = 2, FUN = table)

# con "runif"

TC2 <- 1 * (runif(500) <= 0.5)

# con "sample"
# "sample" permette di simulare (campionare) da una distribuzione finita, con o senza rimpiazzo

TC3 <- sample(x = c(0, 1), size = 500, replace = TRUE)
table(TC3)
TC3 <- sample(x = c("T", "C"), size = 500, replace = TRUE)

# qual'e' il piu' efficiente (veloce)?
# (con 5000000 simulazioni)

nsim <- 5 * 10 ^ 6

system.time(rbinom(n = nsim, size = 1, prob = 0.5))
system.time(1 * (runif(n = nsim) < 0.5))
system.time(sample(x = c("T", "C"), replace = TRUE, size = nsim))

# 9) uso di sample: simulare 50 estrazioni da una distribuzione discreta -2, 1, 5 con prob 1/4, 1/3, 5/12

X <- sample(x = c(-2, 1, 5), size = 50000, replace = TRUE, prob = c(1 / 4, 1 / 3, 5 / 12))
X
table(X) / 50000



# 10) scrivere una funzione per simulare estrazioni di sfere colorate da un'urna con e senza rimpiazzo
# la funzione deve prendere in input i nomi dei colori, il numero di sfere per ogni colore, il numero di estrazioni e se avvengono con o senza rimpiazzo


estrai.urna <- function(nomi.colori = c("R", "B"), numero.sfere = rep(x = 1, length.out = length(nomi.colori)), numero.estrazioni = 1, rimpiazzo = TRUE){
  
  l.n.sfere <- length(numero.sfere)
  if (l.n.sfere == 1) numero.sfere <- rep(x = n.sfere, length.out = length(nomi.colori))
  
  URNA <- rep(nomi.colori, times = numero.sfere)
  URNA <- paste(URNA, 1 : length(URNA))
  
  estrazioni <- sample(x = URNA, size = numero.estrazioni, replace = rimpiazzo)
  
  return(estrazioni)
  }

estrai.urna()

estrai.urna(nomi.colori = c("V", "R", "B"), numero.sfere = c(3, 5, 10), numero.estrazioni = 10, rimpiazzo = FALSE)

debug(estrai.urna)
undebug(estrai.urna)

# 4 colori: 5 Gialle, 10 Verdi, 5 Rosse, 8 Bianche
estrai.urna(nomi.colori = c("G", "V", "R", "B"), numero.sfere = c(5, 10, 5, 8), numero.estrazioni = 20, rimpiazzo = FALSE)

estrai.urna(nomi.colori = c("G", "V", "R", "B"), numero.sfere = c(5, 10, 5, 8), numero.estrazioni = 30, rimpiazzo = FALSE) # errore, troppe estrazioni

# 11) simulare estrazione del lotto su 10+1 ruote

ruote <- c("BA", "CA", "FI", "GE", "MI", "NA", "PA", "Roma", "TO", "VE", "Nazionale")

numeri <- 1 : 90

simula.ruota <- function()sample(x = numeri, size = 5, replace = FALSE)
simula.ruota()

simula.lotto <- function(){
  lotto <- t(replicate(n = length(ruote), expr = simula.ruota()))
  rownames(lotto) <- ruote
  colnames(lotto) <- paste(1 : 5, "o", sep = "")
  return(lotto)
  }

simula.lotto()



# 12) simulare la prima mano a poker con 4 giocatori

ranghi <- c("A", as.character(2 : 10), "J", "Q", "K")

semi <- c("C", "Q", "F", "P")

# "CA" "Q6"... "PK" ... 52 elementi

mazzo <- c(outer(X = ranghi, Y = semi, FUN = paste, sep = ""))

prima.mano.poker <- function(n.giocatori){
  stopifnot(n.giocatori * 5 < 52)

  m <- matrix(data = sample(x = mazzo, size = n.giocatori * 5, replace = FALSE), nrow = n.giocatori, ncol = 5)
  rownames(m) <- paste("Giocatore", 1 : n.giocatori, sep = " ")
  return(m)
  }

prima.mano.poker(n.giocatori = 3)

prima.mano.poker(n.giocatori = 5)

prima.mano.poker(n.giocatori = 12)