2

我正在尝试使用 Colin GillespiepoweRlaw在 R 中的包将对数正态分布拟合到一些计数数据。我知道对数正态分布是连续的并且计数数据是离散的，但是，该包包含连续和离散版本的类和方法的对数正态分布。

当我拟合 xmin（忽略计数值的阈值以下）、记录均值和记录 sd 参数并引导结果以获得p值时，我得到一个向量内存耗尽错误。我发现当包内部函数sample_p_helper试图从拟合分布中生成随机数时会发生这种情况。拟合的 log mean 和 log sd 参数非常低，以至于拒绝采样算法试图生成数十亿个数字来获得高于 xmin 的任何值，因此存在内存问题。

输入：

library(poweRlaw)

counts <- c(54, 64, 126, 161, 162, 278, 281, 293, 296, 302, 322, 348, 418, 511, 696, 793, 1894)

dist <- dislnorm$new(counts)      # Create discrete lnorm object
dist$setXmin(estimate_xmin(dist)) # Get xmin and parameters

bs <- bootstrap_p(dist) # Run bootstrapping

错误信息：

Expected total run time for 100 sims, using 1 threads is 24.6 seconds.
Error in checkForRemoteErrors(val) : 
  one node produced an error: vector memory exhausted (limit reached?)

那么问题就变成了为什么首先要拟合如此低且拟合不佳的 log mean 和 log sd 参数值。

我注意到，如果我拟合对数正态分布的连续版本，则不会出现错误，并且参数值似乎更合理（实际上，p值​​表明数据与对数正态分布兼容）：

dist_cont <- conlnorm$new(counts)
dist_cont$setXmin(estimate_xmin(dist_cont))

bs <- bootstrap_p(dist_cont)

bs

查看包的源代码，我注意到离散和连续对数正态分布的似然函数是不同的。具体来说就是计算联合概率的部分。

连续版本看起来像我期望的那样：

########################################################
#Log-likelihood 
########################################################
conlnorm_tail_ll = function(x, pars, xmin) {
  if(is.vector(pars)) pars = t(as.matrix(pars))
  n = length(x)
  joint_prob = colSums(apply(pars, 1, 
                             function(i) dlnorm(x, i[1], i[2], log=TRUE)))

  prob_over = apply(pars, 1, function(i) 
    plnorm(xmin, i[1], i[2], log.p=TRUE, lower.tail=FALSE))
  joint_prob - n*prob_over

}

然而，在离散版本中，联合概率的计算方式不同：

########################################################
#Log-likelihood 
########################################################
dis_lnorm_tail_ll = function(xv, xf, pars, xmin) {
  if(is.vector(pars)) pars = t(as.matrix(pars))
  n = sum(xf)
  p = function(par) {
    m_log = par[1]; sd_log = par[2]
    plnorm(xv-0.5, m_log, sd_log, lower.tail=FALSE) - 
      plnorm(xv+0.5, m_log, sd_log, lower.tail=FALSE)
  }
  if(length(xv) == 1L) {
    joint_prob = sum(xf * log(apply(pars, 1, p)))
  } else {
    joint_prob = colSums(xf * log(apply(pars, 1, p)))
  }
  prob_over = apply(pars, 1, function(i) 
    plnorm(xmin-0.5, i[1], i[2], 
           lower.tail = FALSE, log.p = TRUE))

  return(joint_prob - n*prob_over)
}

指数分布的离散和连续实现之间存在类似的差异，但离散和连续幂律分布没有。在连续版本中，joint_prob通过相对简单的调用来计算dlnorm，但离散版本调用plnorm。此外，他们调用plnorm了两次，首先是观察到的数据值 -0.5，然后是观察到的值 +0.5，然后从后者中减去前者。

所以，最后，我的问题：

  • 为什么powerlaw在对数正态分布的离散实现中要这样计算联合概率呢？我敢肯定它以这种方式编写是有原因的，这只是我的数学无知，但我并不真正理解它。

  • 即使我的数据是离散的，使用 powerlaw 的连续对数正态分布是否安全，因为无论如何它似乎工作得很好？

  • 在尝试拟合离散对数正态分布时，关于我的数据可能出现什么问题的任何其他线索？我在想某处可能存在扩展问题，但很难解决这个问题。

  • 我的可笑的小数据集会起作用吗？我正在尝试将分布拟合到仅高于 xmin 的 8 个值，这对于最大似然性来说太少了，我知道。

感谢您通过这篇冗长的帖子忍受我。我知道这既是一个统计问题，也是一个编码问题。非常感谢任何朝着正确方向有帮助的推动！干杯。

4

1 回答 1

2

  1. dlnorm()给出概率密度值. 请记住，密度积分为一，但不要总和为一。因此，为了计算离散分布，我们取整数两侧的值。它们也将是一个归一化常数。对于 CTN 情况，对数似然只是 的乘积dlnorm()，这样更容易更快。

  2. “安全”是一个很难定义的词。对于该数据，CTN 和离散在视觉上给出了相同的拟合。但两者都不合适。

  3. 离散分布的估计参数值在非常极端的尾部给出截断的对数正态分布。模拟该地区的数据具有挑战性

  4. 是的，您的数据就是问题所在。但这也是模型不起作用时的问题；）

于 2019-08-30T18:34:15.207 回答