# 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