R语言log()生成NaN错误及蒙特卡洛模拟无结果问题求助

Hey there! Let's work through these issues step by step to get your BXIIWG Monte Carlo simulation running smoothly.

1. Fixing NaN Errors in Logarithms #

The NaN warnings you're seeing come from invalid inputs to the log() function—either values ≤ 0 being passed in. Here's why this happens and how to fix it:

Root Causes: #

• log(1 - theta) NaNs: When the estimated theta approaches or hits 1, 1 - theta gets very close to 0 (or even negative due to numerical precision), causing log() to return NaN or negative infinity.
• log(1 - theta*(1+x^rho)^(-lambda)*exp(-beta*x^alpha)) NaNs: If theta is near 1 and the product (1+x^rho)^(-lambda)*exp(-beta*x^alpha) is close to 1, their product can exceed 1, making 1 - ... negative.

Solutions: #

Adjust Parameter Bounds

Modify the upper bound for theta to a value slightly less than 1 (e.g., 0.999) to avoid boundary issues, and set non-zero lower bounds for all parameters to prevent zero-value inputs:

lower=c(alpha=1e-5,beta=1e-5,rho=1e-5,lambda=1e-5,theta=1e-5), 
upper=c(alpha=Inf,beta=Inf,rho=Inf,lambda=Inf,theta=0.999)

Add Safeguards to the Log-Likelihood Function

Include checks to ensure only valid values are passed to log(), returning infinity (a large penalty) if invalid parameters are tested:

BXIIWG_neglogl=function(alpha,beta,rho,lambda,theta){
  n = length(x)
  
  # Reject invalid parameter ranges immediately
  if (alpha <= 0 || beta <= 0 || rho <=0 || lambda <=0 || theta <=0 || theta >=1) {
    return(Inf)
  }
  
  term1 = alpha*beta*x^(alpha-1)*(1+x^rho) + rho*lambda*x^(rho-1)
  # Avoid log(0) or negative values
  if (any(term1 <= 0)) return(Inf)
  log_term1 = sum(log(term1))
  
  log_term2 = (lambda+1)*sum(log(1+x^rho))
  
  log_term3 = n*log(1-theta)
  
  sum_term4 = beta*sum(x^alpha)
  
  inner_term = theta*(1+x^rho)^(-lambda)*exp(-beta*x^alpha)
  # Ensure 1 - inner_term is positive
  if (any(inner_term >= 1)) return(Inf)
  log_term5 = 2*sum(log(1 - inner_term))
  
  # Return negative log-likelihood
  -log_term1 + log_term2 - log_term3 + sum_term4 + log_term5
}

2. Resolving Missing Mean/RMSE/Bias Outputs #

Your code has several indexing and assignment errors that prevent results from being calculated correctly. Here are the key fixes:

1. Correct Matrix Assignment #

You defined the results matrix coef1 but accidentally assigned values to an undefined coef matrix. Fix it to:

coef1[nsamp,] = coef(x1_BXIIWG)

2. Fix Bias Assignment Operator #

You used a subtraction (Bias-as.vector(...)) instead of an assignment (<-). Update it to:

Bias <- as.vector(sapply(1:5,function(x){
  Bias[(length(size)*(x-1)+1):(length(size)*x)] = Mean[(length(size)*(x-1)+1):(length(size)*x)] - par1[x]
}))

3. Simplify Index Logic with Matrices #

Using matrices to store intermediate results makes indexing clearer and avoids errors. Replace vector-based storage with matrices:

# Initialize result matrices
Mean_mat = matrix(NA, nrow=length(size), ncol=5, dimnames=list(size, c('alpha','beta','rho','lambda','theta')))
RMSE_mat = matrix(NA, nrow=length(size), ncol=5, dimnames=list(size, c('alpha','beta','rho','lambda','theta')))
Bias_mat = matrix(NA, nrow=length(size), ncol=5, dimnames=list(size, c('alpha','beta','rho','lambda','theta')))

# Inside the iter_size loop:
Mean_mat[iter_size,] = apply(coef1,2,mean,na.rm=TRUE)
RMSE_mat[iter_size,] = apply((coef1 - matrix(rep(par1, samp), ncol=5, byrow=TRUE))^2, 2, function(x){sqrt(mean(x,na.rm=TRUE))})
Bias_mat[iter_size,] = Mean_mat[iter_size,] - par1

# Flatten into vectors for final output
Mean = as.vector(t(Mean_mat))
RMSE = as.vector(t(RMSE_mat))
Bias = as.vector(t(Bias_mat))
samplesize = as.vector(t(mapply(rep, size, 5)))

4. Define n in the Log-Likelihood #

Your BXIIWG_neglogl function uses n but never defines it. Add this line at the start of the function:

n = length(x)

5. Correct Function Call #

When invoking BXIIWG_simulation, explicitly specify the size parameter to avoid mismatched arguments:

BXIIWGsim1 <- BXIIWG_simulation(size=c(25), samp=1000, par1=par1)

3. Full Corrected Code Example #

