如何从非相同分布的伯努利随机变量中矢量化采样？答案

【问题标题】：How to vectorise sampling from non-identically distributed Bernoulli random variables?如何从非相同分布的伯努利随机变量中矢量化采样？
【发布时间】：2017-06-10 03:12:20
【问题描述】：

给定一系列独立但非同分布的伯努利试验，成功概率由向量给出，例如：

x <- seq(0, 50, 0.1)
prob <- - x*(x - 50)/1000 # trial probabilities for trials 1 to 501

从每次试验中获取随机变量的最有效方法是什么？我假设矢量化是要走的路。

我知道两个给伯努利随机变量的函数：

rbernoulli 来自包 purr，它不接受成功概率向量作为输入。在这种情况下，可以将函数包装在 apply 类型的操作中。
rbinom 带参数size = 1 给出伯努利随机变量。它还接受一个概率向量，因此：

rbinom(n = length(prob), size = 1, prob = prob)

给出正确长度的输出。但是，我不完全确定这实际上是我想要的。帮助文件?rbinom 中似乎相关的位是：

结果的长度由 rbinom 的 n 决定，并且是另一个数值参数的最大长度功能。

n 以外的数值参数被回收到结果。仅使用逻辑参数的第一个元素。

但是，n 是一个没有默认值的参数，所以我不确定第一句是什么意思。我认为第二句话意味着我得到了我想要的，因为只有 size = 1 应该被回收。但是this thread 似乎暗示这种方法不起作用。

这个blog post 还提供了一些其他方法。一位评论员使用rbinom 提到了我的建议。

【问题讨论】：

你可以在 purrr 中迭代和简化：library(purrr); map_lgl(prob, ~rbernoulli(1, p = .x))
不过，这应该等同于（并且比）as.logical(rbinom(length(prob), 1, prob))。你可以看到它适用于rbinom(500, 1, seq(0,1, len=500))。
我可以看到 rbinom 生成的东西看起来像我想要的，但这让我担心互联网上的其他地方似乎有很多其他方法的参考，我恐怕rbinom 的文档对我来说不是很清楚。
n 是您返回的数字数量。 size 是试验次数，其中成功意味着将 1 添加到结果中，因此如果 size = 1，您将针对返回的每个数字运行一次试验。 prob 是每个单独试验的成功机会，并且是矢量化的（如果相关，可以回收）。你已经得到了答案；你只需要相信它。或者测试一下，如果有帮助的话。

标签： r random probability

【解决方案1】：

如果您不信任互联网上的随机陌生人并且不了解文档，也许您可以通过测试来说服自己。只需设置随机种子即可获得可重现的结果：

x <- seq(0, 50, 0.1)
prob <- - x*(x - 50)/1000

#501 seperate draws of 1 random number
set.seed(42)
res1 <- sapply(prob, rbinom, n = 1, size = 1)

#501 "simultaneous" (vectorized) draws
set.seed(42)
res2 <- rbinom(501, 1, prob)

identical(res1, res2)
#[1] TRUE

【讨论】：

【解决方案2】：

另一种测试rbinom 是否为prob 向量化的方法，利用N 个伯努利随机变量之和是分母为N 的二项式随机变量这一事实：

x <- seq(0, 50, 0.1)
prob <- -x*(x - 50)/1000

n <- rbinom(prob, size=1000, prob)

par(mfrow=c(1, 2))
plot(prob ~ x)
plot(n ~ x)

【讨论】：