# esempio di generatore con pdf data sqrt(1-x^2) nell'intervallo [1,1] # (istruzioni ridotte, assumendo familiare rlineare_hit-miss.R) n <- 100000 x <- runif(n) y <- runif(n) hit <- x^2 + y^2 < 1 xr <- x[hit] hist(xr, nc=100) nr <- length(xr) # nr <- sum(hit) plot(x, y, pch=19, cex=0.1, asp=1) points(x[hit], y[hit], pch=19, cex=0.1, col='red') # stima integrale int = nr / n # ma, siccome sappiamo che l'integrale vale pi/4, possiamo # usare il risultato della simulazine per "stimare" pi greco pi.mc <- int*4 # Ovviamente la stima di pi greco aumenta con il nr di estrazioni pir <- numeric() # vettore vuoto for (pn in 1:7) { n <- 10^pn print(n) x <- runif(n) y <- runif(n) pir[pn] <- length(x[x^2 + y^2 < 1]) / n * 4 } plot(10^(1:7), pir, log='x', xlab='n', ylab='pi') abline(h=pi, lty=2)