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

Ambos lados da revisão anterior Revisão anterior
Próxima revisão
Revisão anterior
disciplinas:ce089-2014-02 [2014/11/10 21:22]
walmes
disciplinas:ce089-2014-02 [2014/12/15 14:16] (atual)
walmes
Linha 22: Linha 22:
 ---- ----
  
 +/*
 ==== Histórico das Aulas do Curso ==== ==== Histórico das Aulas do Curso ====
  
Linha 67: Linha 67:
  
 ---- ----
 +
 +*/
  
 <code R> <code R>
Linha 137: Linha 139:
  
 plot(profile(m0)) plot(profile(m0))
 +
 +
 +##​--------------------------------------------
 +## Implementação conforme sugestão do Wikipedia.
 +
 +require(rootSolve)
 +
 +n <- rbinom(1, size=200, p=0.8)
 +y <- c(rpois(n, lambda=exp(2)),​ rep(0, 200-n))
 +length(y)
 +
 +barx <- mean(y)
 +p <- sum(y==0)/​length(y)
 +
 +f <- function(lambda,​ L){
 +    L$barx*(1-exp(-lambda))-lambda*(1-L$p)
 +}
 +
 +L <- list(barx=barx,​ p=p)
 +gradient(f, x=2, L=L)
 +
 +curve(f(x, L=L), 0, 15); abline(h=0)
 +
 +## Newton-Raphson.
 +maxiter <- 50; i <- 1          ## Número máximo de iterações e contador.
 +tol <- 1e-5; error <- 100*tol ​ ## Tolerância e erro inicial.
 +theta <- matrix(NA, nrow=1, ncol=1)
 +theta[1,] <- barx
 +while(i <= maxiter & error>​tol){
 +    theta <- rbind(theta,​ rep(NA,1))
 +    G <- f(theta[i,​],​ L)
 +    H <- gradient(f=f,​ x=theta[i,​],​ L=L)
 +    theta[i+1,] <- theta[i,​]-solve(H)%*%G
 +    error <- sum(abs((theta[i+1,​]-theta[i,​])/​theta[i+1,​]))
 +    i <- i+1
 +    ## print(c(theta[i,​],​ G))
 +    print(cbind(H,​ G, theta[i,]))
 +}
 +
 +lam <- theta[i,​] ​       ## Estimativa do lambda.
 +pii <- 1-barx/​theta[i,​] ## Estimativa do pi.
 +
 +##​-----------------------------------------------------------------------------
 +## Vendo os contornos da verossimilhança.
 +
 +llmax <- ll(th=c(log(pii/​(1-pii)),​ log(lam)), y=y)
 +
 +th1 <- seq(-7, 5, l=50)
 +th2 <- seq(-1, 4, l=50)
 +lla <- outer(th1, th2, llv, y=y)
 +
 +contour(th1,​ th2, lla,
 +        ## levels=seq(from=llmax-30,​ to=llmax, by=10),
 +        nlevels=20,
 +        xlab="​theta1:​ logit(p)",​
 +        ylab="​theta2:​ log(lambda)"​)
 +abline(v=log(pii/​(1-pii)),​ h=log(lam), col=2)
  
 </​code>​ </​code>​

QR Code
QR Code disciplinas:ce089-2014-02 (generated for current page)