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

[Romero Luiz M. Sales Filho]

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/20141117/e151f788/attachment.html>