# qualche banale simulazione:
 
# stimiamo pi greco:   -----------------------------------------------
 
n <- 1000000
x <- runif(n)
y <- runif(x)        # usa la lunghezza di x  
 
hit <- length(x[x^2 + y^2 <= 1])
 
pi  <- hit/n * 4     # attenzione: va specificato  
                     # meglio il tipo di inferenza che si sta facendo 
                     # basata su P(hit) = (pi/4) / 1
 
 
# Teorema del limite centrale da funzione 'strana' -------------------
# costante in [0 , 1/4] e  [3/4, 1], 
 
# diveri metodi per estrarre secondo tale funzione
 
# 0) metodo logicamente piu' pulito, ma lento (tante estrazioni)
myr <- function(n) {
    r <- numeric()
    i <- 0
    while (length(r) < n) {
      i <- i+1
      if (runif(1) > 0.5) {
        r[i] <- runif(1, 0.75, 1)
      }
      else {
        r[i] <- runif(1, 0., 0.25)
      } 
    }
    return(r) 
}
 
# 1) metodo veloce usiando recursione
# buttiamo il centro ed estraiamo i numeri mancanti
myr1 <- function(n) {
    r <- runif(n)                 
    r <- r[r <= 0.25 | r >=0.75] 
    if (length(r) < n) r <- c(r, myr(n-length(r)) ) 
    return(r) 
}
 
 
# 2) recuperiamo i numeri buttati ridistrubuendoli a destra e sinistra
myr2 <- function(n) {
    r <- runif(n)                 
    r[r < 0.25 & r <=0.5] <-  r[r <= 0.25 & r <=0.5] - 0.25
    r[r < 0.25 & r <0.55] <-  r[r <= 0.25 & r <=0.5] + 0.25
    return(r) 
}
# comunque, la myr1 e' quella piu' veloce
 
# 3) decidiamo con una multinomiale numeri generare nella parte bassa
# e quanti nella parte alta, quindi estraiamo secondo uniformi    
myr3 <- function(n) {
    ab <- as.vector(rmultinom(1, n, c(0.5,0.5)))
    r <- c( runif(ab[1],0,0.25),    runif(ab[2],0.75,1) )
    return( sample(r) )   # sample() 'mischia' a caso r
}
 
# myr3 e' decisamente il metodo piu' veloce!
 
# ora facciamo la funzione che fa anche le somme:
sumyr <- function(n, m=10000) {
   rs <- rep(0,m)
   for (i in 1:n) {
      rs <- rs + myr3(m)
   }
   return(rs)  
}
 
# e questo e' lo script che la chiama, ad esempio
x2 <- sumyr(2)        # oppure aumentando le estrazioni
x2 <- sumyr(2,100000)
hist(x2,ncl=100,prob=T)
 
# riproduciamo la figura della dispensa:
layout(matrix(1:5, 5, 1)) # splittiamo lo schermo
x1 <- sumyr(1) ;  hist(x1,ncl=100,prob=TRUE)
x2 <- sumyr(2) ;  hist(x2,ncl=100,prob=TRUE)
x3 <- sumyr(3) ;  hist(x3,ncl=100,prob=TRUE)
x10 <- sumyr(10) ;  hist(x10,ncl=100,prob=TRUE)
x50 <- sumyr(50) ;  hist(x50,ncl=100,prob=TRUE)
m50 <- mean(x50)
s50 <- sd(x50)
curve(dnorm(x,m50,s50),add=TRUE)
 
layout(1)   # rimettiamo a posto lo schermo