模拟R中直方图二项式起伏的零计数箱

Question

我有一个直方图，上面有观察到的垃圾桶数量。我想根据观察到的计数运行模拟，看看相同数量的观察如何以不同的方式发生。我将直方图转换为向量，并将观察到的计数作为向量的一个元素。我使用从二项式分布(来自rbinom(n, size, prob))生成的随机数来模拟每个bin，这些分布具有基于bin频率的概率。

我的问题是模拟观察计数为零的垃圾箱。当仓位计数为零时，prob=0，因此该仓位的模拟计数始终为零。这是非物质的，不是我想要的。目前，我通过用bin计数1覆盖零bin计数来处理这个问题。我不确定这样做的效果，所以我不知道我的模拟是否超出了我的容忍度。我正在寻找一个比我的特别方法更好或更优雅的解决方案。

有什么想法吗？谢谢。

下面是我的相关代码：

sim.vector <- function(x, n = length(x)) {
  sum.vector <- round(sum(x), 0)  # the number of observations
  x.dummy <- x
  x.dummy[x.dummy == 0] <- 1  # override bins with zero counts
  f <- x.dummy / sum(x) # the frequency of each bin
  x.sim <- rep(0, n)
  while(sum.vector != sum(x.sim)) {  # make sure the simulation has the same
                                     # number of total counts as the observation
    for (i in 1:n) {
      p.target <- f[i]  # set the probability of each bin to the frequency
      x.sim[i] <- rbinom(1, sum.vector, p.target)  # create a random binomial
    }
  }
  return(x.sim)
}

Answer 1

您尝试做的事情听起来很像bootstrapping，如果您从一个包含n个值的数组开始，您将从该数组中随机选择n个值并进行替换。正如您已经注意到的，bootstrapping不会为您提供您尚未拥有的值。

您设置零箱的方法是可以的。生态学家使用的一种技术是将零值设置为他们可能产生的最小测量误差。例如，如果计数树，最小错误计数将是1。如果保持分布的平均值对您很重要，请确保将0更改为1不会使平均值增加太多。

另一种选择是使用参数分布来拟合回收箱。你所在领域的人是否有他们使用的典型分布，或者数据是否建议分布？如果仍然需要存储箱，可以将拟合的参数分布存储在存储箱中。祝好运。

这里有一个更快的方法来生成你的垃圾桶计数：

# you have three bins which we label 1,2,3
# change to 1:n where n = number of your bins
binIndex = 1:3
# you have 2 values in bin 1, 3 values in bin 2, 1 value in bin 3
count1 = c(2,3,1)

# create a vector with all bin labels
unbin = mapply(function(bin1,ct1){
    rep(bin1,times=ct1)
},binIndex,count1)

unbin = unlist(unbin)
print(unbin)

# generate "bootstrapBinCount" of bootstrapped bins
bootstrapBinCount = 10

# in general, use lapply instead of while/for loops in R - it's orders of magnitude faster
# newBins is a list of binCounts
# to access the first bin count, try newBins[[1]]
newBins = lapply(1:bootstrapBinCount,function(x){
    # generate a bootstrap from the list of bin labels by sampling with replacement
    bootstrap1 = sample(unbin, size=length(unbin), replace = TRUE)
    # count the number of times each bin label shows up
    rebin = table(bootstrap1)
    # get the names of the labels that showed up
    names1 = as.integer(names(rebin))
    # fill in the counts from rebin, but also include cases where there are no values     for a given bin label
    rebinPlusZeroes = rep(0,times=length(binIndex))
    rebinPlusZeroes[names1] = as.integer(rebin)

    return(rebinPlusZeroes)
})

print(str(newBins))

Answer 2

在频率估计中解决零计数的问题有时被称为平滑(特别是在自然语言处理社区中，他们处理数千个“bin”，因此零计数很常见)；您可以使用pseudocounts而不是使用计数。

与您所做的类似的一种简单方法称为Laplace's "rule of succession"：只需在每个bin中添加1(而不仅仅是0)。更一般地，您可以在additive smoothing中添加一个不同的数字(通常小于1)。就贝叶斯而言，这相当于在二项式概率上使用Dirichlet先验。

一种更复杂的方法是Good-Turing smoothing。

如果您的数据具有不同的概率模型，并且在您的设置中更有意义，您也可以使用该模型(如Eric所建议的)。