R中模拟广义帕累托分布时出现无效尺度参数错误求助

在R中研究广义帕累托分布（GPD），计算偏差、均方误差和均方根误差时，使用mle函数出现如下错误：

Error in dgpd(x = samples, scale = scale, shape = shape, log = TRUE) : 
  invalid scale

相关代码如下：

library(stats4)
library(evd)
gpd_simulation<-function(N,n,scale,shape){
  conc<-numeric(N)
  for(i in 1:N){
    samples=rgpd(n,scale=scale,shape=shape)
    gpdlik<-function(scale,shape){
      -sum(dgpd(x=samples,scale=scale,shape=shape,log=TRUE))
    }
    result<-mle(minuslogl=gpdlik,start=list(scale=11.63599901,shape=0.1),lower=c(0,0),method="L-BFGS-B")
    conc[i]<-result@coef
  }
  conc
}

summarygpd1<-function(conc){
  bias=mean(conc)-3;bias
  mse=var(conc)+bias^2;mse
  rmse=sqrt(var(conc)+bias^2);rmse
  return(c(bias,mse,rmse))
}

simulate40<-gpd_simulation(1000,40,11.63599901,0.1);simulate40
summarygpd1(simulate40)

现象：

• 使用自有数据集时无此错误
• 样本量增大到80、100、120时代码正常运行，样本量为40、60时错误重现

问题原因与解决方法 #

原因分析 #

错误和rgpd无关，核心问题是小样本下L-BFGS-B优化器容易试探到scale≤0的取值，而dgpd要求scale必须严格大于0。小样本时GPD似然函数形态更复杂，优化过程易触碰边界值，触发无效参数错误。另外代码存在两处逻辑错误：

1. conc[i]<-result@coef仅提取了第一个参数（scale），但result@coef是包含scale和shape的向量，会导致赋值异常
2. summarygpd1中偏差计算用mean(conc)-3，但真实scale是11.63599901，此处取值错误

修改后的代码 #

library(stats4)
library(evd)

gpd_simulation<-function(N,n,true_scale,true_shape){
  # 分别存储scale和shape的估计值
  scale_est <- numeric(N)
  shape_est <- numeric(N)
  
  for(i in 1:N){
    samples <- rgpd(n, scale = true_scale, shape = true_shape)
    
    gpdlik<-function(scale, shape){
      # 提前过滤无效scale值，强制优化器避开
      if(scale <= 0) return(Inf)
      -sum(dgpd(x = samples, scale = scale, shape = shape, log = TRUE))
    }
    
    # 用样本均值作为scale初始值，小样本下比固定值更稳定
    start_scale <- mean(samples)
    result <- mle(minuslogl = gpdlik,
                  start = list(scale = start_scale, shape = true_shape),
                  lower = c(1e-6, 0),  # 给scale设极小下界，避免试探0值
                  method = "L-BFGS-B")
    
    scale_est[i] <- result@coef[["scale"]]
    shape_est[i] <- result@coef[["shape"]]
  }
  # 返回两个参数的估计结果
  data.frame(scale_est, shape_est)
}

summarygpd1<-function(estimates, true_value){
  bias <- mean(estimates) - true_value
  mse <- var(estimates) + bias^2
  rmse <- sqrt(mse)
  return(c(bias = bias, mse = mse, rmse = rmse))
}

# 运行模拟
simulate40 <- gpd_simulation(1000, 40, 11.63599901, 0.1)
# 计算scale的偏差、MSE、RMSE
summarygpd1(simulate40$scale_est, 11.63599901)
# 计算shape的指标
summarygpd1(simulate40$shape_est, 0.1)

关键修改点 #

• 在似然函数中加入if(scale <=0) return(Inf)，强制优化器避开无效scale取值
• 将scale的下界从0改为1e-6，避免优化器试探到0值
• 用样本均值作为scale初始值，小样本下比固定初始值更适配
• 修正参数提取和偏差计算的逻辑错误，同时分别存储scale和shape的估计结果

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