从 MATLAB 中的一个非常大的数组中按索引选择 n 个加权元素答案

【问题标题】：Select n weighted elements by index from a very large array in MATLAB从 MATLAB 中的一个非常大的数组中按索引选择 n 个加权元素
【发布时间】：2012-02-24 14:54:55
【问题描述】：

假设我有一个非常大的方阵 M(i, j)，这样矩阵中的每个元素都代表了在加权随机选择中选择该元素的概率。我需要从矩阵中采样 n 个元素（通过 (i, j) 索引）并进行替换。权重会在主循环的每次迭代中发生变化。

目前，我正在使用类似以下内容：

for m = 1:M_size
    xMean(m) = mean(M(:, m));
end

[~, j_list] = histc(rand(n, 1), cumsum([0; xMean'./sum(xMean)']));
for c = 1:n
    [~, i_list(c)] = ...
      histc(rand(1, 1), cumsum([0;, M(:, j_list(c))./sum(M(:, j_list(c)))]));
end

但这似乎是一个相当笨重的方法，由于 for 循环，这也需要很长时间。有没有更有效的方法？也许如果我以某种方式对矩阵进行矢量化？

*编辑我应该提到我无权访问统计工具箱

非常感谢。

【问题讨论】：

在尝试对矩阵进行矢量化以便只调用 histc 和 cumsum 一次之后，似乎这样做实际上需要更长的时间......
你知道M矩阵的总和吗（没有实际求和）？
我不知道没有求和的总和，因为概率是任意估价的。这是因为它被用作强化学习算法的一部分，它只会增加矩阵中少数元素的值。我想我可以简单地计算一次总和，然后在每次传递时添加新值。

标签： matlab matrix random-sample weighted

【解决方案1】：

randsample (docs) 是你的朋友。我将使用以下方法转换为索引然后返回下标：

selected_indexes = randsample(1:numel(M), n, true, M(:));
[sub_i, sub_j] = ind2sub(size(M), selected_indexes);

您可能需要对 M 进行一些转置才能获得适当的尺寸。

【讨论】：

感谢您的回复，但不幸的是，randsample 是统计工具箱的一部分，我无权访问（我的机构是小气鬼）。
哦，对不起，我像你一样读了它。我应该更加小心。

【解决方案2】：

% M is ixj
xMean = transpose(mean(M,1));
%xMean is jx1, so i hope n == j
[~, j_list] = histc(rand(n, 1), cumsum([0; xMean./sum(xMean)]));
% j_list is not used? but is j x 1
cumsumvals = cumsum([zeros(1,jj);, M(:,j_list(1:n))./kron(sum(M(:,j_list(1:n))),ones(ii,1))],1),1)
% cumsumvals is i+1 x j, so looks like it should work
% but histc won't work with a matrix valued edge parameter
% you'll need to look into hist3 for that
for c = 1:n
    [~, i_list(c)] = ...
      histc(rand(1, 1), cumsumvals(:,c));
end

所以它更接近，但你需要hist3 才能完全矢量化。

【讨论】：

【解决方案3】：

我想我实际上可以通过取消矢量化来解决这个问题。也就是说，删除所有高级调用和昂贵的操作，并将其剥离到基本要素，仅使用预定义的数组和简单的操作。

算法的核心是：

确定权重之和
选择n个介于0和权重之和之间的随机数，对它们进行排序。
手动实现一个 cumsum 循环。但是，不是存储所有累积和，而是存储累积和从小于当前随机数跳到大于当前随机数的索引。

在代码中（带有一点计时装置），看起来像这样：

tic
for ixTiming = 1:1000

    M = abs(randn(50));
    M_size = size(M, 2);
    n = 8;
    total = sum(M(:));

    randIndexes = sort(rand(n,1) * total);

    list = zeros(n,1);
    ixM = 1;
    ixNextList = 1;
    curSum = 0;
    while ixNextList<=n  && ixM<numel(M)
        while curSum<randIndexes(ixNextList) && ixM<=numel(M)
            curSum = curSum+M(ixM);
            ixM = ixM + 1;
        end
        list(ixNextList) = ixM;
        ixNextList = ixNextList+1;
    end
    [i_list, j_list] = ind2sub(size(M),list);

end
toc; %0.216 sec. on my computer

将此与原始问题中的代码进行比较：

tic
for ixTiming = 1:1000
    M = abs(randn(50));
    M_size = size(M, 2);
    n = 8;

    for m = 1:M_size
        xMean(m) = mean(M(:, m));
    end

    [~, j_list] = histc(rand(n, 1), cumsum([0; xMean'./sum(xMean)']));
    for c = 1:n
        [~, i_list(c)] = ...
            histc(rand(1, 1), cumsum([0;, M(:, j_list(c))./sum(M(:, j_list(c)))]));
    end
end
toc;  %1.10 sec on my computer

注意事项和优化。

我没有对此进行广泛的测试。随机数操作很难用于正确的随机行为。在大量蒙特卡罗集上运行一些测试用例，以确保行为符合预期。尤其要注意一对一类型的错误。
配置文件，然后在任何缓慢的步骤中寻找其他改进。一些可能性。
- 在更改M 时保持total 值，因此您无需重新计算。
- 对照0 和total 检查randIndexes 的最小值和最大值。如果randIndexes(1) is larger thantotal-randIndexes(end), then incrementixMfromnumel(M)to1, rather than from1tonumel(M)`.

【讨论】：