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 Ambos lados da revisão seguinte | ||
disciplinas:ce089-2014-02 [2014/10/13 18:53] walmes |
disciplinas:ce089-2014-02 [2014/11/10 21:22] walmes |
||
---|---|---|---|
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)) | ||
+ | |||
+ | </code> | ||
==== Links úteis ==== | ==== Links úteis ==== |