Você está aqui: start » disciplinas » ce718 » atividades2011 » pi

Diferenças

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

Link para esta página de comparações

Ambos lados da revisão anteriorRevisão anterior
Próxima revisão		Revisão anterior
disciplinas:ce718:atividades2011:pi [2011/06/10 09:25] – paulojusdisciplinas:ce718:atividades2011:pi [2011/06/21 00:56] (atual) – [section 4] fernandomayer
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>
 +
  

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