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

Question

如果我想有效地生成 N 个概率总和为 1 的随机离散概率分布，我可以接受 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" 参数，比如说 perturbation = 0.05 这意味着：(a) 我们让每个概率 i 在 apriori[i] +- perturbation 范围或 (b) 我们让每个概率 i 都在 apriori[i] +- perturbation * apriori[i] 范围内（即 plus/minus 该先验概率的 5%，而不是绝对的 5%）

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

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

我什至还没有实现它，因为那是非常低效的（比如我想要 100K 个这样的发行版......）。想法？

Answer 1

并为每个概率使用正态分布？

perturbation <- 0.05
plouf <- sapply(apriori,function(x){max(rnorm(1,mean = x, sd = perturbation*x),0)})
plouf <- plouf/sum(plouf)
> plouf
[1] 0.2020629 0.3057111 0.0994482 0.3927778

Answer 2

我有一个解决办法，但最终会出现抽签正常的情况。我认为您可以做类似的事情来绘制均匀分布。在这方面没有太多经验，但我会倾向于拒绝类型的政策，您可以快速绘制很多东西，然后拒绝那些不符合您标准的东西

rm(list = ls())

library(parallel)
library(data.table)
library(tictoc)

# set up the distribution informatoin
P <- 4
values <- 1:P
dist_scores <- data.table(param = values,
                          prob = c(0.2, 0.3, 0.1, 0.4), key = "param")
perturbation <- 0.05
method = "a"

switch (method,
  "a" = {dist_scores[, min := prob - perturbation]
    dist_scores[, max := prob + perturbation]},
  "b" = {dist_scores[, min := prob * (1-perturbation)]
    dist_scores[, max := prob * (1+perturbation)]}
)

# turn this in to a set of data that can be sampled
N <- 10000
v <- unlist(sapply(values, FUN = function(x){
  rep(x, round(dist_scores$prob[x]*N, 0))
}))
table(v)/N

# set number of samples, and number of draws for each iteration
sams <- 10000
reps <- 200

tic()
# loop through and draw reps from the sample. Rejection policy will remove
# ones that dont meet the conditions
new_iters <- mclapply(1:sams, FUN = function(x){
  y <- data.table(param = sample(v, reps, replace = TRUE))
  out <- y[, .(val = .N/reps), keyby = param]
  out <- dist_scores[out,]
  if(out[,all(val >= min & val <= max)]){
    return(out[, c("param", "val"), with = FALSE])
  }else{
    return(NULL)
  }
})
reject_rate <- sum(sapply(new_iters, is.null))/sams
# number of samples
sams - reject_rate*sams
toc()

out <- rbindlist(new_iters)

par(mfrow = c(2,2))
for(i in values){
  hist(out[param == i, val])
}enter code here

Answer 3

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

总和 _i c_i = 1
每个c_i在[0...1]范围内

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

给定 C_i，并查看分布定义（检查 link），剩下的唯一可用参数是

a₀ = 总和 _i a_i

和每个a_i = C_i * a₀

这样选择 a_i 将（再次自动）提供适当的平均值 E[c_i] = C_我.

更大的 a₀ - c_i 在 C_i 周围会更窄。方差大致来说就是Var[c_i] ~ C_i/a₀,所以为5%你可以尝试使用 a₀ of 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)

...

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

generating a discrete random probability distribution, by perturbing an existing one

random

algorithm

r

probability