3

Estou interessado em estimar dois parâmetros de máxima verossimilhança usando a função nlm.

**-------- FUNÇÃO RBETAMIFI -------**

rbetamifi <- function(n,mi,fi){  
p <- mi*fi  
q <- (1-mi)*fi  
return(rbeta(n,p,q))  
}

**----------------------------------**

aqui estão meus dados:

amostra10 <- matrix(NA,3,10)  
amostra10[1,] <- rbetamifi(10,0.4,13)  
amostra10[2,] <- rbetamifi(10,0.4,13)  
amostra10[3,] <- rbetamifi(10,0.4,13)  

Estimando mi e phi de cada amostra pelo métodos dos momentos (que será usado para o chute inicial):

mi1 <- mean(amostra10[1,])  
mi2 <- mean(amostra10[2,])  
mi3 <- mean(amostra10[3,])   

fi1 <- (mi1*(1-mi1)/var(amostra10[1,]))-1  
fi2 <- (mi2*(1-mi2)/var(amostra10[2,]))-1  
fi3 <- (mi3*(1-mi3)/var(amostra10[3,]))-1

Função de Verossimilhança:

menoslogV1 <- function(par){  
          logV <- length(amostra10[1,])*log(gamma(par[2])/(gamma(par[1]*par[2])*gamma((1-par[1])*par[2])))+  
            ((par[1]*par[2])-1)*sum(log(amostra10[1,]))+(((1-par[1])*par[2])-1)*sum(log(1-amostra10[1,]))  
          return(-logV)  
        }  

menoslogV2 <- function(par){
          logV <- length(amostra10[2,])*log(gamma(par[2])/(gamma(par[1]*par[2])*gamma((1-par[1])*par[2])))+
            ((par[1]*par[2])-1)*sum(log(amostra10[2,]))+(((1-par[1])*par[2])-1)*sum(log(1-amostra10[2,]))
          return(-logV)
        }

menoslogV3 <- function(par){
          logV <- length(amostra10[3,])*log(gamma(par[2])/(gamma(par[1]*par[2])*gamma((1-par[1])*par[2])))+
            ((par[1]*par[2])-1)*sum(log(amostra10[3,]))+(((1-par[1])*par[2])-1)*sum(log(1-amostra10[3,]))
          return(-logV)
        }

Estimando:

es1 <- nlm(menoslogV1, par <- c(mi1,fi1), hessian=TRUE)  
es2 <- nlm(menoslogV2, par <- c(mi2,fi2), hessian=TRUE)  
es3 <- nlm(menoslogV3, par <- c(mi3,fi3), hessian=TRUE)  

O problema é que quando vou estimar usando o NLM algumas vezes acontece esse erro:

Warning messages:
1: In log(gamma(par[2])/(gamma(par[1] * par[2]) * gamma((1 - par[1]) * : NaNs produced
2: In nlm(menoslogV2, par <- c(mi3, fi3), hessian = TRUE) : NA/Inf substituido pelo máximo valor positivo

Ele estima, mas estou com medo de que esses warnings estejam influenciando, como resolver?

4
  • 1
    nlm não permite NA (Not Avaiable) e NaN (Not a Number), não sei se seus dados tem estes valores ou são gerados nas formulas e funções, princialmente pq tem log no meio. Mas vi que vc colocou amostra10 <- matrix(NA,3,10), tenta assim amostra10 <- matrix(numeric(0),3,10) Commented 27/02/2015 às 7:10
  • Não tinha Na, esse código com os NAS era so pra criar a matriz, logo embaixo eu preenchi ela com as amostras. Commented 27/02/2015 às 18:45
  • @GabrielIanhezPereira você pode colocar seu comentário como uma resposta Commented 4/03/2015 às 0:59
  • @CarlosCinelli feito, obg! Commented 5/03/2015 às 3:33

1 Resposta 1

1

Meus professores descobriram o erro, era um erro numérico quando ele calculava o log(gamma), foi resolvido substituindo por uma função só que é a lgamma (que calcula direto o log da gamma):

menoslogV1 <- function(par){ 
  logV <- length(amostra10[1,])*(lgamma(par[2]) - 
    lgamma(par[1]*par[2]) - lgamma((1-par[1])*par[2]))+                                                 
    ((par[1]*par[2])-1)*sum(log(amostra10[1,])) + 
    (((1 - par[1])*par[2])-1)*sum(log(1-am‌​ostra10[1,])) 
 return(-logV) 
} 

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .