# 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