Error in EM algorithm: Non-finite Function Value in Likelihood Estimation with Integrals

Question

I am encountering a persistent issue in my R program while performing Monte Carlo simulations for Maximum Likelihood Estimation (MLE). The simulation involves drawing samples from a modified Weibull distribution and estimating parameters through the EM algorithm. The likelihood function includes integrals, and despite thorough testing of individual components, I consistently face the error message

Error in integrate(ll1, x1e[j], Inf, rel.tol = 1e-06) :non-finite function value

Here is my program:

Code:

###### R Codes for Monte Carlo Simulation####
##########################
##########################
##########################
# Install and load required libraries
#
library(pracma)
library('bbmle') 
library(numDeriv)

##################################################################
N= 15;n=300; m=100;tau=.4;tau_2=3;lambda=1.5;phi=0.5;paii=0.2;gamma=3;delta=.5; initial_values1 <- c(1.2, 1.6, 2, 1.5);p=.4; k1 <- c(-.5,.5)  
 
inverse_transform_sampling <- function(n, pdf, lower_bound, upper_bound) {
  # Generate uniform random numbers
  u <- runif(n);    x <- numeric(n)
  for (t in 1:n) {
    inverse_cdf <- function(x) integrate(pdf, lower_bound, x)$value - u[t]
    x[t] <- uniroot(inverse_cdf, lower = lower_bound, upper = upper_bound)$root
  }
  return(x)
}
pdf_modified_weibull <- function(x) {
  phi * x^(gamma - 1) * (gamma + delta * x) * exp(-phi * x^gamma * exp(delta * x) + delta * x)
}
T   <- inverse_transform_sampling(n, pdf_modified_weibull, 0, 1000)
###############################################

drawn_observations <- r <-  numeric() ###vector("list", length = m)
removed_observations <- vector("list", length = m)
# Perform the iterations
for (j in 1:m) {
  
  observation <- min(T) #sample(y[y > max(c(prev_observation, -Inf)) ], 1)
  drawn_observations[j] <- observation
  
  # Remove r observations at each iteration (except the last one)
  if (j < m) {
    rr <- rbinom(1, 5 , p)  # Draw the number of observations to remove
    r[j] <- rr
    removed_obs <- sample(setdiff(T, observation), rr, replace = FALSE)
    removed_observations[[j]] <- removed_obs
    rm<- c(observation, removed_obs)
    T <- T[!(T %in% rm)]# setdiff(y, removed_obs);y[y!=observation]
  } else {  # On the last iteration, remove all remaining observations
    removed_observations[[j]] <- T
    r[j] <- length(T)
    T <- NULL  # Clear y as all observations have been drawn
  }
  
  # Update the previously drawn observation
  prev_observation <- observation
}

x1 <- drawn_observations[drawn_observations <= tau ]
xm <- drawn_observations[!(drawn_observations %in% x1)]
m <- length(r); n1 <- length(x1); m1 <- length(xm)
rn1 <- r[1:n1]; rm <- r[ (n1+1) : m]; yc <- max(xm)    
##################################################################
rn1e <- rn1#[rn1 != 0]
rme <- rm#[rm != 0]
x1e <- x1#[rn1 != 0]
xme <- xm#[rm != 0]

  # EM Algorithm
    # Define the functions for the log-likelihood and its gradient
    loglik <- function(params) {
      # Extract parameters
      gamm <- params[1]
      delt <- params[2]
      ph <- params[3]
      lambd <- params[4]

    rn1e <- rn1#[rn1 != 0]
    rme <- rm#[rm != 0]
    x1e <- x1#[rn1 != 0]
    xme <- xm#[rm != 0]
    # Calculate intermediate values length(llik2)   ph; lambd=lamd

      gdx1 <- gamm + delt * x1e
      tplxm <- tau + lambd * (xme -tau)
      gdtplxm <- gamm + delt * (tplxm)
      edtplxm <- exp(delt * (tplxm))

      # Calculate the log-likelihood
      loglik <- m * log(ph) + (m - n1) * log(lambd) + sum((gamm - 1) * log(x1e) + log(gdx1)) -
        ph * sum(x1e ^ gamm * exp(delt * x1e)) + delt * sum(x1e) +
        sum((gamm - 1) * log(tplxm) + log(gdtplxm)) -
        ph * sum(tplxm ^ gamm * edtplxm) + delt * sum(tplxm)  

   llik1 <- c(); llik2 <- c();
      for (j in 1:length(x1e)) {

        ll1 <- function(z) {
          lll11 <- (2 * gamm - 1)* log(z) +log(gamm + delt * z) 
          lll12 <- (-ph * z^gamm * exp(delt * z) + 2 * delt * z)
          return(ifelse(z==0, 0, exp(lll11+lll12)))
        }

        llik1[j] <- (ph * rn1e[j] * integrate(ll1, x1e[j], 1.3, rel.tol = 1e-6)$value) / exp(-ph * x1e[j] ^ gamm * exp(delt * x1e[j]))
      }
  for (j in 1:length(xme) ) {  # ; j=12
          ll2<- function(z){
             lll21 <-   (2 * gamm - 1)*log(tau + lambd * (z - tau)) + log(gamm + delt * (tau + lambd * (z - tau))) 
          lll22 <-   (-ph * (tau + lambd * (z - tau)) ^ gamm * exp(delt * (tau + lambd * (z - tau))) + 2 * delt * (tau + lambd * (z - tau)))
            return(ifelse(z==0, 0, exp(lll21+lll22)))
          }

        #    j=1
        llik2[j] <-  (ph * rme[j] * integrate(ll2, xme[j], 30, rel.tol = 1e-6)$value) / exp(-ph * tplxm[j] ^ gamm * exp(delt * tplxm[j]))   
           }
      logli <- loglik - sum(llik1)- sum(llik2)
      return(logli)
}
    # Estimation under EM Algorithm using optim
    opt_result <- nlminb(c(gamma, delta, phi,lambda), loglik  , lower = 0, upper = Inf)
