EM算法报错：含积分的似然估计中出现非有限函数值

我在R中实现基于蒙特卡洛模拟的极大似然估计（MLE），从修正威布尔分布抽样，用EM算法估计参数。似然函数包含积分，所有组件单独测试都正常，但运行时始终报以下错误：

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

完整代码 #

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

单独测试积分可行的代码 #

当我用真实参数单独计算这些积分时，能得到正常结果：

# 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

求帮忙定位问题或提供替代解决方法，感谢！

内容的提问来源于stack exchange，提问作者Qasim Ramzan