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 | Última revisão Ambos lados da revisão seguinte | ||
|
disciplinas:ce089-2014-02 [2014/11/10 21:22] walmes |
disciplinas:ce089-2014-02 [2014/11/10 22:04] walmes |
||
|---|---|---|---|
| Linha 137: | Linha 137: | ||
| 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> | ||
