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.