[R-br] Fwd: Não consigo estimar por máxima verossimilhança (problemas com logaritmo de argumento negativo)

[Romero Luiz M. Sales Filho]

Prezados, mais uma vez peço desculpas... Não consigo responder diretamente
ao tópico dando o reply na mensagem que chega até a minha caixa de e-mail.
Será que estou fazendo algum procedimento errado???

Como preciso muito solucionar meu caso, estou reenviando uma nova mensagem
com o mesmo título do tópico anteriormente criado. Se existir a
possibilidade de mesclar esse tópico com os outros dois criados eu ficaria
agradecido. Se alguém puder me orientar tbm em como responder uma mensagem
no R-BR eu agradeço, pois já participei de alguns tópicos aqui e acredito
que sempre respondia ao e-mail que chegava na minha conta...  Neste tópico
vou colocar as mensagens postadas anteriormente para pelo menos seguir uma
linha de raciocínio, e novamente peço para me orientarem como responder
nesse tópico atual. Obrigado.

###########################################

A seguir peço que leiam a situação atual do problema:

###########################################

Caros colegas,

estou com um problema que não consigo resolver... Na verdade acredito que
já identifiquei o problemas, mas não sei como contornar esse tipo de
situação no R. É o seguinte estou maximizando a função de
log-verossimilhança de uma dist. de probabilidade, com o intuito de estimar
os parâmetros da dist. Entretanto a otimização não está convergindo para o
parâmetro conhecido. Essa dist. possui dois parâmetros, onde um deles o "z"
assume valores menores que 2, e o outro parâmetro "eta" assume valores
maiores que 0. Consegui simular situações para os casos que tenho 1<z<2 e
para 0<z<1, entretanto para valores negativos de z, não estou conseguindo
convergência para o parâmetro. Estou usando a função optim() do R,
selecionando nessa função o método de Nelder-Mead (já usei outros métodos e
nenhum convergiu). O problema que ocorre é que para z negativo, durante as
iterações ocorrem o cálculo de logaritmos cujos argumentos estão com
valores negativos, tenho quase certeza que o problema da não convergência
está aí, aí eu gostaria de saber se tem como colocar alguma restrição
durante o processo de iteração para que os argumentos das funções
logarítmicas nunca sejam valores negativos. A seguir coloco o código que
estou usando:

##############################################

##### Fornecendo os parâmetros
##############################
n<-50
z<--1.95
eta<-0.1
z.est=vector()
eta.est=vector()

p=1
while (p<1001) {

###### Gerando valores da varíavel cujos parâmetros são z e eta

u<-runif(n)
y = eta*((1-(u)^((1-z)/(2-z)))/(1-z))

#### Obtendo as estimativas para z e eta
########################################

vero <- function(par,x){
z = par[1]
eta = par[2]
saida<-((sum(log(1 - ((1 - z)*x*(1/eta)))))/(1 - z)) + ((n)*log(2 - z)) +
(n*log(1/eta))
return(-saida)
}

saida1<-optim(par=c(-1,4),fn=vero,
      method="Nelder-Mead",x=y
      )
z.est[p]<-saida1[1]$par[1]
eta.est[p]<-saida1[1]$par[2]
p=p+1
}
mean(z.est)
mean(eta.est)

##############################################

Como estou replicando o procedimento mil vezes, irei utilizar a média das
1000 estimativas encontradas para z e para eta... Então fazendo isso tenho
a seguinte saída:

> mean(z.est)
[1] -30177882

> mean(eta.est)
[1] 964132.7

Conforme observamos as estimativas estão muito distantes dos parâmetros
reais definidos no início do código (z = -1.95 e eta = 0.1).


Dessa forma, desconfio que o problema seja decorrente dos vários erros de
logaritmos, que se repetem durante a interação e são da seguinte forma:

1: In log(1 - ((1 - z) * x * (1/eta))) : NaNs produzidos


Então tem alguma forma de controlar esse problema que ocorre durante a
interação??? Tenho que fazer de alguma forma com que os argumentos fiquem
sempre positivos, ou seja: (1-((1-z)*x*(1/eta))) tem que ser sempre
positivo...

Por favor, caso vcs possam me ajudar, eu agradeço imensamente.

Att,

Romero.






###########################################
###########################################
Mensagem anterior enviada por Wagner bonat:
**************************************


Uma outra opcao é vc restringir o espaco de busca da optim() usando o
metodo L-BFGS-B. Vc precisa identificar para quais valores do seu parametro
o erro acontece e se certificar que o algoritmo nao caminhe nesta direcao.

Uma segunda opcao é se for possivel calcule pelo menos a primeira derivada
da funcao de log-verossimilhanca e a forneca ao algoritmo BFGS ou L-BFGS-B
isso ajuda o algoritmo ficar sempre no espaco parametrico, porque a
derivada nao deixa ele tentar valores meio que absurdos.

Para uma implementacao definitiva eu aconselharia usar um Fisher-score
muito eficiente e normalmente partindo de valores iniciais razoaveis vc
sempre converge. 'E claro que para isso vc precisa calcular a matriz de
derivadas segundas e tomar a esperanca, dependendo da situacao pode ser
razoavel.


###########################################
###########################################
Mensagem anterior enviada por Romero (rlmsf)
**************************************


Olá Walmes e Vinícius,

obrigado pelas dicas, pelo menos já estou buscando encontrar alguma forma
de resolver, pois confesso, que já não estava vendo mais como sair do
dilema.

Vinicius, já tinha tentado fazer a otimização da função de verossimilhança
(sem o log), entretanto, não sei o que acontece, mas as estimativas ficam
todas iguais aos pontos iniciais, ou seja, para as mil replicações que
estou fazendo, todos os resultados são iguais aos pontos iniciais... Com
certeza, fazendo dessa forma, não deve estar chegando nem na segunda
iteração (não sei o motivo).

