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 
% 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 に答える
0

ベクトル化を解除することで実際にこれを解決できると思います。つまり、事前定義された配列と単純な操作のみを使用して、高レベルの呼び出しと高価な操作をすべて削除し、本質的なものにまで落とします。

アルゴリズムのコアは次のようになります。

  1. 重みの合計を決定する

  2. 0 から重みの合計までの n 個の乱数を選択し、並べ替えます。

  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

警告と最適化。

  • 私はこれを広範囲にテストしていません。乱数操作は、適切なランダム動作を実現するのが困難です。多数のモンテカルロ セットに対していくつかのテスト ケースを実行して、動作が期待どおりであることを確認します。特に、off-by-one タイプのエラーに注意してください。

  • プロファイリングし、遅いステップで追加の改善を探します。いくつかの可能性。

    • total変更しても値を維持するMため、再計算する必要はありません。

    • とrandIndexesに対しての最低値と最高値を確認します。もしtotal-randIndexes(end) ixM numel(M) 1 1 numel(M)`.0totalrandIndexes(1) is larger than, then incrementfromto, rather than fromto

于 2012-02-24T18:10:48.850 に答える