r - R 的 powerlaw 包中连续与离散对数正态分布的似然函数差异

标签 r distribution data-fitting power-law

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

当我拟合 xmin(低于该阈值的计数值将被忽略)、记录均值和记录 sd 参数并引导结果以获得 p 值时,我收到向量内存耗尽错误。我发现当包内部函数 sample_p_helper 尝试从拟合分布生成随机数时会发生这种情况。拟合的对数均值和对数标准差参数非常低,以至于拒绝采样算法试图生成数十亿个数字以获得高于 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?)

问题就变成了为什么首先要拟合如此低且拟合不良的对数均值和对数标准差参数值。

我注意到,如果我拟合对数正态分布的连续版本,则不会发生错误,并且参数值看起来更合理(事实上,p值表明数据与对数正态分布):

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

bs <- bootstrap_p(dist_cont)

bs

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

The continuous version看起来如我所料:

########################################################
#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

}

但是,在 the discrete version ,联合概率的计算方式不同:

########################################################
#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 个值,我知道,这对于最大似然来说太少了,不可靠。

感谢您耐心看完这篇冗长的文章。我知道这既是一个统计问题,也是一个编码问题。非常感谢任何朝着正确方向有帮助的插入!干杯。

最佳答案

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

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

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

  4. 是的,您的数据就是问题所在。但这也是模型不起作用时的问题;)

关于r - R 的 powerlaw 包中连续与离散对数正态分布的似然函数差异,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/57721406/

相关文章:

r - 列表中列的数据类型 - R

r - 将facet_wrap 与for 循环一起使用

r - “finalfit”在处理连续因变量时丢失标签

r - 由于地理单位不同,我在向 ggmap 添加 shapefile 时遇到问题

r - 声明一个用户定义的分布

c++ - 在 C/C++ 中按照正态分布生成随机数

r - 如何从 glm.nb() 中获取 "prob"参数?

gnuplot - 我想拟合一个带有整数参数的函数

python - 使用 Symfit : typestructure of dataset 进行全局拟合

java - 分发我的 JavaGameEngine(使用 Jython)