Взвешенные случайные числа в MATLAB
Как случайным образом выбрать N чисел из вектора a с весом, присвоенным каждому числу?
Скажем:
a = 1:3; % possible numbers
weight = [0.3 0.1 0.2]; % corresponding weights
В этом случае вероятность забрать 1 должна быть в 3 раза выше, чем для поднятия 2.
Сумма всех весов может быть любой.
Ответы
Ответ 1
R = randsample([1 2 3], N, true, [0.3 0.1 0.2])
randsample входит в панель инструментов статистики
В противном случае вы можете использовать какой-то процесс аналогичный вопрос (хотя и не специфичный для MATLAB). Здесь моя однострочная реализация:
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 get и интервал длины 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]
Ответ 2
amro дает хороший ответ (который я оценил), но он будет очень интенсивным, если вы хотите генерировать множество чисел из большого набора. Это связано с тем, что операция bsxfun может генерировать огромный массив, который затем суммируется. Например, предположим, что у меня был набор из 10000 значений для образца, все с разными весами? Теперь создайте 1000000 номеров из этого образца.
Это потребует некоторой работы, поскольку оно будет генерировать массив 10000x1000000 внутри, с 10 ^ 10 элементами в нем. Это будет логический массив, но при этом необходимо выделить 10 гигабайт памяти.
Лучшим решением является использование 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.
По общему признанию, моя версия занимает 2 строки для записи. Операция индексации должна выполняться во второй строке, так как она использует второй вывод histc. Также обратите внимание, что я использовал способность новой версии matlab с оператором тильды (~) в качестве первого аргумента histc. Это приводит к немедленному сбросу первого аргумента в ведро бит.
Ответ 3
TL; DR
Для максимальной производительности, если вам нужен только один образец, используйте
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) является самым быстрым.
В любом случае randsample, конечно, выглядит неплохим выбором. Также показано, что если алгоритм может быть организован для генерации всех случайных величин заранее, то он будет работать намного лучше (обратите внимание, что на три порядка меньше значений, сгенерированных в случае N=1 за аналогичное время выполнения).
Код доступен здесь.
Ответ 4
У Amro действительно хороший ответ для этой темы. Тем не менее, может потребоваться супер-быстрая реализация для выборки из огромных PDF файлов, где домен может содержать несколько тысяч. Для таких сценариев было бы очень утомительно использовать bsxfun и cumsum очень часто. Мотивированный из ответ Gnovice, было бы разумно реализовать алгоритм колесика рулетки с помощью схемы кодирования длины пробега. Я выполнил бенчмарк с решением 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 больше копируется по сравнению с нечастыми значениями. В конце дня мы используем индекс для взвешенной таблицы образцов, используя равномерное распределение и используем соответствующее значение.
Это интенсивный объем памяти, но при таком подходе даже можно масштабировать до PDF длиной в сотни тысяч. Следовательно, доступ очень быстрый.
