######################
## 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(seed = 1)
runif(n = 5)

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

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

# 3) controllare graficamente l'uniformita del generatore casuale
# usare 10000 simulazioni e i comandi plot e hist / density
U1 <- runif(n = 10000)
hist(U1)
hist(U1, freq = FALSE, breaks = 50)
abline(h = 1)

plot(density(x = U1))

hist(U1, freq = FALSE, breaks = 50)
lines(density(x = U1), col = "red")

# 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(nsim, rate = 1)
E2 <- function()-log(1 - runif(nsim))
E3 <- function()-log(runif(nsim))

system.time(E1())
system.time(E2())
system.time(E3())


# 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 + ... + Yn^2 dove le Yi sono normali standard
# usare "apply" e "colSums
set.seed(0)
nsim <- 10 ^ 6

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

C2 <- function(df){
  stopifnot(df %% 2 == 0)
  nexp <- df / 2
  U <- matrix(data = runif(nsim * nexp), nrow = nexp, byrow = FALSE) # matrice con nexp righe, nsim colonne, riempita con U(0,1)

  X <- -log(U) #matrice con esponenziali di parametro 1
  X <- 2 * apply(X = X, MARGIN = 2, FUN = sum)
  return(X)
  }

system.time(C1(6))
system.time(C2(6))

system.time(C1(20))
system.time(C2(20))

# si velocizza se usiamo "colSums"?

C3 <- function(df){
  stopifnot(df %% 2 == 0)
  nexp <- df / 2

  X <- 2 * colSums(-log(matrix(data = runif(nsim * nexp), nrow = nexp, byrow = FALSE)))
  return(X)
}

system.time(C1(20))
system.time(C2(20))
system.time(C3(20))

C4 <- function(df){
  X <- colSums(matrix(data = rnorm(nsim * df) ^ 2, nrow = df))
  return(X)
  }

system.time(C1(20))
system.time(C2(20))
system.time(C3(20))
system.time(C4(20))


# se si aumentano i gradi di liberta':


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

TC1

# verificare distribuzione casuale usando "table" e "replicate" (valutazioni ripetute di una data espressione)

table(TC1)

TC1.1 <- replicate(expr = rbinom(n = 500, size = 1, prob = 0.5), n = 10)
dim(TC1.1)
apply(X = TC1.1, MARGIN = 2, FUN = table)

# con "runif"

set.seed(0)
TC2 <- runif(n = 500) < 0.5
head(TC2)

TC2.1 <- 1 * (runif(n = 500) < 0.5)
head(TC2.1)


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

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

sample(x = c("T", "C"), size = 500)

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

nsim <- 5 * 10 ^ 7

set.seed(1234)
system.time(rbinom(n = nsim, size = 1, prob = 0.5))

set.seed(1234)
system.time(1 * (runif(n = nsim) < 0.5))

set.seed(1234)
system.time(sample(x = c("T", "C"), size = nsim, replace = TRUE))

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

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

P1.1 <- sample(x = c(-2, 1, 5), size = 5000, replace = TRUE, prob = c(1 / 4, 1 / 3, 5 / 12))

table(P1.1) / 5000
c(1 / 4, 1 / 3, 5 / 12)

# 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

estrazioni.urna <- function(nomi.colori = c("R", "B"),
                            numero.sfere = rep(x = 1, length.out = length(nomi.colori)), # default: 1 pallina per colore
                            numero.estrazioni = 1,
                            rimpiazzo = TRUE){

  l.n.sfere <- length(numero.sfere)
  if (l.n.sfere == 1) numero.sfere <- rep(x = numero.sfere, length.out = length(nomi.colori)) else if (l.n.sfere != length(nomi.colori)) stop("nomi.colori e numero.sfere devono avere la stessa lunghezza")

  URNA <- rep(nomi.colori, times = numero.sfere)

  estrazione <- sample(x = URNA, size = numero.estrazioni, replace = rimpiazzo)

  return(estrazione)

  }

estrazioni.urna(nomi.colori = c("R", "B", "N"), numero.sfere = 10, numero.estrazioni = 15)

estrazioni.urna(nomi.colori = c("R", "B", "N"), numero.sfere = c(50, 10, 1), numero.estrazioni = 15)

estrazioni.urna(nomi.colori = c("R", "B", "N"), numero.sfere = c(5, 10, 1), numero.estrazioni = 20)

estrazioni.urna(nomi.colori = c("R", "B", "N"), numero.sfere = c(5, 10, 1), numero.estrazioni = 20, rimpiazzo = FALSE)

debug(estrazioni.urna)

estrazioni.urna(nomi.colori = c("R", "B", "N"), numero.sfere = c(5, 10, 1), numero.estrazioni = 20, rimpiazzo = FALSE)


# 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 una sola ruota

simula.ruota()

simula.lotto <- function(){
  lotto <- t(replicate(expr = simula.ruota(), n = 11))
  rownames(lotto) <- ruote
  colnames(lotto) <- c("1o", "2o", "3o", "4o", "5o")
  return(lotto)
  }

set.seed(1234)
simula.lotto()

replicate(expr = simula.lotto(), n = 3)

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

ranghi <- c("A", 2:10, "J", "Q", "K")
semi <- c("C", "Q", "F", "P")

paste(ranghi[11], semi[3], sep = "")

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

mano.poker <- function(n.giocatori){

  mano <- sample(x = mazzo, size = 5 * n.giocatori, replace = FALSE)

  mano <- matrix(data = mano, nrow = n.giocatori, ncol = 5)

  rownames(mano) <- paste("G", 1:n.giocatori, sep = ".")
  return(mano)
  }

set.seed(1234)
mano.poker(n.giocatori = 4)