使用Matlab生成给定概率的随机数

Question

我想生成一个给定概率的随机数，但不确定如何实现：

我需要一个介于1和3之间的数字。

num = ceil(rand*3);

但我需要不同的值具有生成不同概率的能力。

0.5 chance of 1
0.1 chance of 2
0.4 chance of 3

我相信这很简单，但我想不出该如何做。

- Eamonn McEvoy

---

这个回答解决了你的问题吗？在Matlab中从预先指定的概率质量函数中绘制随机数 - SecretAgentMan

Answer 1

简单的解决方案是使用均匀分布生成一个数字（使用rand），并稍微操纵一下：

r = rand;
prob = [0.5, 0.1, 0.4];
x = sum(r >= cumsum([0, prob]));

或者说成一行代码：

x = sum(rand >= cumsum([0, 0.5, 0.1, 0.4]));

说明

这里的r是一个在0到1之间均匀分布的随机数。为了生成一个介于1和3之间的整数，技巧是将[0, 1]范围划分为3个段，每个段的长度与其对应的概率成比例。在您的情况下，您将拥有：

• 段[0, 0.5)，对应数字1。
• 段[0.5, 0.6)，对应数字2。
• 段[0.6, 1]，对应数字3。

r落入任何一段的概率与你想要的每个数字的概率成比例。sum(r >= cumsum([0, prob]))只是一种将整数映射到其中一个段的花哨方式。

扩展

如果您希望创建随机数的向量/矩阵，则可以使用循环或arrayfun：

r = rand(3); % # Any size you want
x = arrayfun(@(z)sum(z >= cumsum([0, prob])), r);

当然，也有一种向量化的解决方案，只是我太懒了，不想写。

- Eitan T

---

谢谢你的帮助，当概率为[0,0,1]时，有什么我可以做的吗？在这种情况下，我需要答案是3，但一直得到1。 - Eamonn McEvoy

PS：您忘记在此行中添加0 -> x = sum(rand >= cumsum([0.5, 0.1, 0.4, 0])); - Eamonn McEvoy

向量化解决方案：sum(bsxfun(@ge, r, cumsum([0, prob]), 2))，其中 r 是列向量，prob 是行向量。 - Oleg

@OlegKomarov 谢谢您提供的向量化解决方案 ;) - Eitan T

@yashar 谢谢，已修复。 - Eitan T

Answer 2

到目前为止的答案都是正确的，但对于大量输入来说速度较慢：O(m*n)，其中n是值的数量，m是随机样本的数量。这里有一个O(m*log(n))版本，利用了cumsum结果的单调性和histc中使用的二进制搜索。

% assume n = numel(prob) is large and sum(prob) == 1
r = rand(m,1);
[~,x] = histc(r,cumsum([0,prob]));

- Alec Jacobson

---

小提示，取决于histc的实现方式，这可能是O(n+m*log(n))。但我希望由于第一个输出没有被使用，这不是情况。 - Alec Jacobson

有一种更好的O(n+m)解决方案，使用别名方法。我已经在sample_discrete.m函数中实现了它。 - Alec Jacobson

链接已经失效，供您参考。 - intrepid_em

Answer 3

使用统计与机器学习工具箱中的randsample函数，您可以生成具有指定概率质量函数（pmf）的随机数：

pmf = [0.5, 0.1, 0.4];
population = 1:3;
sample_size = 1;

random_number = randsample(population,sample_size,true,pmf);

我认为这是最简单的方法。

- mamaj

Answer 4

>> c = cumsum([0.5, 0.1, 0.4]);
>> r = rand(1e5, 1);
>> x = arrayfun(@(x) find(x <= c, 1, 'first'), r);
>> h = hist(x, 1:3)

h =

       49953       10047       40000

x 分布按需求。

- Serg

---

@EitanT，我认为sum()并不比find(..., 'first')快。而且，没有必要加零。请进行测试。一般情况下，我只会添加：assert(c(end) == 1);。 - Serg

现在我想起来，我的评论不合适。对不起。 - Eitan T

Answer 5

稍微通用一些的解决方案是：

r=rand;
prob=[.5,.1,.4];
prob=cumsum(prob);
value=[1,2,3];    %values corresponding to the probabilities
ind=find(r<=prob,1,'first');
x=value(ind)

- user1860611

Answer 6

一个使用rand、cumsum和min的向量解决方案。

r = rand(10,1);
p = [0.5 0.1 0.4];
[~, ind] = min(r >= cumsum(p), [], 2)

• 使用rand从0..1中随机抽样r。在这种情况下，我将我的数据放入列向量中。
• 将每个输出索引的概率放入p中。
• r >= cumsum(p)比较r和p的累积概率的每个组合。在这种情况下，结果是一个二维矩阵，其中每行以一系列1开头，以一系列0结尾。第一个0表示随机选择的p元素。
• 对所有行执行min并返回第一个0的列索引。min的第三个输入定义要计算最小值的维度。

如果您想将此扩展到r的n个维度：更改p的形状，使其延伸到比r多一个维度，并将该维度作为min的第三个输入。

r = rand(3, 5, 7);
p = []; 
p(1,1,1,:) = [0.5 0.1 0.4];
[~, ind] = min(r >= cumsum(p), [], 4)

- aimor

Answer 7

当概率是像这样的好数字时，可以进行非常简单和高效的选择。我们重复种群元素，以便均匀选择产生所需的概率分布。在这种情况下，我们创建一个由10个元素组成的种群，其中5次为1（被选中的概率为0.5），等等。

p = [1,1,1,1,1,2,3,3,3,3];
x = p(randi(numel(p));

randi 接受第二个输入参数来确定输出的大小（默认为1），因此从此分布生成多个值很简单。

- Cris Luengo