r - 如何借助 lmomco 函数在 R 中定义自己的 fitdistr 函数分布

Question

我想定义自己的分布，与 fitdistrplus 函数一起使用，以适应从现在起称为“月”的每月降水数据。我正在使用“lmomco”函数来帮助我定义发行版，但无法使其工作。例如，我定义广义极值 (gev) 分布，如下所示:

dgev<-pdfgev   #functions which are included in lmomco
 pgev<-cdfgev
qgev<-quagev

由于“fitdistrplus”需要参数“start”，它由所需分布的初始参数值组成，因此我估计这些初始值如下:

lmom=lmoms(month,nmom=5)     #from lmomco package
para=pargev(lmom, checklmom=TRUE)

现在，我终于尝试使用“fitdist”函数将“month”拟合到gev分布中:

fitgev <- fitdist(month, "gev", start=para[2]) #fitdistrplus

我收到如下错误。无论我在“lmomco”的帮助下定义哪个发行版，我都会得到同样的错误。有人可以提示我我做错了什么吗？谢谢!

fitgev <- fitdist(month, "gev", start=para[2])
[1] "Error in dgev(c(27.6, 97.9, 100.6, 107.3, 108.5, 109, 112.4, 120.9, 137.8,  : \n  unused arguments (para.xi = 196.19347977195, para.alpha = 91.9579520442104, para.kappa = -0.00762962879097294)\n"
attr(,"class")
[1] "try-error"
attr(,"condition")
<simpleError in dgev(c(27.6, 97.9, 100.6, 107.3, 108.5, 109, 112.4, 120.9, 137.8, 138.4, 144.7, 156.8, 163.1, 168.9, 169.1, 171.4, 176.1, 177.1, 178.8, 178.9, 187.2, 190.2, 190.5, 190.8, 191.2, 193.1, 195.2, 198.5, 199.8, 201.7, 206.9, 213.4, 220.7, 240, 253.5, 254.5, 256.1, 256.4, 257.5, 258.3, 261.5, 263.7, 264.7, 279.1, 284.2, 313.1, 314.7, 319.4, 321.6, 328.9, 330.1, 332.2, 358.3, 366.8, 367.9, 403.5, 424.1, 425.9, 457.3, 459.7, 468, 497.1, 508.5, 547.1), para.xi = 196.19347977195, para.alpha = 91.9579520442104,     para.kappa = -0.00762962879097294): unused arguments (para.xi = 196.19347977195, para.alpha = 91.9579520442104, para.kappa = -0.00762962879097294)>
Error in fitdist(month, "gev", start = para[2]) : 
  the function mle failed to estimate the parameters, 
                with the error code 100

Answer 1

tl;dr 这是很挑剔的，而且可能永远都是很挑剔的——将潜在不稳定的分布拟合到极小、嘈杂的数据集，是非常困难的。我在下面概述了一些策略，这些策略将为我们提供答案，但我并不真正相信我得到的任何答案。

对于此处的具体情况，@BelSmek 的答案是最好的:evd::fgev(month) 给出与下面的 mle2/DEoptim 匹配的答案、和给出了更合理的标准误差估计。然而，下面的所有阴谋对于那些试图将参数拟合到一般分布的人来说可能都是有用的东西......

fitdist 需要一个带有命名参数的密度/分布函数，以及更多；我们可以做到这一点，尽管正如我所说，我不相信答案。

library("lmomco")
library("fitdistrplus")
## reproducible:
month <- c(27.6, 97.9, 100.6, 107.3, 108.5,
              109, 112.4, 120.9, 137.8)

设置:

lmom <- lmoms(month,nmom=5)     #from lmomco package
para <- pargev(lmom, checklmom=TRUE)

事实证明，我们需要重新定义 dgev，添加一些额外的管道，以使每个人都满意:

pgev <- function(q, xi, alpha, kappa) {
    if (length(q) == 0) return(numeric(0))
    r <- try(cdfgev(x = q, para = c(xi = xi, alpha = alpha, kappa = kappa)), 
           silent = TRUE)
    if (inherits(r, "try-error")) return(rep(NaN, length(q)))
    r
}
dgev <- function(x,xi,alpha,kappa, minval = 1e-8) {
    r <- pdfgev(x,list(type="gev",para=c(xi,alpha,kappa),source="pargev"))
    r[r==0] <- minval
    r
}

除了将参数从向量更改为列表之外，这里最重要的可能是拦截密度函数下溢到零的情况并将其替换为一个小值。这是一个并不总是有效的技巧:更原则的方法是让密度函数直接计算对数密度(我将在下面尝试这个，尽管在这种情况下它没有多大帮助)。

fitgev <- fitdist(month, "gev", start=as.list(para[[2]]))

我们得到了答案...

Parameters:
        estimate   Std. Error
xi    104.060486 0.0004131185
alpha  39.227041 0.0004150259
kappa   1.162644 0.0004105323

...但我根本不相信这一点，因为标准误差低得不切实际(为什么我们认为在将 3 参数模型拟合到 9 个数据点时可以如此精确地估计参数......？)

另一种方法将 bbmle::mle2 与 evd::dgev 结合使用 - 后者确实有一个 log 参数...

## clean up
rm(dgev)
detach("package:lmomco")
## get new packages
library(evd)
library(bbmle) 

(一般来说，最好在这里开始一个新的 R session ......)

我再次必须包装 dgev 函数来替换不可能的值(即使我们现在正在使用对数刻度，所以事情更加稳定......)

dgev <- function(..., log = FALSE, minval = 1e-8) {
    r <- evd::dgev(..., log = log)
    if (log) {
        r[r == -Inf] <- log(minval)
    }
    r
}
fit2 <- mle2(month ~ dgev(loc = xi, scale = alpha, shape = kappa), 
     data = data.frame(month),
     start = as.list(para[[2]]))
summary(fit2)

请注意，标准误差现在稍微更合理，但仍然小得惊人，而且这些答案与我们从 fitdistrplus< 得到的答案完全不同/.

Coefficients:
        Estimate Std. Error z value     Pr(z)    
xi    99.6720328  0.0765906 1301.36 < 2.2e-16 ***
alpha 30.7447099  0.3027090  101.57 < 2.2e-16 ***
kappa -0.7763013  0.0076273 -101.78 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-2 log L: 82.063 

作为最终的强力方法，我们将尝试差分进化

dgev_lik <- function(pars, minval = 1e-8) {
    r <- evd::dgev(month, pars[1], pars[2], pars[3], log = TRUE)
    r[r == -Inf] <- log(minval)
    -1*sum(r)
}

library(DEoptim)
set.seed(101)
d1 <- DEoptim(dgev_lik, lower = c(90, 10, -2),
        upper = c(130, 50, 2),
        control = DEoptim.control(NP = 1000, itermax = 1000))

d1$optim
$bestmem
      par1       par2       par3 
99.6299712 30.7704978 -0.7762563 

$bestval
[1] 41.03149

这与 mle2 得到的答案基本相同。看看 fitgev 的内部结构，它声称比 mle2 具有更好的对数似然性 (logLik(fitgev) 为 -36.9，而 mle2/DEoptim 为 -41)，但它似乎正在计算不可比较的值:插入 fitgev 将参数直接输入到我们的对数似然函数中会给出更更差的答案(对于负对数似然，值越高越差......)

dgev_lik(fitgev$estimate) ## 57.39