我正在研究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错误将再次发生。
1条答案
按热度按时间vmpqdwk31#
当我运行你的代码时，当
mle()函数尝试一个恰好为零的刻度时，它会失败。它从不选择一个接近零的MLE，所以最简单的方法可能是将参数的下限设置为一个小的正数，而不是使用零，例如，将下限设置为0.01对我来说是有效的，有警告但没有错误。警告来自
conc[i] <- result@coef行。此时result@coef的长度为2，conc[i]只能保存一个数字。如果您只想保存比例，请使用conc[i] <- result@coef[1]。此代码包含两个修复程序：
字符串
一个更一般的解决方案是计算出当scale趋于零时极限中应该发生什么，并让
gpdlik在scale正好为零时返回该值。