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/10/13 18:53]
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>
 +##​-----------------------------------------------------------------------------
 +
 +require(bbmle)
 +
 +##​-----------------------------------------------------------------------------
 +
 +y <- rnorm(100, mean=0, sd=1)
 +crt <- 1.5
 +r <- y<crt
 +y[!r] <- crt
 +r <- as.integer(r)
 +plot(ecdf(y))
 +table(r)
 +
 +ll <- function(theta,​ y, r){
 +    ys <- split(y, r)
 +    ll1 <- dnorm(ys[["​1"​]],​
 +                 ​mean=theta[1],​
 +                 ​sd=exp(theta[2]),​
 +                 ​log=TRUE)
 +    ll2 <- pnorm(ys[["​0"​]],​
 +                 ​mean=theta[1],​
 +                 ​sd=exp(theta[2]),​
 +                 ​lower.tail=FALSE,​
 +                 ​log=TRUE)
 +    ll12 <- sum(ll1)+sum(ll2)
 +    ## Tem que retornar o negativo da soma. A mle2() minimiza.
 +    return(-ll12)
 +}
 +
 +##​-----------------------------------------------------------------------------
 +## Estimação com a bbmle.
 +
 +parnames(ll) <- c("​mu",​ "​lsigma"​)
 +m0 <- mle2(minuslogl=ll,​ start=list(mu=0,​ lsigma=log(1)),​
 +           ​vecpar=TRUE,​ data=list(y=y,​ r=r), method="​BFGS"​)
 +
 +coef(m0)
 +summary(m0)
 +ci <- confint(m0, method="​quad"​)
 +ci
 +
 +htheta <- coef(m0)
 +logL <- c(logLik(m0))
 +
 +##​-----------------------------------------------------------------------------
 +
 +llv <- Vectorize(
 +    FUN=function(th1,​ th2, y, r){
 +        -ll(c(th1, th2), y=y, r=r)
 +    },
 +    vectorize.args=c("​th1",​ "​th2"​))
 +
 +th1 <- seq(-0.4, 0.4, l=30)
 +th2 <- seq(-0.1, 0.4, l=30)
 +
 +lla <- outer(th1, th2, llv, y=y, r=r)
 +
 +contour(th1,​ th2, lla,
 +        xlab=expression(mu),​
 +        ylab=expression(log(sigma)))
 +abline(v=htheta[1],​ h=htheta[2],​ lty=2)
 +contour(th1,​ th2, lla, add=TRUE,
 +        levels=(logL-0.5*qchisq(0.95,​ df=length(htheta))),​
 +        col=2)
 +abline(v=ci[1,​],​ h=ci[2,], lty=3, col=2)
 +
 +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>​
  
 ==== Links úteis ==== ==== Links úteis ====

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