Walmes, a sua dica é boa, alguém já tinha falado em algo do tipo, mas eu
não tinha entendido bem... acho que agora com sua explicação ficou mais
claro, resumindo, por favor veja se é isso: preciso encontrar funções cujos
possíveis resultado sejam pertencentes ao domínio dos meus dois parâmetros
de interesse, com a exigência de que o argumento dessas funções sejam
sempre valores pertencentes aos reais, certo??

Pensando assim, como o parâmetro z é sempre menor que 2 e o parâmetro eta é
sempre um número positivo, então teoricamente posso usar a seguinte
transformação:

z= 2-exp(x) (argumento x real e resultado da função sempre menor que 2)   e
   eta=exp (k) (argumento k real e resultado da função sempre maior que
zero). É isso mesmo, vcs concordam???

Obrigado pela atenção,

Romero.


###########################################
###########################################
Mensagem anterior enviada por Vinicius Brito
**************************************


Romero,

faz tempo que não tenho esse tipo de problema na minha vida.

A sugestão do Walmes me parece a mais correta.

Outra abordagem que me parece ser razoável e tentar maximizar sem a função
de log, pois deveria obter o mesmo resultado.

Ou também usar uma abordagem bayesiana. Determinar uma distribuição a
posteriori e então determinar as distribuições condicionais completas para
os seus parâmetros. Sem dúvida esse tipo de approach vai te dar muito mais
trabalho.

Abs

Vinicius


###########################################
###########################################
Mensagem anterior enviada por Walmes
**************************************

 Em problemas de otimização o primeiro conselho é ter uma função objetivo
com argumentos no domínio dos reais. Tente colocar reparametrizações na sua
função objetivo para garantir isso. Por exemplo, no caso da distribuição
normal a função de log-verrossimilhança tem que ser escrita com mu=mu e
sigma=exp(theta), não se estima diretamente sigma pois sigma>0 mas
exp(theta)>0 para todos theta nos reais. A mesma idéia serve para binomial,
onde p está no (0,1), não se otimiza em p mas em theta sendo que
p=1/(1+exp(theta)) que vem da função de ligação log(p/(1-p))=theta.

À disposição.
Walmes.


###################################################
###################################################
Mensagem anterior enviada por Romero (rlmsf) (1º Tópico)
*********************************************

Caros colegas,

estou com um problema que não consigo resolver... Na verdade acredito que
já identifiquei o problemas, mas não sei como contornar esse tipo de
situação no R. É o seguinte estou maximizando a função de
log-verossimilhança de uma dist. de probabilidade, com o intuito de estimar
os parâmetros da dist. Entretanto a otimização não está convergindo para o
parâmetro conhecido. Essa dist. possui dois parâmetros, onde um deles o "z"
assume valores menores que 2, e o outro parâmetro "eta" assume valores
maiores que 0. Consegui simular situações para os casos que tenho 1<z<2 e
para 0<z<1, entretanto para valores negativos de z, não estou conseguindo
convergência para o parâmetro. Estou usando a função optim() do R,
selecionando nessa função o método de Nelder-Mead (já usei outros métodos e
nenhum convergiu). O problema que ocorre é que para z negativo, durante as
iterações ocorrem o cálculo de logaritmos cujos argumentos estão com
valores negativos, tenho quase certeza que o problema da não convergência
está aí, aí eu gostaria de saber se tem como colocar alguma restrição
durante o processo de iteração para que os argumentos das funções
logarítmicas nunca sejam valores negativos. A seguir coloco o código que
estou usando:

##############################################

##### Fornecendo os parâmetros
##############################
n<-50
z<--1.95
eta<-0.1
z.est=vector()
eta.est=vector()

p=1
while (p<1001) {

###### Gerando valores da varíavel cujos parâmetros são z e eta

u<-runif(n)
y = eta*((1-(u)^((1-z)/(2-z)))/(1-z))

#### Obtendo as estimativas para z e eta
########################################

vero <- function(par,x){
z = par[1]
eta = par[2]
saida<-((sum(log(1 - ((1 - z)*x*(1/eta)))))/(1 - z)) + ((n)*log(2 - z)) +
(n*log(1/eta))
return(-saida)
}

saida1<-optim(par=c(-1,4),fn=vero,
      method="Nelder-Mead",x=y
      )
z.est[p]<-saida1[1]$par[1]
eta.est[p]<-saida1[1]$par[2]
p=p+1
}
mean(z.est)
mean(eta.est)

##############################################

Como estou replicando o procedimento mil vezes, irei utilizar a média das
1000 estimativas encontradas para z e para eta... Então fazendo isso tenho
a seguinte saída:

> mean(z.est)
[1] -30177882

> mean(eta.est)
[1] 964132.7

Conforme observamos as estimativas estão muito distantes dos parâmetros
reais definidos no início do código (z = -1.95 e eta = 0.1).


Dessa forma, desconfio que o problema seja decorrente dos vários erros de
logaritmos, que se repetem durante a interação e são da seguinte forma:

1: In log(1 - ((1 - z) * x * (1/eta))) : NaNs produzidos


Então tem alguma forma de controlar esse problema que ocorre durante a
interação??? Tenho que fazer de alguma forma com que os argumentos fiquem
sempre positivos, ou seja: (1-((1-z)*x*(1/eta))) tem que ser sempre
positivo...

Por favor, caso vcs possam me ajudar, eu agradeço imensamente.

Att,

Romero.
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20141118/029ae1fe/attachment.html>