Você está aqui: start » disciplinas » ce089-2014-02

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/09/16 17:00] – 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 53: Linha 53:
   - 15/09:   - 15/09:
     * Avaliação via simulação do poder de testes paramétricos e não paramétricos na presença e ausência dos pressupostos.     * Avaliação via simulação do poder de testes paramétricos e não paramétricos na presença e ausência dos pressupostos.
 +  - 17/09:
 +    * Avaliação via simulação do poder de testes paramétricos e não paramétricos na presença e ausência dos pressupostos.
 +  - 22/09:
 +    * Nível de cobertura de intervalos de confiança - intervalos baseados na deviance a na sua aproximação (Wald).
 +  - 24/09:
 +    * Nível de cobertura de intervalos de confiança - efeito da parametrização.
 +  - 29/09:
 +    * Desenvolvimento do trabalho 03.
 +  - 01/10:
 +    * Desenvolvimento do trabalho 03.
 +  - 13/10:
 +    * Método Newton-Raphson.
  
 ---- ----
 +
 +*/
 +
 +<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)