由 `fitdistrplus` 包中的 `fitdistr()` 函数拟合的幂律

Question

我使用rplcon()包中的函数生成一些随机变量poweRlaw

data <- rplcon(1000,10,2)

现在，我想知道哪些已知分布最适合数据。对数范数？经验？伽玛？幂律？指数截止的幂律？

所以我fitdist()在包中使用函数fitdistrplus：

fit.lnormdl <- fitdist(data,"lnorm")
fit.gammadl <- fitdist(data, "gamma", lower = c(0, 0))
fit.expdl <- fitdist(data,"exp")

由于幂律分布和具有指数截止的幂律不是根据CRAN Task View: Probability Distributions的基本概率函数，所以我根据示例 4 编写了幂律的 d,p,q 函数?fitdist

dplcon <- function (x, xmin, alpha, log = FALSE) 
{
    if (log) {
        pdf = log(alpha - 1) - log(xmin) - alpha * (log(x/xmin))
        pdf[x < xmin] = -Inf
    }
    else {
        pdf = (alpha - 1)/xmin * (x/xmin)^(-alpha)
        pdf[x < xmin] = 0
    }
    pdf
}
pplcon <- function (q, xmin, alpha, lower.tail = TRUE) 
{
    cdf = 1 - (q/xmin)^(-alpha + 1)
    if (!lower.tail) 
        cdf = 1 - cdf
    cdf[q < round(xmin)] = 0
    cdf
}
qplcon <- function(p,xmin,alpha) alpha*p^(1/(1-xmin))

最后，我使用下面的代码来获取参数xmin和alpha幂律：

fitpl <- fitdist(data,"plcon",start = list(xmin=1,alpha=1))

但它抛出一个错误：

<simpleError in optim(par = vstart, fn = fnobj, fix.arg = fix.arg, obs = data,     ddistnam = ddistname, hessian = TRUE, method = meth, lower = lower,     upper = upper, ...): function cannot be evaluated at initial parameters>
Error in fitdist(data, "plcon", start = list(xmin = 1, alpha = 1)) : 
  the function mle failed to estimate the parameters, 
                with the error code 100

我尝试在google和stackoverflow中搜索，出现了这么多类似的错误问题，但是经过阅读和尝试，我的问题没有解决方案，我该怎么做才能正确完成以获取参数？感谢所有帮过我的人！

Answer 1

这是一个有趣的发现，我对这个发现并不完全满意，但我会告诉你我的发现，看看它是否有帮助。

在调用该fitdist函数时，默认情况下它想mledist从同一个包中使用。这本身会导致调用stats::optim通用优化函数。在它的返回值中，它给出了一个收敛错误代码，?optim有关详细信息，请参阅。在100你看到的是不返回的人之一optim。因此，我将代码拆开mledist并fitdist找出该错误代码的来源。不幸的是，它在不止一种情况下被定义并且是一个通用的陷阱错误代码。如果您分解所有代码，fitdist则此处尝试执行以下操作，并事先进行各种检查等。

fnobj <- function(par, fix.arg, obs, ddistnam) {
  -sum(do.call(ddistnam, c(list(obs), as.list(par), 
                           as.list(fix.arg), log = TRUE)))
}

vstart = list(xmin=5,alpha=5)
fnobj <- function(par, fix.arg obs, ddistnam) {
  -sum(do.call(ddistnam, c(list(obs), as.list(par), 
                           as.list(fix.arg), log = TRUE)))
}
ddistname=dplcon
fix.arg = NULL
meth = "Nelder-Mead"
lower = -Inf
upper = Inf
optim(par = vstart, fn = fnobj, 
      fix.arg = fix.arg, obs = data, ddistnam = ddistname, 
      hessian = TRUE, method = meth, lower = lower, 
      upper = upper)

如果我们运行这段代码，我们会发现一个更有用的错误“无法在初始参数处评估函数”。如果我们查看函数定义，这是有道理的。具有xmin=0或alpha=1将产生 的对数似然-Inf。好的，所以想尝试不同的初始值，我尝试了一些随机选择，但都返回了一个新错误，“非有限有限差分值1 ”。

optim进一步搜索这两个错误的来源，它们不是 R 源本身的一部分，但是有一个.External2调用，所以我只能假设错误来自那里。非有限错误意味着某处的函数评估之一给出了非数字结果。该函数dplcon将在alpha <= 1或时这样做xmin <= 0。fitdist允许您指定传递给的mledist其他参数或其他参数（取决于您选择的方法，默认为 mle），其中lower之一用于控制要优化的参数的下限。所以我尝试施加这些限制并再次尝试：

fitpl <- fitdist(data,"plcon",start = list(xmin=1,alpha=2), lower = c(xmin = 0, alpha = 1))

令人讨厌的是，这仍然给出了错误代码 100。跟踪它会产生错误“L-BFGS-B 需要 'fn' 的有限值”。当您指定边界时，优化方法已从默认的 Nelder-Mead 更改，并且在外部 C 代码调用的某处出现此错误，可能接近任一xmin或alpha我们接近无穷大时数值计算的稳定性很重要的限制。

我决定做分位数匹配而不是最大似然来尝试找出更多

fitpl <- fitdist(data,"plcon",start = list(xmin=1,alpha=2),
    method= "qme",probs = c(1/3,2/3))
fitpl
## Fitting of the distribution ' plcon ' by matching quantiles 
## Parameters:
##          estimate
## xmin   0.02135157
## alpha 46.65914353

这表明 的最佳值xmin接近 0，这是极限。我不满意的原因是我无法得到分布的最大似然拟合，fitdist但希望这个解释有帮助，分位数匹配提供了一个替代方案。

编辑：

在学习了更多关于幂律分布的一般知识后，这并不像您期望的那样工作是有道理的。参数 power 参数有一个似然函数，它可以在给定的 xmin 条件下最大化。然而，xmin 不存在这样的表达式，因为似然函数在 xmin 中增加。通常，xmin 的估计来自 Kolmogorov--Smirnov 统计，请参阅此mathoverflow 问题和 poweRlaw 包的 d_jss_paper 小插图以获取更多信息和相关参考资料。

poweRlaw包本身具有估计幂律分布参数的功能。

m = conpl$new(data)
xminhat = estimate_xmin(m)$xmin
m$setXmin(xminhat)
alphahat = estimate_pars(m)$pars
c(xmin = xminhat, alpha = alphahat)