Maximum likelihood estimation using nlm

Asked

Viewed 333 times

3

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

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

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

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

here is 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)

Cherishing mi and phi of each sample by the timing methods (which will be used for the initial kick):

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 estimate using NLM sometimes this error happens:

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 replaced by maximum positive value

He does, but I’m afraid these warnings are influencing, how to solve?

r
  • 1

    nlm does not allow NA (Not Avaiable) and Nan (Not a Number), do not know if your data has these values or are generated in the formulas and functions, but pq has log in between. But I saw that you put sample10 <- Matrix(NA,3,10), so try sample10 <- Matrix(Numeric(0),3,10)

  • There was no Na, this code with the NAS was just to create the matrix, right below I filled it with the samples.

  • @Gabrielianhezpereira you can put your comment as a response

  • @Carloscinelli done, obg!

1 answer

1


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

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) 
}

Browser other questions tagged r

You are not signed in. Login or sign up in order to post.