3

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 任务视图：概率分布的基本概率函数，所以我根据示例 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中搜索，出现了很多类似的错误问题，但是在阅读和尝试之后，我的问题没有解决方案，我应该怎么做才能正确完成以获得参数？感谢所有帮助我的人！

4

1 回答 1

5

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

在调用该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 问题和 powRlaw 包的 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)
于 2016-05-11T08:40:02.980 回答