基于修正的伯努利分布生成整数序列

如何解决基于修正的伯努利分布生成整数序列

我想用 R 随机生成一个整数序列，每个整数都是从整数池 (0,1,2,3....,k) 中挑选出来的，并带有替换。 k 是预先确定的。 (0,k) 中每个整数 k 的选择概率是 p^k(1-p) 其中 p 是预先确定的。也就是说，与 k 相比，1 被选中的概率要高得多，而且我的最终整数序列可能比 k 有更多的 1。我不知道如何在 R 中实现这个数字选择过程。

解决方法

解决此类问题的通用方法是：

计算每个整数的 p^k * (1-p)
在表 t 中创建这些的累积总和。
用 range(t) 绘制一个均匀分布的数字
测量该数字落入 t 的距离并检查对应于哪个整数。
一个整数的概率越大，它覆盖的范围就越大。

这是快速而肮脏的示例代码：

draw <- function(n=1,k,p) {
    v <- seq( 0,k )
    pr <- (p ** v) * (1-p)
    t <- cumsum(pr)
    r <- range(t)
    x <- runif( n,min=min(r),max=max(r) )
    f <- findInterval( x,vec=t )
    v[ f+1 ] ## first interval is 0,and it will likely never pass highest interval
}

请注意，建议的解决方案并不关心您的密度函数加起来是否为 1。根据您的描述，在现实生活中很可能会如此。但这对于解决方案来说并不重要。

Sirius 的回答很好。但据我所知，您所描述的内容类似于截断的 geometric distribution。

我应该注意到几何分布在不同的作品中定义不同（例如见MathWorld），所以我们使用如下定义的分布：

P(X = x) ~ p^x * (1 - p)，其中 x 是 [0,k] 中的整数。

我对 R 不是很熟悉，但解决方案是调用 rgeom(1,1 - p) 直到结果为 k 或更少。

或者，您可以使用通用拒绝采样器，因为概率是已知的（这里最好称为权重，因为它们不需要总和为 1）。拒绝抽样描述如下：

将权重存储在列表中。计算最高权重，称其为max。然后，使用拒绝采样在区间 [0,k] 中选择一个整数：

在区间 [0,i] 中选择一个统一的随机整数 k。
以概率 weights[i]/max（在您的情况下为 weights[i] = p^i * (1-p)），返回 i。否则，请转到第 1 步。

给定每个项目的权重，除了拒绝抽样或Sirius答案中的解决方案之外，还有许多其他方法可以做出加权选择；看我的note on weighted choice algorithms。