Взвешенные случайные числа в MATLAB

Как случайным образом выбрать N чисел из вектора a с весом, присвоенным каждому числу?

Предположим:

 a = 1:3; % possible numbers weight = [0.3 0.1 0.2]; % corresponding weights 

В этом случае вероятность подобрать 1 должна быть в 3 раза выше, чем подобрать 2.

Сумма всех весов может быть любой.

     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]

    amro дает хороший ответ (который я оценил), но он будет очень интенсивным, если вы хотите генерировать множество чисел из большого набора. Это связано с тем, что операция bsxfun может генерировать огромный массив, который затем суммируется. Например, предположим, что у меня был набор из 10000 значений для образца, все с разными весами? Теперь создайте 1000000 номеров из этого образца.

    Это потребует некоторой работы, так как оно будет генерировать массив 10000×1000000 внутри, с 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. Это приводит к немедленному сбросу первого аргумента в ведро бит.

    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) , 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) , 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) , 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) , 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 версии bsxfun ) является самым быстрым.

    В любом случае, randsample безусловно, выглядит неплохим выбором. Также показано, что если алгоритм может быть организован для генерации всех случайных величин заранее, то он будет работать намного лучше (обратите внимание, что на три порядка меньше значений, сгенерированных в случае N=1 в аналогичное время выполнения).

    Код доступен здесь .

    У Amro есть действительно хороший ответ для этой темы. Тем не менее, может потребоваться супер-быстрая реализация для выборки из огромных PDF-файлов, где домен может содержать несколько тысяч. Для таких сценариев было бы очень утомительно использовать bsxfun и cumsum очень часто. Мотивированный ответом Гновице , было бы разумно реализовать алгоритм колесика рулетки с помощью схемы кодирования длины прогона. Я выполнил бенчмарк с решением 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-страниц сотен тысяч. Следовательно, доступ очень быстрый.

    Interesting Posts

    Импорт файлов Excel в R, xlsx или xls

    AngularJS Upgrade (от 1.5 до 1.6.1.7) Делает привязку видимости области действия неопределенной

    Восстановление Windows XP Active Desktop для пользователей с ограниченным доступом

    Что такое сырой тип и почему мы не должны его использовать?

    Изменение отображения по умолчанию для строки «не анализировано» в Elasticsearch

    Mac: Как записывать телефонный звонок VoIP (микрофон и выход одновременно)?

    Как вредоносная программа влияет на маршрутизатор?

    Какова цель анонимных блоков {} в языках стиля C?

    Как начать сеанс консоли с моего маршрутизатора Cisco 5 долларов США через порт RS-232?

    Является ли это совместимым расширением компилятора для обработки стандартных библиотечных функций, отличных от constexpr, как constexpr?

    Как разбить строку в оболочке и получить последнее поле

    C / C ++: арифметика указателей

    Когда stream Java жив?

    Использование примеров и примеров шаблона декоратора GoF для IO

    Класс Android Parcelable с ArrayList

    Давайте будем гением компьютера.