r - 通过扰乱现有概率分布来生成离散随机概率分布

Question

如果我想有效地生成一个总和为 1 的 N 个概率的随机离散概率分布，我可以接受 Hadley 的评论 here :

prop.table(runif(N))

如果我重复多次，N 个元素中每个元素的平均概率应该是 ~1/N。

如果我希望 N 个元素中的每一个的平均概率不是 1/N 而是一个指定的数字a priori怎么办？

例如N = 4 元素，我有 apriori 分布:

apriori <- c(0.2, 0.3, 0.1, 0.4)

我想要随机分布基于这个先验，例如:

c(0.21, 0.29, 0.12, 0.38)
c(0.19, 0.29, 0.08, 0.44)
c(0.19, 0.33, 0.1, 0.38)

等等

我们遵循以下任一规则的地方:

1) 平均每个元素的概率将是(近似)它在先验分布中的概率

2) 有一个“扰动”参数，比如说 perturbation = 0.05 这意味着:(a) 我们让每个概率 i 都在apriori[i] +- perturbation 范围或 (b) 我们让每个概率 i 都在 apriori[i] +- perturbation *先验[i] 范围(即先验概率的正负 5%，不是绝对的 5%)

我不知道如何在遵守规则 1 的同时做到这一点。

关于规则 2，我最初的低效想法是以随机允许的数量扰动前 N - 1 个元素中的每一个，将最后一个元素设置为 1 - sum(N-1_probs) 并且用 while 循环包装它，直到最后一个元素也是合法的。

我什至还没有实现它，因为那是非常低效的(比如我想要 100K 的这样的分布......)。想法？

Answer 1

根据 prof.Bolker 的建议，您应该查看 Dirichlet distribution .让我们用大写字母 C_i 表示平均先验值，用小写字母 c_i 表示采样值。它将自动从分发属性中为您提供两个功能:

1. 求和 _i c_i = 1

2. 每个c_i都在[0...1]范围内

所以您可以立即将它们用作概率。

给定 C_i，并查看分布定义(检查链接)，剩下的唯一自由参数是

a₀ = _i a_i

之和

并且每个 a_i = C_i * a₀

a_i 的这种选择将(再次自动地)提供适当的平均值 E[c_i] = C_i。

越大的 a₀ - c_i 在 C_i 周围会更窄。方差大致是 Var[c_i] ~ C_i/a₀，所以对于 5% 你可以尝试使用 a ₀，共 50 个。

一些R代码

library(MCMCpack)

apriori <- c(0.2, 0.3, 0.1, 0.4) # your C_i
a0 <- 50
a <- a0*apriori

set.seed(12345)
# sample your c_i and use it, for example, to throw uneven dice
ci <- rdirichlet(1, a)
dice <- rmultinom(1, 1, ci)

# another dice throw
ci <- rdirichlet(1, a)
dice <- rmultinom(1, 1, ci)

...