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:ce089-2014-02 [2014/11/10 21:22] walmesdisciplinas: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)