1

我有一个带有观察到的 bin 计数的直方图。我想根据观察到的计数运行模拟，以了解相同数量的观察结果可能如何不同地发生。我将直方图转换为向量，并将观察到的计数作为向量的一个元素。rbinom(n, size, prob)我使用基于 bin 频率的概率从二项式分布（来自 ）生成的随机数来模拟每个 bin 。

我的问题是模拟观察计数为零的箱。当 bin 计数为零时，prob=0，因此该 bin 的模拟计数始终为零。这是非物质的，不是我想要的。目前，我通过用 1 的 bin 计数覆盖零 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)
}
4

2 回答 2

0

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

与您正在做的类似的一种简单方法称为拉普拉斯的“继承规则”：只需在每个 bin 中添加一个（不仅是零）。更一般地，您可以在附加平滑中添加不同的数字（通常小于 1）。贝叶斯，这对应于在二项式概率上使用狄利克雷先验。

更复杂的方法是Good-Turing 平滑。

如果您的数据有不同的概率模型，在您的设置中更有意义，您也可以使用它（如 Eric 建议的那样）。

于 2012-04-24T18:05:45.620 回答
0

您尝试做的事情听起来很像自举，如果您从一个包含 n 个值的数组开始，您会从该数组中随机选择 n 个值并进行替换。正如您所注意到的，引导不会给您一个您还没有的价值。

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

另一种选择是使用参数分布来拟合您的 bin。您所在领域的人是否有他们使用的典型分布，或者数据是否暗示了分布？如果您仍然需要 bin，则可以对拟合的参数分布进行 bin。祝你好运。

这是生成垃圾箱计数的更快方法：

# 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))
于 2012-04-24T17:41:04.417 回答