r – Maximum likelihood estimation using nlm

Question:

I am interested in estimating two maximum likelihood parameters using the nlm function.

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

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

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

here are my data:

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)  

Estimating mi and phi of each sample by the method of moments (which will be used for the kickoff):

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

Likelihood Function:

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)  

The problem is that when I go to estimate using NLM this error sometimes happens:

Warning messages:
1: In log(gamma(par[2])/(gamma(par[1] * par[2]) * gamma((1 – par[1]) * : NaNs produced
2: In nlm(minuslogV2, par <- c(mi3, fi3), hessian = TRUE) : ​​NA/Inf replaced by the maximum positive value

He appreciates it, but I'm afraid that these warnings are influencing, how to solve it?

Answer:

My teachers discovered the error, it was a numerical error when he calculated the log(gamma), it was solved by replacing it with a function only that is the lgamma (which calculates the log of the gamma directly):

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) 
} 
Scroll to Top