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.