Here’s the complete code with all fixes applied:

library(rootSolve)
library(Matrix)
library(bbmle)

# True parameters
alpha=4
beta=0.3
rho=2
lambda=0.5
theta=0.2
samp=1000
par1=c(alpha,beta,rho,lambda,theta)

####Define BXIIWG quantile
BXIIWG_quantile=function(alpha,beta,rho,lambda,theta,u){
  f=function(x){
    beta*x^alpha + lambda*log(1+x^rho) + log(1-u)
  }
  x=uniroot(f,c(0,100),tol=0.0001)$root
  return(x)
}

####Define BXIIWG log-likelihood with safeguards
BXIIWG_neglogl=function(alpha,beta,rho,lambda,theta){
  n = length(x)
  
  # Reject invalid parameter ranges immediately
  if (alpha <= 0 || beta <= 0 || rho <=0 || lambda <=0 || theta <=0 || theta >=1) {
    return(Inf)
  }
  
  term1 = alpha*beta*x^(alpha-1)*(1+x^rho) + rho*lambda*x^(rho-1)
  # Avoid log(0) or negative values
  if (any(term1 <= 0)) return(Inf)
  log_term1 = sum(log(term1))
  
  log_term2 = (lambda+1)*sum(log(1+x^rho))
  
  log_term3 = n*log(1-theta)
  
  sum_term4 = beta*sum(x^alpha)
  
  inner_term = theta*(1+x^rho)^(-lambda)*exp(-beta*x^alpha)
  # Ensure 1 - inner_term is positive
  if (any(inner_term >= 1)) return(Inf)
  log_term5 = 2*sum(log(1 - inner_term))
  
  # Return negative log-likelihood
  -log_term1 + log_term2 - log_term3 + sum_term4 + log_term5
}

####Define simulation process of BXIIWG
BXIIWG_simulation=function(size=c(25,50,100,200,400,800),samp,par1){
  # Initialize result matrices
  Mean_mat = matrix(NA, nrow=length(size), ncol=5, dimnames=list(size, c('alpha','beta','rho','lambda','theta')))
  RMSE_mat = matrix(NA, nrow=length(size), ncol=5, dimnames=list(size, c('alpha','beta','rho','lambda','theta')))
  Bias_mat = matrix(NA, nrow=length(size), ncol=5, dimnames=list(size, c('alpha','beta','rho','lambda','theta')))
  
  for (iter_size in 1:length(size)){
    current_n = size[iter_size]
    coef1=matrix(NA,samp,5)
    colnames(coef1)=c('alpha','beta','rho','lambda','theta')
    
    for (nsamp in 1:samp){
      tryCatch(
        {
          # Generate sample
          q=runif(current_n,0,1)
          x1=sapply(q,BXIIWG_quantile, alpha=par1[1],beta=par1[2],rho=par1[3],lambda=par1[4],theta=par1[5])
          
          # Fit model
          x1_BXIIWG<-mle2(BXIIWG_neglogl, 
                          start=list(alpha=par1[1],beta=par1[2],rho=par1[3],lambda=par1[4],theta=par1[5]), 
                          method="L-BFGS-B",
                          data=list(x=x1), 
                          lower=c(alpha=1e-5,beta=1e-5,rho=1e-5,lambda=1e-5,theta=1e-5), 
                          upper=c(alpha=Inf,beta=Inf,rho=Inf,lambda=Inf,theta=0.999),
                          use.ginv=TRUE)
          coef1[nsamp,]=coef(x1_BXIIWG)
        },error=function(e){}
      )
    }
    
    # Calculate metrics
    Mean_mat[iter_size,] = apply(coef1,2,mean,na.rm=TRUE)
    RMSE_mat[iter_size,] = apply((coef1 - matrix(rep(par1, samp), ncol=5, byrow=TRUE))^2, 2, function(x){sqrt(mean(x,na.rm=TRUE))})
    Bias_mat[iter_size,] = Mean_mat[iter_size,] - par1
  }
  
  # Flatten into vectors for final output
  Mean = as.vector(t(Mean_mat))
  RMSE = as.vector(t(RMSE_mat))
  Bias = as.vector(t(Bias_mat))
  samplesize = as.vector(t(mapply(rep, size, 5)))
  
  return(data.frame(samplesize,Mean,RMSE,Bias))
}

# Run simulation
BXIIWGsim1<-BXIIWG_simulation(size=c(25), samp=1000, par1=par1)
Mean1<-BXIIWGsim1$Mean
RMSE1<-BXIIWGsim1$RMSE
Bias1<-BXIIWGsim1$Bias

# View results
print(BXIIWGsim1)

Additional Tips #

• Test with Smaller Samples: For faster debugging, set samp=100 instead of 1000 to verify the code works before running full-scale simulations.
• Debug Failed Fits: Add a print statement inside the error function of tryCatch to see why some iterations fail (e.g., error=function(e){print(paste("Iteration", nsamp, "failed:", e))}).

内容的提问来源于stack exchange，提问作者S.A. Osagie