Ajustes
Distribuição Normal
Função de log-verossimilhança
# ll.n: log-verossimilhança da normal
ll.n <- function(mu, sigma.l, y){
saida <- sum(dnorm(y,
mean = mu,
sd = exp(sigma.l),
log = TRUE))
return(-saida)}Estimando os parâmetros e a máxima verossimilhança
# aval.ll.n: avaliando a função de log-verossimilhança da normal
aval.ll.n <- mle2(ll.n,
start = list(mu = mean(y),
sigma.l = log(sd(y))),
method = "BFGS",
data = list(y = y))Estimativas
list(mu = coef(aval.ll.n)[[1]],
sigma = exp(coef(aval.ll.n)[[2]]))$mu
[1] 9.125
$sigma
[1] 4.883531Comparando com a função fitdistr()
fitdistr(y, "normal")$estimate mean sd
9.125000 4.883531 Estimativa da verossimilhança
logLik(aval.ll.n)'log Lik.' -60.09614 (df=2)Comparando com a função fitdistr()
fitdistr(y, "normal")$loglik[1] -60.09614Intervalos de confiança (IC)
list(mu.ic = confint(aval.ll.n,
method = "quad")[1, ],
sigma.ic = exp(confint(aval.ll.n,
method = "quad")[2, ]))$mu.ic
2.5 % 97.5 %
6.984738 11.265262
$sigma.ic
2.5 % 97.5 %
3.582178 6.657646 Perfis de verossimilhança
perf.n <- profile(aval.ll.n)
plot(perf.n,
xlab = c(expression(mu), expression(log(sigma))))
Intervalos de confiança (IC)
list(mu.ic = confint(perf.n)[1, ],
sigma.ic = exp(confint(perf.n)[2, ]))$mu.ic
2.5 % 97.5 %
6.877702 11.372298
$sigma.ic
2.5 % 97.5 %
3.687679 6.899465 Distribuição Normal com transformação na resposta
Estimando os parâmetros e a máxima verossimilhança
# aval.ll.nt: avaliando a função de log-verossimilhança da normal com resposta transformada
aval.ll.nt <- mle2(ll.n,
start = list(mu = mean(y.bc),
sigma.l = log(sd(y.bc))),
method = "BFGS",
data = list(y = y.bc))Estimativa da verossimilhança
logLik(aval.ll.nt)'log Lik.' -26.40036 (df=2)Comparando com a função fitdistr()
fitdistr(y.bc, "normal")$loglik[1] -26.40036Tornando a verossimilhança comparável
\[ l(\theta; y) = l(\theta; y^{*}) - \sum_{i=1}^{n} \log(J_{i}), \]
onde
\(l(\theta; y)\) é a log-verossimilhança na escala da variável original;
\(l(\theta; y^{*})\) é a log-verossimilhança do modelo ajustado com a variável transformada;
\(\sum_{i=1}^{n} \log(J_{i})\) é a soma dos log-Jacobianos para as observações individuais.
Aqui o log-Jacobiano utilizado é dado por:
\[ \log(J) = \left\{\begin{array}{ll} (\frac{1}{\lambda}-1) \log(\lambda \cdot Y^{*} + 1) & \mbox{ se } \lambda \neq 0 \\ Y^{*} & \mbox{ se } \lambda = 0 \\ \end{array}\right.. \]
# J.log: log-Jacobianos
J.log <- (1/(lambda) - 1) * log(lambda * y.bc + 1)
(aval.ll.nt <- logLik(aval.ll.nt) - sum(J.log))'log Lik.' -59.12448 (df=2)Distribuição Log-Normal
Função de log-verossimilhança
# ll.ln: log-verossimilhança da log-normal
ll.ln <- function(mu.log, sigma.log, y){
saida <- sum(dlnorm(y,
meanlog = exp(mu.log),
sdlog = exp(sigma.log),
log = TRUE))
return(-saida)}Estimando os parâmetros e a máxima verossimilhança
# aval.ll.ln: avaliando a função de log-verossimilhança da log-normal
aval.ll.ln <- mle2(ll.ln,
start = list(mu.log = log(mean(y)),
sigma.log = log(sd(y))),
method = "BFGS",
data = list(y = y))Estimativas
exp(coef(aval.ll.ln)) mu.log sigma.log
2.0452607 0.6416122 Comparando com a função fitdistr()
fitdistr(y, "lognormal")$estimate meanlog sdlog
2.045257 0.641607 Estimativa da verossimilhança
logLik(aval.ll.ln)'log Lik.' -60.40833 (df=2)Comparando com a função fitdistr()
fitdistr(y, "lognormal")$loglik[1] -60.40833Intervalos de confiança (IC)
exp(confint(aval.ll.ln,
method = "quad")) 2.5 % 97.5 %
mu.log 1.7825409 2.3467014
sigma.log 0.4706356 0.8747027Perfis de verossimilhança
perf.ln <- profile(aval.ll.ln)
plot(perf.ln,
xlab = c(expression(log(mu)), expression(log(sigma))))
Intervalos de confiança (IC)
exp(confint(perf.ln)) 2.5 % 97.5 %
mu.log 1.7499072 2.340571
sigma.log 0.4844939 0.906464Distribuição Gamma
Função de log-verossimilhança
# ll.g: log-verossimilhança da gamma
ll.g <- function(forma.l, escala.l, y){
saida <- sum(dgamma(y,
shape = exp(forma.l),
scale = exp(escala.l),
log = TRUE))
return(-saida)}Estimando os parâmetros e a máxima verossimilhança
Como definir chutes iniciais?
Vamos chamar o parâmetro de forma de \(k\), e o parâmetro de escala de \(\theta\).
A esperança (média) e variância da variável aleatória são dados por:
\(E[X] = k \cdot \theta\);
\(Var[X] = k \cdot \theta^{2}\).
Logo,
\[ \text{média} = k \cdot \theta \quad \rightarrow \quad \theta = \frac{\text{média}}{k} \]
\[ \text{variância} = k \cdot \theta^{2} \rightarrow \text{variância} = k \cdot \left( \frac{\text{média}}{k} \right)^{2} = \frac{\text{média}^{2}}{k} \]
e
\[ k = \frac{\text{média}^{2}}{\text{variância}} \]
k <- (mean(y)**2) / var(y)
theta <- mean(y) / k
# aval.ll.g: avaliando a função de log-verossimilhança da gamma
aval.ll.g <- mle2(ll.g,
start = list(forma.l = log(k),
escala.l = log(theta)),
method = "BFGS",
data = list(y = y))Estimativas
list(forma = exp(coef(aval.ll.g)[[1]]),
escala = exp(coef(aval.ll.g)[[2]]))$forma
[1] 3.173311
$escala
[1] 2.875546Comparando com a função fitdistr()
fitdistr(y, "gamma")$estimate shape rate
3.1733266 0.3477619 # A fitdistr() usa o parâmetro rate, que é igual a 1 / scale
1 / exp(coef(aval.ll.g)[[2]])[1] 0.3477601Estimativa da verossimilhança
logLik(aval.ll.g)'log Lik.' -58.77988 (df=2)Comparando com a função fitdistr()
fitdistr(y, "gamma")$loglik[1] -58.77988Intervalos de confiança (IC)
list(forma.ic = exp(confint(aval.ll.g,
method = "quad")[1, ]),
escala.ic = exp(confint(aval.ll.g,
method = "quad")[2, ]))$forma.ic
2.5 % 97.5 %
1.758834 5.725330
$escala.ic
2.5 % 97.5 %
1.517232 5.449902 Perfis de verossimilhança
perf.g <- profile(aval.ll.g)
plot(perf.g,
xlab = c("forma ( k )", expression("escala ("~theta~")")))
Intervalos de confiança (IC)
list(forma = exp(confint(perf.g)[1, ]),
escala = exp(confint(perf.g)[2, ]))$forma
2.5 % 97.5 %
1.663249 5.449225
$escala
2.5 % 97.5 %
1.617788 5.924838