Não foi possível enviar o arquivo. Será algum problema com as permissões?
Diferenças
Aqui você vê as diferenças entre duas revisões dessa página.
| 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 ==== | ||
