R optim（）{fExtremes}得到0 hessian矩阵

-1

我正在使用R {fExtremes}为我的数据（一个向量）找到GEV分布的最佳参数。但出现以下错误消息R optim（）{fExtremes}得到0 hessian矩阵

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

我追溯，以适应$麻袋，发现我的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

我更新了我的数据集上的谷歌文档： https://docs.google.com/spreadsheets/d/1IRRpjmdrrJPhNmfiLism_P0efV_Ot4HlEsa6kwMnljc/edit?usp=sharing

来源

2016-08-23 D.deng

很难说出了什么问题，因为你的例子不可重现：请看这里。 http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example 由于“样本”数据集的特定特征（例如重复值），可能会导致失败。如果您不发布至少一部分数据，或者可能出现问题出现的其他数据集，我不会做太多工作。 –

嗨，我在这里更新了我的数据！ –

这将是一个很长的故事，可能更适合https://stats.stackexchange.com/。

======第1部分 - 的问题======

这是产生误差的序列：

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

我们面临lack-的典型原因在使用optim()和朋友时收敛性不足：优化的起始值不足。要看看出了什么问题，让我们使用PWM估计器（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）：xi的PWM估计值-0.46非常接近。现在，PWM估计是由gevFit()内部用作起始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。

来源

2016-08-29 14:37:37

非常感谢！解释很清楚，这真的有帮助！ –

R optim（）{fExtremes}得到0 hessian矩阵

回答

相关问题