我编写了一个函数来计算广义帕累托分布的 MLE 估计值。当我将它与任何数据一起使用时,虽然我遇到了这样的错误
1: In log(beta * ksi) : NaNs produced
2: In nlm(loglik, theta, stepmax = 5000, iterlim = 1000) :
NA/Inf replaced by maximum positive value
我想知道是否有人可以发现我的代码有任何错误?
MLGPD<-function(data){
xi0 <- 1
beta0 <- 360
theta <- c(xi0, beta0)
excess <- data
assign("tmp", excess)
loglik <- function(theta){
ksi <- theta[1]
beta <- theta[2]
y <- ((tmp - 0.1)/beta)
f <- ((1/ksi)+1)*sum(log(1+y)) + length(tmp) * log(beta*ksi)
f
}
fit <- nlm(loglik, theta, stepmax = 5000, iterlim= 1000)
return(fit)
par.ests <- fit$x
return(par.ests)
}
#Checking our MLE algorithm works:
rgpd<-function(n,ksi, beta){
10000+beta*(((1-runif(n, min=0, max=1))^-ksi)-1)
}
rgpd1 <- rgpd(100, 1, 2.5)
MLGPD(rgpd1)
谢谢!