1

我正在做一些极值分析。出于各种原因，我不想使用 fevd 包（首先我希望能够调整一些我无法做到的事情）。我编写了自己的代码。它大多非常简单，我以为我已经解决了所有问题。但是对于某些参数组合，我的对数似然分析（基于optim）得出的 Hessian 将不正确。

一步一步走过去。我的代码 - 或它的选定部分 - 看起来像这样：

# routines for non stationary
Log_lik_GEV <- function(dataIN,scaleIN,shapeIN,locationIN){
    # simply calculate the negative log likelihood value for a set of X and parameters, for the GPD
    #xi, mu, sigma  - xi is the shape parameter, mu the location parameter, and sigma is the scale parameter.
    # shape = xi
    # location = mu
    # scale = beta
    library(fExtremes)
    #dgev   Density of the GEV Distribution, dgev(x, xi = 1, mu = 0, sigma = 1)

    LLvalues <- dgev(dataIN, xi = shapeIN, mu = locationIN, beta = scaleIN) 
    NLL <- -sum(log(LLvalues[is.finite(LLvalues)]))
    return(NLL)
}


function_MLE <- function(par , dataIN){
    scoreLL <- 0
    shape_param <- par[1]
    scale_param <- par[2]
    location_param <- par[3]
    scoreLL <- Log_lik_GEV(dataIN, scale_param, shape_param, location_param)
    if (abs(shape_param) > 0.3) scoreLL <- scoreLL*10000000
    if ((scale_param) <= 0) {
        scale_param <- abs(scale_param)
        par[2] <- abs(scale_param)
        scoreLL <- scoreLL*1000000000
    }
    sum(scoreLL) 
}

kernel_estimation <- function(dati_AM, shape_o, scale_o, location_o) {

     paramOUT <- optim(par = c(shape_o, scale_o, location_o), fn = function_MLE, dataIN = dati_AM, control = list(maxit = 3000, reltol = 0.00000001), hessian = TRUE)

     # calculation std errors
     covmat <- solve(paramOUT$hessian)
     stde <- sqrt(diag(covmat))
     print(covmat)

     print('')

     result <- list(shape_gev =paramOUT$par[1], scale_gev = paramOUT$par[2],location_gev =paramOUT$par[3], var_covar = covmat)

     return(result)
}

在某些情况下，一切都很好。如果我运行我的例程和 fevd 例程，我会得到完全相同的结果。在某些情况下（在我的具体情况下，当 shape=-0.29 如此强烈地为负/weibull 时），我的例程会给出负方差和时髦的粗麻布。它并不总是错误的，但某些参数组合显然没有给出有效的 hessian （注意：参数仍然估计正确，含义与 fevd 结果相同，但协方差矩阵完全关闭）。

我发现这篇文章比较了两个程序的粗麻布，确实 optim 似乎很不稳定。但是，如果我只是在我的例程中替换 maxLik，它根本不会收敛（即使在收敛发生的情况下）。

 paramOUT = maxLik(function_MLE, start =c(shape_o, scale_o, location_o),   
                          dataIN=dati_AM, method ='NR' )

我试图给出不同的初始值——即使是正确的值——但它就是不收敛。

我没有提供数据，因为我认为在我的示例中正确使用了 optim 例程。简而言之，对于某些参数组合，数值结果并不稳定。我的问题是：

1) 我在使用 maxLik 的方式上是否遗漏了什么？

2) 除了 maxLik 之外，还有其他优化例程可以从中提取粗麻布吗？

谢谢

4

0 回答 0