我想定义我自己的分布，以与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"函数将"月"拟合到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
tl;博士这很挑剔，而且可能总是挑剔 - 将潜在的不稳定分布拟合到极小、嘈杂的数据集中，简直很难。我在下面概述了一些策略，这些策略将得到我们的答案，但我并不真正相信我得到的任何答案。
对于这里的具体情况，@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 会话......
我再次不得不包装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
确保累积函数中的参数具有变量q: pgev(q, par1, par2)而不是pgev(x, par1, par2)
因为错误消息本质上告诉您它找不到变量 q。
关键点是使用：x作为 pdf 输入;q作为 CDF 输入;p作为反向 CDF 输入
例如，拟合您自己定义的 Gumble 分布
# Data
x1 <- c(6.4,13.3,4.1,1.3,14.1,10.6,9.9,9.6,15.3,22.1,13.4,
13.2,8.4,6.3,8.9,5.2,10.9,14.4)
# Define pdf, cdf , inverse cdf
dgumbel <- function(x,a,b) 1/b*exp((a-x)/b)*exp(-exp((a-x)/b))
pgumbel <- function(q,a,b) exp(-exp((a-q)/b))
qgumbel <- function(p,a,b) a-b*log(-log(p))
# Fit with MLE
f1c <- fitdist(x1,"gumbel",start=list(a=10,b=5))
# Goodness of Fit
gofstat(f1c, fitnames = 'Gumbel MLE')
参考： https://www.rdocumentation.org/packages/fitdistrplus/versions/0.2-1/topics/fitdist
如果提供的示例不再有效，这是另一种解决方案：
library(evd)
fitgev <- fgev(month)
# e.g. extract log-likelihood
logLik(fitgev)[[1]]