Maximum likelihood estimation using nlm

Question

Maximum likelihood estimation using nlm

Navigation

#1 by (1 votes)

3

I'm 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 methods of the moments (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'm estimating using the NLM sometimes this error happens:

Warning messages:
  1: In log (gamma (par [2]) / (gamma (par [1] * par [2]) * gamma ((1 -     NaNs produced
  2: In nlm (lesslogV2, par

r

asked by anonymous 26.02.2015 / 23:25

1 answer

Route system in PHP using the Google Maps API Slow response PHP API

score 1 · Answer 1

My teachers discovered the error, it was a numerical error when he calculated the log (gamma), it was solved by substituting a function only that is the lgamma (which calculates the gamma log 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) 
}