[R-br] EMV
Walmes Zeviani
walmes em ufpr.br
Segunda Junho 23 12:24:49 BRT 2014
Ana Paula,
A optim() otimiza qualquer função. Nem precisa ser uma função de
verossimilhança. Ela pode ser usada para definir delineamentos ótimos,
máximo de funções como curvas de lactação, enfim, serve à propósitos
gerais. Sendo assim, depende do usuário escrever a função objetivo e
passar. Abaixo seguem exemplos da gama com a optim() e da weibull com
optim(), nlm(), survival::survreg() e bbmle::mle2().
##-----------------------------------------------------------------------------
## Gama.
n <- 100;
shape <- 5; scale <- 3
curve(dgamma(x, shape, scale), 0, 10)
y <- rgamma(n, shape=shape, scale=scale)
ll <- function(theta, y){
sum(dgamma(y, shape=theta[1], scale=theta[2], log=TRUE))
}
start <- c(5,3)
op <- optim(start, fn=ll, y=y, hessian=TRUE,
method="BFGS", control=list(fnscale=-1))
op$par
op$value
op$hessian
##-----------------------------------------------------------------------------
## Weibull.
n <- 100;
shape <- 5; scale <- 3
curve(dweibull(x, shape, scale), 0, 10)
y <- rweibull(n, shape=shape, scale=scale)
ll <- function(theta, y){
sum(dweibull(y, shape=theta[1], scale=theta[2], log=TRUE))
}
start <- c(5,3)
op <- optim(start, fn=ll, y=y, hessian=TRUE,
method="BFGS", control=list(fnscale=-1))
op$par
op$value
op$hessian
##-----------------------------------------------------------------------------
## Usando a nlm (só minimiza).
ll <- function(theta, y){
-sum(dweibull(y, shape=theta[1], scale=theta[2], log=TRUE))
}
start <- c(5,3)
op <- nlm(p=start, f=ll, y=y, hessian=TRUE)
op
##-----------------------------------------------------------------------------
## Usando a survreg (para a gama seria a glm(..., family=gamma)).
require(survival)
m0 <- survreg(Surv(time=y)~1, dist="weibull")
summary(m0)
exp(coef(m0)) ## = scale.
1/m0$scale ## = shape.
##-----------------------------------------------------------------------------
## Usando a bblme.
require(bbmle)
## Forma de especificar I.
ml0 <- mle2(y~dweibull(shape=sh, scale=sc),
start=list(sh=5, sc=3), data=list(y=y))
class(ml0)
## Vantagem é a disponibilidade de métodos.
summary(ml0)
vcov(ml0)
plot(profile(ml0))
confint(ml0)
confint(ml0, method="quad")
## Forma de específicar II.
ll <- function(shape, scale, y){
-sum(dweibull(y, shape=shape, scale=scale, log=TRUE))
}
start <- list("shape"=5, "scale"=3)
ml1 <- mle2(minuslogl=ll, start=start, method="BFGS", data=list(y=y))
summary(ml1)
vcov(ml1)
plot(profile(ml1))
confint(ml1)
confint(ml1, method="quad")
##-----------------------------------------------------------------------------
À
disposição
Walmes.
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20140623/cb2e13ec/attachment.html>
Mais detalhes sobre a lista de discussão R-br