# Estimation under EM Algorithm using optim
print(opt_result) 

It's worth noting that when I calculate these integrals separately, they yield results as follows:

    # Calculate intermediate values
      gdx1 <- gamma + delta * x1e
      tplxm <- tau + lambda * (xme - tau)
      gdtplxm <- gamma + delta * (tplxm)
      edtplxm <- exp(delta * (tplxm))
      
      # Calculate the log-likelihood
      loglik <- m * log(phi) + (m - n1) * log(lambda) + sum((gamma - 1) * log(x1e) + log(gdx1)) -
        phi * sum(x1e ^ gamma * exp(delta * x1e)) + delta * sum(x1e) +
        sum((gamma - 1) * log(tplxm) + log(gdtplxm)) -
        phi * sum(tplxm ^ gamma * edtplxm) + delta * sum(tplxm)  

            llik1 <- c(); llik2 <- c(); 
      for (j in 1:length(x1e)) {
        
        ll1 <- function(z) {
          
          lll11 <- (2 * gamma - 1)* log(z) +log(gamma + delta * z) 
          lll12 <- (-phi * z^gamma * exp(delta * z) + 2 * delta * z)
          return(ifelse(z==0, 0, exp(lll11+lll12)))
        }
         
        llik1[j] <- (phi * rn1e[j] * integrate(ll1, x1e[j], 3, rel.tol = 1e-6)$value) / exp(-phi * x1e[j] ^ gamma * exp(delta * x1e[j]))
        
      }
      
      for (j in 1:length(xme) ) {  # ; j=12
       
          ll2<- function(z){
             lll21 <-   (2 * gamma - 1)*log(tau + lambda * (z - tau)) + log(gamma + delta * (tau + lambda * (z - tau))) 
            lll22 <-   (-phi * (tau + lambda * (z - tau)) ^ gamma * exp(delta * (tau + lambda * (z - tau))) + 2 * delta * (tau + lambda * (z - tau)))
            return(ifelse(z==0, 0, exp(lll21+lll22)))
          }
        #    j=1
        llik2[j] <-  (phi * rme[j] * integrate(ll2, xme[j], 3, rel.tol = 1e-6)$value) / exp(-phi * tplxm[j] ^ gamma * exp(delta * tplxm[j]))   
           }
     

logli <- loglik - sum(llik1)- sum(llik2)
logli

Can someone help me identify the issue or suggest an alternative approach?

Thank you in advance for your assistance!

Answer 1

If you set the option error=utils::recover you can inspect the state of the object at the time of an error inside a function.

options(error = recover)

opt_result <- nlminb(c(gamma, delta, phi,lambda), loglik  , lower = 0, upper = Inf)
------ result ----------
Error in integrate(ll1, x1e[j], 1.3, rel.tol = 1e-06) : 
  non-finite function value

Enter a frame number, or 0 to exit   

1: nlminb(c(gamma, delta, phi, lambda), loglik, lower = 0, upper = Inf)
2: objective(.par, ...)
3: #34: integrate(ll1, x1e[j], 1.3, rel.tol = 1e-06)

Selection: 3    #choosing to look at the state at time of error at innermost location
Called from: objective(.par, ...)
Browse[1]> j
[1] 83

Browse[1]> x1e
[1] 0.2307633 0.2406504 0.2893463 0.2944903 0.3232987 0.3649652 0.3796235 0.3845039

#---------
# when it's time to exit debug-mode, exit the debugger with "0"'s and then:

options(error = NULL)

Now it's your job to figure out why you are trying to use the 83rd item in an 8-element vector. You might also pick up the habit of using set.seed() when constructing simulated random draws. I can't tell if that is the issue here.