Неверный гессиан из Оптима в R

Question

Я делаю анализ экстремальных значений. Я не хочу использовать пакет fevd по разным причинам (во-первых, я хочу иметь возможность настроить некоторые вещи, которые я не могу сделать иначе). Я написал свой собственный код. В основном это очень просто, и я подумал, что все решил. Но для некоторых комбинаций параметров Hessian, полученный из моего логарифмического анализа правдоподобия (основанного на optim), не будет правильным.

Переходя на один шаг в то время. Мой код - или выбранная его часть - выглядит так:

# 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, я получаю точно такие же результаты. В некоторых случаях (в моем конкретном случае, когда форма =-0,29 так сильно отрицательна / вейбулла), моя рутина будет давать отрицательные отклонения и фанки-гессианы. Это не всегда неправильно, но некоторые комбинации параметров явно не дают действительного гессиана (Примечание: параметры по-прежнему оцениваются правильно, то есть идентичны результатам fevd, но ковариационная матрица полностью отключена).

Я нашел этот пост, в котором сравнивались гессианы из двух процедур, и, действительно, оптимизм кажется странным. Однако, если я просто заменю maxLik в своей процедуре, он просто не сходится вообще (даже в тех случаях, когда происходила конвергенция).

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

Я пытался дать разные начальные значения - даже правильные - но они просто не сходятся.

Я не предоставляю данные, потому что я думаю, что в моем примере правильно используется подпрограмма optim. Просто численные результаты не являются стабильными для некоторой комбинации параметров. Мой вопрос:

1) Я что-то упускаю при использовании maxLik?

2) Есть ли другие процедуры оптимизации, кроме maxLik, из которых я могу извлечь гессиан?

Спасибо