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)
## Define um n comum
n1 <- 1:1e+4
n2 <- seq(0, 1e+6, 1000)[-1]
##----------------------------------------------------------------------
## Tempo Eder com n1
res.eder.n1 <- matrix(NA, ncol=2, nrow=length(n1))
con <- 1
eder.n1 <- system.time(
for (i in n1){
res.eder.n1[con,] <- buffon.eder(i)
con <- con+1
}
)
## Tempo Walmes com n1
walmes.n1 <- system.time(
bf <- buffon.walmes(max(n1))
)
res.walmes.n1 <- bf$rho/(cumsum(bf$R)/c(1:length(bf$R)))
## Tempo Fernando com n1
fernando.n1 <- system.time(
res.fernando.n1 <- buffon.fernando(n1)
)
## Tempo Stuart com n1
stuart.n1 <- system.time(
res.stuart.n1 <- buffon.stuart(max(n1))
)
##----------------------------------------------------------------------
## Tempo Eder com n2
res.eder.n2 <- matrix(NA, ncol=2, nrow=length(n2))
con <- 1
eder.n2 <- system.time(
for (i in n2){
res.eder.n2[con,] <- buffon.eder(i)
con <- con+1
}
)
## Tempo Walmes com n2
walmes.n2 <- system.time(
bf <- buffon.walmes(max(n2))
)
res.walmes.n2 <- bf$rho/(cumsum(bf$R)/c(1:length(bf$R)))
## Tempo Fernando com n2
fernando.n2 <- system.time(
res.fernando.n2 <- buffon.fernando(n2)
)
## Tempo Stuart com n2
stuart.n2 <- system.time(
res.stuart.n2 <- buffon.stuart(max(n2))
)
#### Parei em 3261 seg. ~ 55 min.
## Usando n1 = 1:1e4
##----------------------------------------------------------------------
## Comparacao de tempos de execucao
tempo.n1 <- c(eder.n1[3], walmes.n1[3], fernando.n1[3], stuart.n1[3])
names(tempo.n1) <- c("eder", "walmes", "fernando", "stuart")
tempo.n1 <- sort(tempo.n1)
## barchart
require(lattice)
barchart(tempo.n1, xlab = "Tempo (s)",
panel = function(...){
panel.barchart(...)
panel.text(x = tempo.n1, y = 1:4,
labels = do.call(as.character,
list(round(tempo.n1, 2))), pos = 2)
})
## Cria um data.frame com todos os resultados
res.n1 <- data.frame(n = n1,
eder = res.eder.n1[,2],
walmes = res.walmes.n1,
fernando = res.fernando.n1$pi.est,
stuart = res.stuart.n1)
## modifica o data.frame para o lattice
require(reshape)
res.n1 <- melt(res.n1, id = 1)
## Comparacao grafica
xyplot(value ~ n | variable, data = res.n1, type = "l", as.table = TRUE,
xlab = "Número de jogadas da agulha",
ylab = expression(paste("Estimativa de ", pi)),
panel = function(...){
panel.xyplot(...)
panel.abline(h = pi, lty = 2)
}, scales = list(relation = "free"))
## Usando n2 = seq(0, 1e+6, 1000)[-1]
##----------------------------------------------------------------------
## Comparacao de tempos de execucao
tempo.n2 <- c(eder.n2[3], walmes.n2[3], fernando.n2[3])
names(tempo.n2) <- c("eder", "walmes", "fernando")
tempo.n2 <- sort(tempo.n2)
## barchart
barchart(tempo.n2, xlab = "Tempo (s)",
panel = function(...){
panel.barchart(...)
panel.text(x = tempo.n2, y = 1:4,
labels = do.call(as.character,
list(round(tempo.n2, 2))), pos = 2)
})
## Cria um data.frame com todos os resultados
## Stuart nao entra pq nao terminou a execução
## Walmes fica de fora pq vai ate 1e6
res.n2 <- data.frame(n = n2,
eder = res.eder.n2[,2],
fernando = res.fernando.n2$pi.est)
## modifica o data.frame para o lattice
res.n2 <- melt(res.n2, id = 1)
## Comparacao grafica
xyplot(value ~ n | variable, data = res.n2, type = "l", as.table = TRUE,
xlab = "Número de jogadas da agulha",
ylab = expression(paste("Estimativa de ", pi)),
panel = function(...){
panel.xyplot(...)
panel.abline(h = pi, lty = 2)
})#, scales = list(relation = "free"))
## Plot separado do Walmes
xyplot(res.walmes.n2 ~ 1:max(n2), type = "l",
xlab = "Número de jogadas da agulha",
ylab = expression(paste("Estimativa de ", pi)),
panel = function(...){
panel.xyplot(...)
panel.abline(h = pi, lty = 2)
})