Não foi possível enviar o arquivo. Será algum problema com as permissões?
Diferenças

Diferenças

Aqui você vê as diferenças entre duas revisões dessa página.

Link para esta página de comparações

Próxima revisão
Revisão anterior
disciplinas:ce718:atividades2011:pi [2011/06/10 09:23]
paulojus criada
disciplinas:ce718:atividades2011:pi [2011/06/21 00:56] (atual)
fernandomayer [section 4]
Linha 1: Linha 1:
-====== Estimação de <m>pi</​m> ​======+====== Estimação de pi por simulação ​======
  
-===== Simulação em quadrado/​circulo ​=====+==== Simulação em quadrado/​circulo ====
  
 ==== Agulha de Buffon ==== ==== Agulha de Buffon ====
  
-=== Comparação dos algorítmos ===+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 ===
 + 
 +=== 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>​ 
  

QR Code
QR Code disciplinas:ce718:atividades2011:pi (generated for current page)