logo标签列表

在R中生成三个随机数，使它们相加等于1

rrandom
13

我希望创建三个(非负)准随机数，它们相加为一，并且可以一遍又一遍地重复。

基本上，我正在尝试在多次试验中将某物分成三个随机部分。

虽然我知道 

a = runif(3,0,1)

我认为在下一个 runif 中可以使用 1-a 作为最大值，但这看起来有点混乱。

当然，这些数的总和并不为一。有没有什么想法呢，聪明的stackoverflow用户们？

2
生成随机数后重新归一化是否是一个选项？ - Anders Gustafsson
生成两个随机数a和b怎么样？然后a+b+c = 1 => c = 1 - (a+b) - Frank Schmitt
如果a和b的总和大于1呢？ - mmann1123
你想在间隔上使用什么分布？例如，第一个间隔的长度是否应该在（0,1）上均匀选择？ - ALiX
你所遇到的问题似乎被称为“区间的随机分割”。这是统计学中的一个经典问题，但在谷歌上搜索它时，出现的漂亮图片却很少... - Christian
7个回答
16

这个问题涉及到比一开始看起来更微妙的问题。在查看以下内容后，您可能需要仔细考虑您正在使用这些数字所代表的过程：

## My initial idea (and commenter Anders Gustafsson's):
## Sample 3 random numbers from [0,1], sum them, and normalize
jobFun <- function(n) {
    m <- matrix(runif(3*n,0,1), ncol=3)
    m<- sweep(m, 1, rowSums(m), FUN="/")
    m
}

## Andrie's solution. Sample 1 number from [0,1], then break upper 
## interval in two. (aka "Broken stick" distribution).
andFun <- function(n){
  x1 <- runif(n)
  x2 <- runif(n)*(1-x1)
  matrix(c(x1, x2, 1-(x1+x2)), ncol=3)
}

## ddzialak's solution (vectorized by me)
ddzFun <- function(n) {
    a <- runif(n, 0, 1)
    b <- runif(n, 0, 1)
    rand1 = pmin(a, b)
    rand2 = abs(a - b)
    rand3 = 1 - pmax(a, b)
    cbind(rand1, rand2, rand3)
}
    
## Simulate 10k triplets using each of the functions above
JOB <- jobFun(10000)
AND <- andFun(10000)
DDZ <- ddzFun(10000)

## Plot the distributions of values
par(mfcol=c(2,2))
hist(JOB, main="JOB")
hist(AND, main="AND")
hist(DDZ, main="DDZ")

在此输入图像描述

不错，我正在考虑绘制结果，但你已经做了这个。有趣的是，显然没有一个解决方案真正做到了直觉上想要的事情。有趣的是，在这些图中，你实际上看不出DDZ根据平均值做了正确的事情，而AND甚至没有做到这一点。 - Christian
Stefan Jelkovich -- 你能否修改你建议的编辑中的图表，使它们的x轴从0到1，并在帖子正文中指出哪部分是你添加的内容？然后重新提交，我会接受这些编辑。 （你可以通过点击“编辑”按钮并查看它显示的更改历史记录来恢复你已经进行的编辑。）谢谢。 - Josh O'Brien
经过仔细查看，我撤销了您的编辑，因为您提供的结果并没有回答这篇文章所要解决的问题。（请查看您图表的x轴上的概率，以了解您的解决方案存在何种问题。） - Josh O'Brien
11

从（0，1）中随机选择两个数字，如果我们假设它们为a和b，则得到：

rand1 = min(a, b)
rand2 = abs(a - b)
rand3 = 1 - max(a, b)
此外，如果 a == b（这应该是非常罕见的情况），您必须重新生成第二个数字。 - ddzialak
@user 如果a=0.85，b=0.99，那么你得到的数字是：0.85、0.14、0.01（对我来说，这是从0到1中非常好的3个随机数）。 - ddzialak
3
产生的分布似乎并不是很平凡：http://www.jstor.org/discover/10.2307/2983572?uid=2129&uid=2&uid=70&uid=4&sid=21100849643501，以及一篇后来可以免费访问的论文http://doc.utwente.nl/70657/1/Sleutel67random.pdf - Christian
9

