MATLAB中的加权随机数

R = randsample([1 2 3], N, true, [0.3 0.1 0.2])

randsample

否则你可以用一些 roulette-wheel selection 过程。看到这个了吗 similar question

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]

amro给出了一个很好的答案(我对此进行了评价),但是如果您希望从一个大的集合中生成许多数字,那么它将是高度密集的。这是因为bsxfun操作可以生成一个巨大的数组,然后对其求和。例如,假设我有一组10000个值要从中采样,所有这些值都具有不同的权重?现在,从这个样本中生成1000000个数字。

这需要做一些工作,因为它将在内部生成一个10000x1000数组,其中包含10^10个元素。它将是一个逻辑阵列,但即便如此,也必须分配10gb的ram。

更好的解决方法是使用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.

诚然,我的版本需要两行字。索引操作必须发生在第二行,因为它使用histc的第二个输出。还要注意的是,我使用了新的matlab版本的功能,将tilde(~)操作符作为histc的第一个参数。这将导致第一个参数立即转储到位桶中。

热释光;博士

为了获得最佳性能,如果您只需要一个样本,请使用

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

无论如何, 从各方面来看,这的确是个糟糕的选择。它还表明,如果一个算法可以预先生成所有随机变量,那么它将执行 许多的 类似的执行时间)。

here .

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长度。因此访问速度非常快。