Não foi possível enviar o arquivo. Será algum problema com as permissões?
Diferenças
Aqui você vê as diferenças entre duas revisões dessa página.
Ambos lados da revisão anterior Revisão anterior Próxima revisão | Revisão anterior | ||
disciplinas:ce718:atividades2011:pi [2011/06/10 09:25] paulojus |
disciplinas:ce718:atividades2011:pi [2011/06/21 00:56] (atual) fernandomayer [section 4] |
||
---|---|---|---|
Linha 4: | Linha 4: | ||
==== Agulha de Buffon ==== | ==== Agulha de Buffon ==== | ||
+ | |||
+ | Diferentes abordagens para a estimativa de π através do experimento da agulha de Buffon, realizadas pelos participantes do curso: | ||
+ | |||
+ | === Éder === | ||
+ | |||
+ | <code R> | ||
+ | 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') | ||
+ | </code> | ||
+ | |||
+ | === Walmes === | ||
+ | |||
+ | <code R> | ||
+ | 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) | ||
+ | </code> | ||
+ | |||
+ | === Fernando === | ||
+ | |||
+ | <code R> | ||
+ | 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)) | ||
+ | </code> | ||
+ | |||
+ | A função ''plot.buffon()'' está aqui: | ||
+ | |||
+ | <code R> | ||
+ | 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") | ||
+ | } | ||
+ | </code> | ||
+ | |||
+ | === Stuart === | ||
+ | |||
+ | Para comparação, o código feito pelos autores da apostila (e também disponível no pacote CIM) está abaixo: | ||
+ | |||
+ | <code R> | ||
+ | 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) | ||
+ | </code> | ||
+ | |||
==== Comparação dos algorítmos ==== | ==== Comparação dos algorítmos ==== | ||
+ | |||
+ | === Executando e armazenando tempos e resultados === | ||
+ | |||
+ | <code R> | ||
+ | ## 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. | ||
+ | </code> | ||
+ | |||
+ | === Comparando performances === | ||
+ | |||
+ | <code R> | ||
+ | ## 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) | ||
+ | }) | ||
+ | </code> | ||
+ | |||