[R-br] [Dúvida] Hessiana pela função optim

[Wagner Bonat]

Esse problema pode ocorrer por muitos motivos incluindo apenas uma má
aproximação do hessiano.
Eu recomendo algumas coisas:
1 - Calcule a verossimilhança perfilhada de cada parâmetro. Isso vai de dar
certeza que o algoritmo convergiu e te mostrar a cara da função. De acordo
com o comportamento da profile likelihood vc pode tentar alguma
reparametrização pra melhorar o comportamento do algoritmo numerico. Eu
usaria o pacote bbmle que vai fazer isso automaticamente pra vc.
2 -  Calcule ao menos a primeira derivada da sua log-verossimilhança e use
no BFGS.
3 - Calcule o hessiano numericamente porém separado, vc pode usar por
exemplo o pacote numDeriv.

Att

Em 4 de agosto de 2016 22:39, Pedro Rafael via R-br <
r-br em listas.c3sl.ufpr.br> escreveu:

> Caros,
>
> sabemos que a matriz Hessiana multiplicada por -1 da função de
> log-verossimilhança nos fornecem a matriz de informação observada que
> converge assintoticamente para a matriz de informação esperada. A diagonal
> principal da matriz inversa da informação esperada nos dá as variâncias os
> estimadores de máxima verossimilhança.
>
> Um fato importante para o o questionamento que vou fazer é que estou
> fazendo uso da matriz de informação observada de modo a ter ao menos uma
> aproximação da variância dos estimadores de máxima verossimilhança obtidos
> numericamente pelo método BFGS.
>
> O problema é que estou obtendo em alguns caros que a
> solve(-diag(hessiana)) negativa o que não era de se esperar. Vejam se
> possível o código.
>
> # PDF
> pdf_ekww <- function(par,x){
>   a = par[1]
>   b = par[2]
>   c = par[3]
>   alpha = par[4]
>   beta = par[5]
>   g = dweibull(x = x, shape = alpha, scale = beta, log = FALSE)
>   G = pweibull(q = x, shape = alpha, scale = beta, lower.tail = TRUE,
> log.p = FALSE)
>   a * b * c * g * G^(a-1) * (1-G^a)^(b-1) * (1-(1-G^a)^b)^(c-1)
> }
>
> # Quantilica
> sample_ekww <- function(n,par){
>   a = par[1]
>   b = par[2]
>   c = par[3]
>   alpha = par[4]
>   beta = par[5]
>   u = runif(n=n,min=0, max=1)
>   qweibull(p = (1-(1-u^(1/c))^(1/b))^(1/a), shape = alpha, scale = beta,
> lower.tail = TRUE, log.p = FALSE)
> }
> set.seed(1987)
>
> vector_par_true = c(1,1,1,1,1)
> data = sample_ekww(n = 1000, par = vector_par_true)
>
> # Função de log-verossimilhança.
> obj_ekww = function(par,x){
>   sum(log(pdf_ekww(par,x)))
> }
>
> # Maximizando log-verossimilhança.
> result = optim(fn = obj_ekww, par = c(0.5,1.4,1,1,1), method = "BFGS", x =
> data, hessian = TRUE, control=list("fnscale"=-1))
>
> diag(solve(-result$hessian))
>
> Observem que o primeiro elemento do vetor logo acima é negativo o que não
> deveria ser verdade. Notem também que houve convergência segundo o critério
> de parada da função optim, em que convergence é igual a zero. O que para
> vocês podem estar provocando esse problema? A função objetivo é complicada
> a ponto de provocar problemas na convergência do algoritmo vindo por sua
> vez acarretar esse tipo de problema?
>
> Digo isso porque essas novas classes de distribuições de probabilidade por
> apresentar diversos parâmetros produzem log-verossimilhanças muito
> complicadas incluindo problemas de regiões aproximadamente planas bem como
> problemas piores como log-verossimilhanças monótonas.
>
> Obrigado desde já,
> Pedro Rafael.
>
>
> _______________________________________________
> R-br mailing list
> R-br em listas.c3sl.ufpr.br
> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça
> código mínimo reproduzível.
>



-- 
Wagner Hugo Bonat
----------------------------------------------------------------------------------------------
Department of Mathematics and Computer Science (IMADA)
University of Southern Denmark (SDU) and
Laboratório de Estatística e Geoinformação (LEG)
Universidade Federal do Paraná (UFPR)
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20160805/4139cbe0/attachment.html>