模拟广义帕累托分布以查找 mse rmse 和偏差时产生无效量表

Question

我正在研究 R 中的广义帕累托分布，我想计算分布的偏差、mse 和 rmse。但是，当我使用 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)

但是，如果我使用自己的数据集，就不会出现以下错误，这样使用rgpd有什么问题吗？

补充一下，如果我增加样本量，代码运行得很好，我将其增加到 80、100 和 120，代码运行没有错误，但 40 和 60 会再次出现错误。

Answer 1

当我运行您的代码时，当

mle()

函数尝试恰好为零的比例时，它会失败。它永远不会选择接近零的 MLE，因此将参数的下限设置为一个小的正数而不是使用零可能是最简单的，例如将下限设置为 0.01 对我来说有效，但有警告但没有错误。

警告来自

conc[i] <- result@coef

行。此时

result@coef

的长度为2，并且

conc[i]

只能容纳单个数字。如果您只想保存比例，请使用

conc[i] <- result@coef[1]

。

此代码包含两个修复：

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.01,0),method="L-BFGS-B")
    conc[i]<-result@coef[1]
  }
  conc
}

一个更通用的解决方案是弄清楚当标度趋于零时在极限内应该发生什么，并在

gpdlik

恰好为零时让

scale

返回该值。