[Dúvida] Hessiana pela função optim
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.
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@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@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)
Obrigado Wagner pela dica. Não conhecia o pacote bbmle. Na verdade sempre fui tentado a sempre implementar tudo. Porém, seria uma boa tentar utilizar o pacote. Att. Em 5 de agosto de 2016 09:39, Wagner Bonat <wbonat@gmail.com> escreveu:
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@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@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)
participantes (2)
-
Pedro Rafael -
Wagner Bonat