Estimação de pi por simulação

Diferentes abordagens para a estimativa de π através do experimento da agulha de Buffon, realizadas pelos participantes do curso:

buffon.eder <- function(n,l=1,a=1){
    if(a<l){cat('Erro: a < l, deve ser a > l\n')}
    if(a>=l){
        theta <- runif(n,0,pi)
        dist <- runif(n,0,a/2)
        inter <- sum(dist <= l/2*sin(theta))
        phi_est <- round((n/inter)*(2*l/a),12)
        cat('Número Simulação',n,'phi_estimado', phi_est,'Erro',
            round(pi-phi_est,12), '\n')
        return(c(n,phi_est))
    }
}
 
## Teste de funcionamento
buffon.eder(1000)
 
## Abordagem final
n <- seq(10000,1000000,by=20000)
res <- matrix(NA,ncol=2,nrow=length(n))
con <- 1
system.time(
            for (i in n){
                res[con,] <- buffon.eder(i)
                con <- con+1
            }
            )
 
## Plot
plot(res,type='l',ylab=expression(pi),xlab='Simulações')
abline(h=pi,col='red')

buffon.walmes <- function(simul, d=1, l=1, n=10){
    ## argumentos da função
    ## simul é o escalar inteiro número de agulhas lançadas
    ## d é o escalar espaçamento entre as linhas da grade
    ## l é o escalar comprimento da agulha
    ## n é o escalar inteiro número de linhas na grade
    malha <- seq(0, 10, by=d)
    linhas <- malha[-c(1,length(malha))]
    range <- range(malha)+c(1,-1)*d
    x <- matrix(runif(2*simul, range[1], range[2]), ncol=2)
    a <- runif(simul, 0, 2*pi)
    y <- x+matrix(c(l*sin(a), l*cos(a)), ncol=2)
    R <- sapply(seq_len(simul),
                function(i){
                    sum(linhas>=x[i,1] & linhas<=y[i,1])>0
                })
    rho <- l/d
    ## função retorna um a lista com
    ## R vetor de TRUE/FALSE se a linha toca ou não a grade
    ## o valor de rho
    ## a matrix X com as coordenadas onde a agulha toca
    return(list(R=R, rho=rho, X=cbind(cabeca=x,ponta=y)))
}
 
## Teste de funcionamento
buffon.walmes(1000)
 
## Abordagem final
bf <- buffon.walmes(1000)
pi0 <- bf$rho/(sum(bf$R)/length(bf$R)); pi0
 
## Plot
plot(bf$rho/(cumsum(bf$R)/c(1:length(bf$R))), type="l")
abline(h=pi)

buffon.fernando <- function(n, a = 1, l = 1){
    buffon.calc <- function(n, ...){
        D.random <- runif(n, 0, a/2)
        theta.random <- runif(n, 0, pi)
        d <- (l/2) * sin(theta.random)
        H <- numeric(n)
        H[D.random <= d] <- 1
        h <- sum(H)
        pi.est <- (2*l*n)/(a*h)
        return(pi.est)
    }
    pi.res <- sapply(n, buffon.calc)
    saida <- data.frame(n = n,
                        pi.est = pi.res,
                        abs.error = abs(pi.res - pi))
    class(saida) <- c("buffon", "data.frame")
    return(saida)
}
 
## Teste
system.time(
            testef <- buffon.fernando(1:1000)
            )
 
## Abordagem final
# ja finalizada
source("plot.buffon.R")
 
## Plot
plot(buffon.fernando(1:1e4))

A função plot.buffon() está aqui:

plot.buffon <- function(x, xlab = "Numero de jogadas da agulha",
                        ylab = expression(paste("Estimativa de ", pi)),
                        ...){
    plot(x$n, x$pi.est, type = "l", xlab = xlab, ylab = ylab,
         main = "Experimento de Buffon", ...)
    abline(h = pi, col = "lightgrey")
}

Para comparação, o código feito pelos autores da apostila (e também disponível no pacote CIM) está abaixo:

require(CIM)
buf(100)
 
## edita a funcao para tirar o plot
buffon.stuart <- function(n){
    ttt <- NULL
    ttt[1] <- 0
    x <- runif(n)
    th <- runif(n, 0, pi)
    st <- sin(th)
    for (i in 1:n) {
        if (st[i] > x[i])
            ttt[i + 1] <- ttt[i] + 1
        else ttt[i + 1] <- ttt[i]
    }
    #ttt
    if (ttt[n + 1] > 0)
        2 * (0:n)[ttt > 0]/ttt[ttt > 0]
    else print("no successes")
}
 
## Teste
buffon.stuart(1000)
 
## Plot
plot(buffon.stuart(1000), type = "l")
abline(h = pi)