如果您想随机生成总和为1（或其他值）的数字，那么您应该查看狄利克雷分布。

gtools包中有一个rdirichlet函数，运行RSiteSearch('Dirichlet')会出现很多结果，这些工具可以帮助您进行操作（对于简单的狄利克雷分布，手动编写代码也不难）。

6
我想这取决于您想要的数字分配方式，但这里有一种方法：
diff(c(0, sort(runif(2)), 1))

使用replicate可以获取所需数量的集合：

> x <- replicate(5, diff(c(0, sort(runif(2)), 1)))
> x
           [,1]       [,2]      [,3]      [,4]       [,5]
[1,] 0.66855903 0.01338052 0.3722026 0.4299087 0.67537181
[2,] 0.32130979 0.69666871 0.2670380 0.3359640 0.25860581
[3,] 0.01013117 0.28995078 0.3607594 0.2341273 0.06602238
> colSums(x)
[1] 1 1 1 1 1
5

我会从均匀分布中随机选择3个数字，然后将它们相加并求商：

n <- 3
x <- runif(n, 0, 1)
y <- x / sum(x)
sum(y) == 1

n可以是任何你喜欢的数字。

2

这个问题以及提出的不同解决方案都让我感到困惑。我对建议的三种基本算法进行了一些测试，并计算出它们对生成的数字的平均值。

choose_one_and_divide_rest
means:                [ 0.49999212  0.24982403  0.25018384]
standard deviations:  [ 0.28849948  0.22032758  0.22049302]
time needed to fill array of size 1000000 was 26.874945879 seconds

choose_two_points_and_use_intervals
means:                [ 0.33301421  0.33392816  0.33305763]
standard deviations:  [ 0.23565652  0.23579615  0.23554689]
time needed to fill array of size 1000000 was 28.8600130081 seconds

choose_three_and_normalize
means:                [ 0.33334531  0.33336692  0.33328777]
standard deviations:  [ 0.17964206  0.17974085  0.17968462]
time needed to fill array of size 1000000 was 27.4301018715 seconds

时间测量需要谨慎对待，因为它们可能更受Python内存管理的影响，而不是算法本身。我太懒了，不想用 timeit 做正确的事情。我在1GHz Atom上做的，这就解释了为什么花了那么长时间。
无论如何，choose_one_and_divide_rest 是Andrie和提问者自己（AND）建议的算法：你选择[0,1]中的一个值a，然后在[a,1]中选择一个值，再看看你剩下了什么。它加起来是一的，但仅此而已，第一次划分比其他两次大两倍。人们可能已经猜到了...
choose_two_points_and_use_intervals是ddzialak（DDZ）的答案。它选取区间[0,1]中的两个点，并使用由这些点创建的三个子区间的大小作为三个数字。像魔法一样工作，平均值都是1/3。
choose_three_and_normalize是Anders Gustafsson和Josh O'Brien（JOB）的解决方案。它只生成[0,1]中的三个数字，并将它们归一化为总和为1。它同样有效，而且在我的Python实现中令人惊讶地更快一些。方差比第二个解决方案略低。
就是这样。我不知道这些解决方案对应哪个beta分布或对应论文中的哪组参数，但也许其他人可以弄清楚。
1
最简单的解决方案是使用 Wakefield 包中的 probs() 函数。
调用 probs(3) 将返回一个包含三个值且总和为1的向量。
可以使用 rep(probs(3),x) 来重复生成向量，其中 x 是“一遍又一遍”。
没有问题。
很棒的包 - 感谢提及，@peter-king。虽然不是主题，但我想知道是否有解决方案可以获得n个随机整数，使它们相加等于给定的整数x。 - Stefan Jelkovich
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接