Выберите n взвешенных элементов по индексу из очень большого массива в MATLAB

Question

Предположим, у меня есть очень большая квадратная матрица 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

Но это, кажется, довольно неуклюжий метод, который также требуеточень долгое время из-за цикла.Есть ли более эффективный метод?Возможно, если я каким-то образом векторизую матрицу?

* Изменить Я должен упомянуть, что у меня нет доступа к набору инструментов статистики

Заранее большое спасибо.

Answer 1

randsample ( документы ) - ваш друг здесь.Я бы использовал следующий метод, который преобразует в индексы, а затем обратно в индексы:

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

Возможно, вам придется сделать несколько транспонирований на M, чтобы получить соответствующие измерения.

Answer 2

Я думаю, что на самом деле решил бы это, не векторизация. То есть удалите все высокоуровневые вызовы и дорогостоящие операции и сведите их к основному, используя только предопределенные массивы и простые операции.

Ядром алгоритма будет:

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 than total-randIndexes (конец) , then increment ixM from цифра (M) to 1 , rather than from 1 to цифра (M) `.

Answer 3

% 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 , чтобы сделать полностью векторизованным.