MATLAB中的加权随机数
Posted
技术标签:
【中文标题】MATLAB中的加权随机数【英文标题】：Weighted random numbers in MATLAB 【发布时间】：2011-02-27 23:34:03 【问题描述】：如何从向量a中随机选取N个数字，并为每个数字分配权重？
假设：
a = 1:3; % possible numbers
weight = [0.3 0.1 0.2]; % corresponding weights
在这种情况下，捡到 1 的概率应该是捡到 2 的 3 倍。
所有权重的总和可以是任何值。
【问题讨论】：
权重值是否介于 0 和 1 之间？ 没有必要。但可以归一化。 结果应该是 a (fe [2 1 3]) 中数字的排列，还是允许 a 中的数字也有重复 (fe [2 2 1]) ? 当然允许重复。 【参考方案1】：R = randsample([1 2 3], N, true, [0.3 0.1 0.2])
randsample 包含在统计工具箱中
否则你可以使用某种roulette-wheel selection 进程。请参阅此similar question（尽管不是特定于 MATLAB）。这是我的单行实现：
a = 1:3; %# possible numbers
w = [0.3 0.1 0.2]; %# corresponding weights
N = 10; %# how many numbers to generate
R = a( sum( bsxfun(@ge, rand(N,1), cumsum(w./sum(w))), 2) + 1 )
说明：
考虑区间 [0,1]。我们为列表中的每个元素（1:3）分配一个与每个元素的权重成比例的长度子区间；因此1 得到长度间隔0.3/(0.3+0.1+0.2)，其他相同。
现在如果我们生成一个在 [0,1] 上均匀分布的随机数，那么 [0,1] 中的任何数字都有相同的概率被选中，因此子区间的长度决定了随机数的概率落在每个区间的数字。
这符合我上面所做的：选择一个数字 X~U[0,1]（更像N 数字），然后以矢量化的方式找到它落入哪个区间..
您可以通过生成足够大的序列N=1000来检查上述两种技术的结果：
>> tabulate( R )
Value Count Percent
1 511 51.10%
2 160 16.00%
3 329 32.90%
或多或少匹配归一化权重w./sum(w)[0.5 0.16667 0.33333]
【讨论】：
非常感谢您提供如此详细的回答。 一个小评论...randsample 不支持无放回加权随机抽样。
@BajajG OP 特别想要更换采样。还有另一个函数datasample 支持无替换加权采样（根据文档，使用 Wong 和 Easton 的算法）【参考方案2】：
amro 给出了一个很好的答案（我对此进行了评价），但如果您希望从一个大集合中生成许多数字，它将非常密集。这是因为 bsxfun 操作可以生成一个巨大的数组，然后对其求和。例如，假设我有一组 10000 个值可供抽样，所有值都具有不同的权重？现在，从该样本中生成 1000000 个数字。
这需要做一些工作，因为它将在内部生成一个 10000x1000000 数组，其中包含 10^10 个元素。这将是一个逻辑数组，但即便如此，也必须分配 10 GB 的内存。
更好的解决方案是使用 histc。于是……
a = 1:3
w = [.3 .1 .2];
N = 10;
[~,R] = histc(rand(1,N),cumsum([0;w(:)./sum(w)]));
R = a(R)
R =
1 1 1 2 2 1 3 1 1 1
但是，对于我上面建议的大问题，它很快。
a = 1:10000;
w = rand(1,10000);
N = 1000000;
tic
[~,R] = histc(rand(1,N),cumsum([0;w(:)./sum(w)]));
R = a(R);
toc
Elapsed time is 0.120879 seconds.
诚然，我的版本需要 2 行来编写。索引操作必须发生在第二行，因为它使用 histc 的第二个输出。另请注意，我使用了新的 matlab 版本的功能，将波浪号 (~) 运算符作为 histc 的第一个参数。这会导致第一个参数立即转储到位桶中。
【讨论】：
@woodcihps，谢谢你的好解决方案。顺便说一句，Amro 建议的 RANDSAMPLE 函数也使用 histc 方法。 @woodschips：感谢改进，看看 RANDSAMPLE 它也在使用 HISTC【参考方案3】：TL;DR
为了获得最佳性能，如果您只需要一个样本，请使用
R = a( sum( (rand(1) >= cumsum(w./sum(w)))) + 1 );
如果您需要多个样本，请使用
[~, R] = histc(rand(N,1),cumsum([0;w(:)./sum(w)]));
避免randsample。预先生成多个样本比生成单个值***个数量级。
性能指标
由于这显示在我的 Google 搜索顶部附近，我只想添加一些性能指标，以表明正确的解决方案在很大程度上取决于 N 的值和应用程序的要求。此外，更改应用程序的设计可以显着提高性能。
对于大的N，或者实际上是N > 1：
a = 1:3; % possible numbers
w = [0.3 0.1 0.2]; % corresponding weights
N = 100000000; % number of values to generate
w_normalized = w / sum(w) % normalised weights, for indication
fprintf('randsample:\n');
tic
R = randsample(a, N, true, w);
toc
tabulate(R)
fprintf('bsxfun:\n');
tic
R = a( sum( bsxfun(@ge, rand(N,1), cumsum(w./sum(w))), 2) + 1 );
toc
tabulate(R)
fprintf('histc:\n');
tic
[~, R] = histc(rand(N,1),cumsum([0;w(:)./sum(w)]));
toc
tabulate(R)
结果：
w_normalized =
0.5000 0.1667 0.3333
randsample:
Elapsed time is 2.976893 seconds.
Value Count Percent
1 49997864 50.00%
2 16670394 16.67%
3 33331742 33.33%
bsxfun:
Elapsed time is 2.712315 seconds.
Value Count Percent
1 49996820 50.00%
2 16665005 16.67%
3 33338175 33.34%
histc:
Elapsed time is 2.078809 seconds.
Value Count Percent
1 50004044 50.00%
2 16665508 16.67%
3 33330448 33.33%
在这种情况下，histc 最快
但是，在可能无法预先生成所有 N 个值的情况下，可能是因为每次迭代都会更新权重，即N=1：
a = 1:3; % possible numbers
w = [0.3 0.1 0.2]; % corresponding weights
I = 100000; % number of values to generate
w_normalized = w / sum(w) % normalised weights, for indication
R=zeros(N,1);
fprintf('randsample:\n');
tic
for i=1:I
R(i) = randsample(a, 1, true, w);
end
toc
tabulate(R)
fprintf('cumsum:\n');
tic
for i=1:I
R(i) = a( sum( (rand(1) >= cumsum(w./sum(w)))) + 1 );
end
toc
tabulate(R)
fprintf('histc:\n');
tic
for i=1:I
[~, R(i)] = histc(rand(1),cumsum([0;w(:)./sum(w)]));
end
toc
tabulate(R)
结果：
0.5000 0.1667 0.3333
randsample:
Elapsed time is 3.526473 seconds.
Value Count Percent
1 50437 50.44%
2 16149 16.15%
3 33414 33.41%
cumsum:
Elapsed time is 0.473207 seconds.
Value Count Percent
1 50018 50.02%
2 16748 16.75%
3 33234 33.23%
histc:
Elapsed time is 1.046981 seconds.
Value Count Percent
1 50134 50.13%
2 16684 16.68%
3 33182 33.18%
在这种情况下，自定义cumsum 方法（基于bsxfun 版本）最快。
无论如何，randsample 看起来确实是个糟糕的选择。它还表明，如果可以安排算法预先生成所有随机变量，那么它的性能会好多更好（请注意，在N=1 案例中生成的值要少三个数量级类似的执行时间）。
代码可用here。
【讨论】：
感谢您的宝贵意见。很有趣。【参考方案4】：Amro 对这个话题有一个非常好的答案。然而，人们可能想要一种超快速的实现来从域可能包含数千个的巨大 PDF 中进行采样。对于这样的场景，经常使用 bsxfun 和 cumsum 可能会很乏味。受Gnovice's answer 的启发，使用游程编码模式实现轮盘赌算法是有意义的。我使用 Amro 的解决方案和新代码进行了基准测试：
%% Toy example: generate random numbers from an arbitrary PDF
a = 1:3; %# domain of PDF
w = [0.3 0.1 0.2]; %# Probability Values (Weights)
N = 10000; %# Number of random generations
%Generate using roulette wheel + run length encoding
factor = 1 / min(w); %Compute min factor to assign 1 bin to min(PDF)
intW = int32(w * factor); %Get replicator indexes for run length encoding
idxArr = zeros(1,sum(intW)); %Create index access array
idxArr([1 cumsum(intW(1:end-1))+1]) = 1;%Tag sample change indexes
sampTable = a(cumsum(idxArr)); %Create lookup table filled with samples
len = size(sampTable,2);
tic;
R = sampTable( uint32(randi([1 len],N,1)) );
toc;
tabulate(R);
上面代码的一些评估是针对 PDF 域包含巨大长度的非常大的数据。
a ~ 15000, n = 10000
Without table: Elapsed time is 0.006203 seconds.
With table: Elapsed time is 0.003308 seconds.
ByteSize(sampTable) 796.23 kb
a ~ 15000, n = 100000
Without table: Elapsed time is 0.003510 seconds.
With table: Elapsed time is 0.002823 seconds.
a ~ 35000, n = 10000
Without table: Elapsed time is 0.226990 seconds.
With table: Elapsed time is 0.001328 seconds.
ByteSize(sampTable) 2.79 Mb
a ~ 35000 n = 100000
Without table: Elapsed time is 2.784713 seconds.
With table: Elapsed time is 0.003452 seconds.
a ~ 35000 n = 1000000
Without table: bsxfun: out of memory
With table : Elapsed time is 0.021093 seconds.
这个想法是创建一个运行长度编码表，与非频繁值相比，PDF 的频繁值被复制得更多。在一天结束时，我们对加权样本表的索引进行抽样，使用均匀分布，并使用相应的值。
它是内存密集型的，但使用这种方法甚至可以扩展到数十万的 PDF 长度。因此访问速度非常快。
【讨论】：
以上是关于MATLAB中的加权随机数的主要内容，如果未能解决你的问题，请参考以下文章