R optim(){fExtremes} 得到 0 hessian 矩阵

Question

我正在使用 R {fExtremes} 为我的数据（向量）寻找 GEV 分布的最佳参数。但得到以下错误信息

Error in solve.default(fit$hessian) : Lapack routine dgesv: system is exactly singular: U[1,1] = 0

我追溯到fit$hessian，发现我的hessian矩阵是一个奇异矩阵，所有元素都是0。 gevFit()的源代码(https://github.com/cran/fExtremes/blob/master/R/GevFit.R)显示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() 和朋友时面临缺乏收敛的典型原因：优化的起始值不足。

要查看问题出在哪里，让我们使用 PWM 估算器 (http://arxiv.org/abs/1310.3222)；这由一个分析公式组成，因此不会引起收敛问题，因为它没有使用 optim():

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

估计尾参数xi为负，对应有界上尾；实际上，拟合分布显示的 "upper tail boundedness" 比样本数据还要多，正如您从右侧分位数-分位数图的 "leveling off" 中看到的那样：

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）：PWM 为 xi 估计的 -0.46 值非常接近。现在，gevFit() 在内部使用 PWM 估计值作为 optim() 的起始值：如果您打印函数 gevFit():

的代码，您可以看到这一点

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

optim的起始值为theta，通过PWM获得。对于手头的具体数据，这个起始值是不够的，因为它导致 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

现在检查一下：PWM 估计的 beta 是 0.1245；通过将其更改为很小的量，使 MLE 收敛：

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

这很有希望清楚地说明，盲目地 optim()ise 工作直到它不工作，然后可能会变成一个非常微妙的努力。因此，将此回复留在此处可能比迁移到 CrossValidated 更有用。

R optim(){fExtremes} 得到 0 hessian 矩阵

R optim(){fExtremes} gets 0 hessian matrix

math

optimization

r

hessian

mle