0

假设我有一个非常大的方阵 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 循环，这也需要很长时间。有没有更有效的方法？也许如果我以某种方式对矩阵进行矢量化？

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

提前谢谢了。

4

3 回答 3

1

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

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

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

于 2012-02-24T15:52:37.460 回答
0

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

算法的核心是：

  1. 确定权重的总和

  2. 选择 n 个介于 0 和权重之和之间的随机数，对它们进行排序。

  3. 手动实现一个 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

警告和优化。

  • 我没有对此进行广泛的测试。随机数操作很难用于正确的随机行为。在大量蒙特卡罗集上运行一些测试用例，以确保行为符合预期。尤其要注意一对一类型的错误。

  • 分析，然后在任何缓慢的步骤中寻找额外的改进。一些可能性。

    • 在更改时保持total值M，因此您不需要重新计算。

    • randIndexes检查反对0和的最低和最高值total。如果randIndexes(1) is larger thantotal-randIndexes(end) , then incrementixM fromnumel(M) to1 , rather than from1 tonumel(M)`。

于 2012-02-24T18:10:48.850 回答
0 
% 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才能完全矢量化。

于 2012-02-24T16:08:09.633 回答