R optim(){fExtremes} получает 0 гессенскую матрицу

Question

Я использую R {fExtremes}, чтобы найти лучшие параметры распределения GEV для моих данных (вектор). но получите следующее сообщение об ошибке

Ошибка в solve.default(fit$hessian): Подпрограмма Лапака dgesv: система точно единственная: U[1,1] = 0

Я проследил до соответствия $hessian, обнаружил, что моя матрица гессиана - это сигулярная матрица, все элементы равны 0. Исходный код ( https://github.com/cran/fExtremes/blob/master/R/GevFit.R) gevFit() показывает, что fit $ hessian рассчитывается с помощью optim(). Выходные параметры имеют точно такое же значение, как и исходные параметры. Мне интересно, какие могут быть проблемы с моими данными, которые вызывают эту проблему? Я скопировал свой код здесь

> min(sample);
[1] 5.240909

> max(sample)
[1] 175.8677

> length(sample)
[1] 6789

> mean(sample)
[1] 78.04107

>para<-gevFit(sample, type = "mle")
Error in solve.default(fit$hessian) : 
  Lapack routine dgesv: system is exactly singular: U[1,1] = 0

fit = optim(theta, .gumLLH, hessian = TRUE, ..., tmp = data)
> fit

   $par

xi   -0.3129225
mu   72.5542497 
beta  16.4450897 

$value
[1] 1e+06

$counts
function gradient 
       4       NA 

$convergence
[1] 0

$message
NULL



$hessian

     xi  mu beta

xi    0    0     0

mu    0    0      0

beta  0     0      0

Я обновил свой набор данных в Документах Google: https://docs.google.com/spreadsheets/d/1IRRpjmdrrJPhNmfiLism_P0efV_Ot4HlEsa6kwMnljc/edit?usp=sharing

Answer 1

Это будет длинная история, возможно, более подходящая для https://stats.stackexchange.com/.

====== Часть 1 - проблема ======

Это последовательность, генерирующая ошибку:

library(fExtremes)
samp <- read.csv("optimdata.csv")[ ,2]
## does not converge
para <- gevFit(samp, type = "mle")

Мы сталкиваемся с типичной причиной отсутствия конвергенции при использовании optim() и друзья: неадекватные начальные значения для оптимизации.

Чтобы увидеть, что идет не так, давайте воспользуемся оценщиком ШИМ ( http://arxiv.org/abs/1310.3222); он состоит из аналитической формулы, следовательно, он не влечет за собой проблем сходимости, так как не использует optim():

para <- gevFit(samp, type = "pwm")
fitpwm<- attr(para, "fit")
fitpwm$par.ests

Расчетный параметр хвоста xi отрицательный, соответствующий ограниченному верхнему хвосту; на самом деле подогнанное распределение отображает даже больше "ограниченности верхнего хвоста", чем данные выборки, как вы можете видеть из "выравнивания" графика квантиль-квантиль справа:

qqgevplot <- function(samp, params){
  probs <- seq(0.1,0.99,by=0.01)
  qqempir <- quantile(samp, probs)
  qqtheor <-  qgev(probs, xi=params["xi"], mu=params["mu"], beta=params["beta"])
  rang <- range(qqempir,qqtheor)
  plot(qqempir, qqtheor, xlim=rang, ylim=rang,
     xlab="empirical", ylab="theoretical",
     main="Quantile-quantile plot")
  abline(a=0,b=1, col=2)
}
qqgevplot(samp, fitpwm$par.ests)

За xi<0.5 оценка MLE не является регулярной ( http://arxiv.org/abs/1301.5611): значение -0,46, оцененное ШИМ для xi очень близко к этому. Теперь оценки ШИМ используются внутренне gevFit() в качестве начальных значений для optim(): вы можете увидеть это, если распечатаете код для функции gevFit():

print(gevFit)
print(.gevFit)
print(.gevmleFit)

Начальное значение для optim theta, полученный ШИМ. Для конкретных данных это начальное значение не является адекватным, так как оно приводит к не сходимости optim(),

====== Часть 2 - решения? ======

Решение 1 заключается в использовании para <- gevFit(samp, type = "pwm") как указано выше. Если вы хотите использовать ML, вы должны указать хорошие начальные значения для optim(), К сожалению, fExtremes Пакет не позволяет сделать это легко. Затем вы можете переопределить свою собственную версию .gevmleFit включить те, например,

.gevmleFit <- function (data, block = NA, start.param, ...) 
{
  data = as.numeric(data)
  n = length(data)
  if(missing(start.param)){
    theta = .gevpwmFit(data)$par.ests
  }else{
    theta = start.param
  }
  fit = optim(theta, .gevLLH, hessian = TRUE, ..., tmp = data)
  if (fit$convergence) 
    warning("optimization may not have succeeded")
  par.ests = fit$par
  varcov = solve(fit$hessian)
  par.ses = sqrt(diag(varcov))
  ans = list(n = n, data = data, par.ests = par.ests, par.ses = par.ses, 
             varcov = varcov, converged = fit$convergence, nllh.final = fit$value)
  class(ans) = "gev"
  ans
}
## diverges, just as above
.gevmleFit(samp)
## diverges, just as above
startp <- fitpwm$par.ests
.gevmleFit(samp, start.param=startp)
## converges
startp <- structure(c(-0.1, 1, 1), names=names(fitpwm$par.ests))
.gevmleFit(samp, start.param=startp)$par.ests

Теперь проверьте это: beta по оценкам ШИМ 0,1245; изменяя это крошечное количество, MLE сделан, чтобы сходиться:

startp <- fitpwm$par.ests
startp["beta"]
startp["beta"] <- 0.13
.gevmleFit(samp, start.param=startp)$par.ests

Надеюсь, это ясно показывает, что слепо optim()Ise работает, пока не сработает, и может превратиться в довольно деликатное дело. По этой причине может быть полезно оставить этот ответ здесь, а не мигрировать в CrossValidated.