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

Question

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

data <- rplcon(1000,10,2)

现在，我想知道哪些已知分布最适合数据。对数范数？经验值？ Gamma ？幂律？具有指数截断的幂律？

所以我使用函数 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，幂律分布和指数截断幂律不是基本概率函数, 所以我根据?fitdist 的例子4写了幂律的d,p,q函数

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 统计量，参见 this 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)