# gaussiana con il trucco delle 12 uniformi
 
my.norm <- function(n) {
    r <- rep(-6,n)
    for (i in 1:12) {
       r <- r + runif(n)
    }
    return(r)
}
 
x<- my.norm(100000)
hist(x,ncl=100,prob=TRUE)
mean(x)
sd(x)
 
# gaussiane a coppie
# se x e y sono Gaussiane standardizzate, modulo del raggio
# e' distribuito come f(r) = r e^(-r/2)
# [si ricorda che dx*dy = dr*(r*dphi) = r*dr*dphi ] 
# ovvero F(r) = 1-e^(-r^2/2).
# [Si tratta di una Rayleigh con parametro s=1.]
# => si estrae r e quindi phi uniforme e si ricostruiscono x e y
# ATT: nel seguito il raggio e' 'raggio', mentre r e' il vettore dei
# numeri casuali:
my.norm1 <- function(n) {
    twopi = 2 * pi
    n <- as.integer(n)           # per sicurezza
    if( n %% 2 ) dispari <- TRUE
    else dispari <- FALSE
    if (dispari) n <- n + 1      # se dispari abbondiamo
    m <- n %/% 2                 # numero di coppie
    raggio <- sqrt(-2*log(runif(m)))
    phi    <- runif(m, 0, twopi)
    r <- raggio * cos(phi) 
    r <- c(r, raggio * sin(phi) )
    if (dispari) r <- r[-n]
    return(r)
}
 
 
# per n non troppo grandi e' confrontabile in velocita' con rnorm()
# (ma per 10 milioni di valori ci mette circa 32 s, contro i 17 di rnorm,
# che e' scritta direttamente in C).
 
x<-my.norm1(10000000)           # per controllare le code a 4 e 5 sigma
length(x[abs(x)>4]) / 10000000  # -> 6.56e-05 (in un campione)
pnorm(-4)*2                     # -> 6.334248e-05
length(x[abs(x)>5]) / 10000000  # -> 5e-07    (in un campione)
pnorm(-5)*2                     # -> 5.733031e-07
 
# confrontiamo con quella piu' rozza: 
x<-my.norm(10000000)            # ci mette 3-4 minuti
length(x[abs(x)>4]) / 10000000  # -> 1.86e-05
length(x[abs(x)>5]) / 10000000  # -> 0      
# per code meno pronunciate:
length(x[abs(x)>3]) / 10000000  # -> 0.0020198
pnorm(-3)*2                     # -> 0.002699796    <= al limite!!
length(x[abs(x)>2]) / 10000000  # -> 0.0445292
pnorm(-2)*2                     # -> 0.04550026     <= quasi OK