Используйте parfor ускорять код Монте-Карло

В этом примере показано, как ускорить код Монте-Карло при помощи parfor- циклы. Методы Монте-Карло найдены во многих полях, включая физику, математику, биологию и финансы. Методы Монте-Карло включают выполнение функции много раз со случайным образом распределенными входными параметрами. С Parallel Computing Toolbox можно заменить for- цикл с parfor- цикл, чтобы легко ускорить код.

Этот пример запускает простую стохастическую симуляцию на основе долларового аукциона. Запустите несколько симуляций, чтобы найти рыночную стоимость для однодолларовой банкноты с помощью метода Монте-Карло. В этом примере долларовый аукцион обработан как функция черного ящика, которая производит выходные параметры, которые зависят от вероятностных процессов. Чтобы узнать больше о модели, смотрите Долларовый Аукцион. Чтобы видеть, как ускорить код Монте-Карло в целом, смотрите Использование цикл parfor, чтобы Оценить Рыночную стоимость.

Долларовый аукцион

Долларовый аукцион является игрой небалансовой суммы, сначала введенной Мартином Шубиком в 1 971. В игре, предложении проигрывателей на однодолларовую банкноту. После того, как проигрыватель делает предложение, любой игрок может принять решение сделать предложение выше, чем предыдущий участник торгов. Аукцион заканчивается, когда больше игроков не решает сделать ставку. Лицо, предлагающее самую высокую цену получает однодолларовую банкноту, однако, в отличие от типичного аукциона и лицо, предлагающее самую высокую цену и второе лицо, предлагающее самую высокую цену дают их предложение аукционисту.

Стохастическая модель

Можно смоделировать игры, похожие на долларовый аукцион с помощью стохастической модели. Состояние (текущая ставка и количество активных игроков) может быть смоделировано с помощью процессов Маркова, и поэтому результаты (рыночная стоимость), как могут ожидать, будут иметь четко определенную статистику. Результаты чертятся от условного распределения, и поэтому долларовый аукцион идеален для анализа Монте-Карло. Рыночная стоимость под влиянием следующих факторов:

  • Число игроков, (nPlayers)

  • Проигрыватели действий берут

В этом примере следующий алгоритм определяет то, что проигрыватели действий берут (предложение цены или выбывающий) в зависимости от состояния.

  1. Установите предложение на предыдущее предложение плюс incr.

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

  3. Если никакие ставки не были ранее сделаны, перейдите в 8.

  4. Если предыдущее предложение меньше 1, сгенерируйте случайное число между 0 и 1. Если случайное число меньше dropoutRate, перейдите в 7.

  5. Вычислите сколько денег gain может быть получен, если проигрыватель делает успешную аукционную заявку.

  6. Вычислите сколько денег loss проигрыватель проигрывает, если они - второе лицо, предлагающее самую высокую цену. Если gain больше loss, перейдите в 8.

  7. Проигрыватель выбывает. Удалите проигрыватель из набора проигрывателей, затем перейдите в 9.

  8. Проигрыватель делает ставку.

  9. Если существует 2 или больше проигрывателя, остающиеся, перейдите к шагу 1.

Функция поддержки dollarAuction симулирует долларовый аукцион. Чтобы просмотреть код, смотрите dollarAuction.m. Функция берет три входных параметров: nPlayers, incr, и dropoutRate. Установите каждое из значений.

nPlayers = 20;
incr = 0.05;
dropoutRate = 0.01;

Запустите случайный сценарий путем выполнения dollarAuction функция. Сохраните выходные параметры bids и dropouts.

[bids,dropouts] = dollarAuction(nPlayers,incr,dropoutRate);

В то время как игра продолжается, некоторые проигрыватели делают ставки, и некоторые выбывают. Если предложение превышает 1, проигрыватели заблокированы на "войне предложений", пока только один проигрыватель не остается.

Таблица dropouts содержит две переменные: Player, уникальный номер присвоен каждому проигрывателю; Epoch, раунд предложения цены, когда Player опущенный. Используйте findgroups сгруппировать dropouts.Epoch, и используйте splitapply получить число игроков, кто выбывает в каждом из уникальных раундов в dropouts.Epoch.

[G,epochs] = findgroups(dropouts.Epoch);
numberDropouts = splitapply(@numel,dropouts.Epoch,G);

Первоначально, нет никаких уволенных. Добавьте эту информацию в epochs и numberDropouts преднезаконченным 1 и 0.

epochs = [1;epochs];
numberDropouts = [0;numberDropouts];

Используйте nPlayers и cumsum вычислить число игроков, остающееся от numberDropouts. Вычислите предложения с помощью incr и epochs. Используйте stairs построить предложение против совокупной суммы numberDropouts.

playersRemaining = nPlayers - cumsum(numberDropouts);
stairs(incr*epochs,playersRemaining)
xlabel('Bid')
ylabel('Number of players remaining')

Оцените рыночную стоимость Используя методы Монте-Карло

Можно оценить рыночную стоимость счета со значением origValue при помощи методов Монте-Карло. Здесь, вы производите модель Monte-Carlo и сравниваете скорость с и ни с чем не сравнимый Вычислительный Тулбокс. Определите номер испытаний nTrials используемый к случайным образом демонстрационному результаты.

nTrials = 10000;

Можно произвести возможные исходы путем выполнения функции поддержки dollarAuction многократно. Используйте for- цикл, чтобы произвести nTrials выборки, храня последнее предложение от каждого испытания в B. Каждый раз вы запускаете dollarAuction функция, вы получаете различные результаты. Однако, когда вы много раз запускаете функцию, результаты, к которым вы приводите от всех запусков, будут иметь четко определенную статистику.

Запишите время, потраченное, чтобы вычислить nTrials симуляции. Чтобы уменьшать статистический шум в прошедшее время, повторите этот процесс пять раз, затем займите минимальное прошедшее время.

t = zeros(1,5);
for j = 1:5
    tic
    B = zeros(1,nTrials);
    for i = 1:nTrials
        bids = dollarAuction(nPlayers,incr,dropoutRate);
        B(i) = bids.Bid(end);
    end
    t(j) = toc;
end
forTime = min(t)
forTime = 21.4323

Используйте histogram построить гистограмму финала предлагает B. Используйте xline накладывать график с исходным значением (один доллар) и средней рыночной стоимостью, данной mean.

histogram(B);
origLine = xline(1,'k','LineWidth',3);
marketLine = xline(mean(B),'k--','LineWidth',3);
xlabel('Market value')
ylabel('Frequency')
legend([origLine, marketLine],{'Original value','Market value'},'Location','NorthEast')

С данным алгоритмом и входными параметрами, средняя рыночная стоимость больше исходного значения.

Используйте parfor- цикл, чтобы оценить рыночную стоимость

Можно использовать Parallel Computing Toolbox, чтобы легко ускорить код Монте-Карло. Во-первых, создайте параллельный пул с четырьмя рабочими, использующими 'local' профиль.

p = parpool('local',4);
Starting parallel pool (parpool) using the 'local' profile ...
Connected to the parallel pool (number of workers: 4).

Замените for- цикл с parfor- цикл. Запишите время, потраченное, чтобы вычислить nTrials симуляции. Чтобы уменьшать статистический шум в прошедшее время, повторите, что этот процесс 5 раз затем занимает минимальное прошедшее время.

t = zeros(1,5);
for j = 1:5
    tic
    parfor i = 1:nTrials
        bids = dollarAuction(nPlayers,incr,dropoutRate);
        B(i) = bids.Bid(end);
    end
    t(j) = toc;
end
parforTime = min(t)
parforTime = 5.9174

С четырьмя рабочими результаты indiciate, что код может, работают в три раза быстрее, когда вы используете parfor- цикл.

Приведите к восстанавливаемым результатам со случайными числами в parfor- циклы

Когда вы генерируете случайные числа в parfor- цикл, каждый запуск цикла может привести к различным результатам. Чтобы создать восстанавливаемые результаты, каждая итерация цикла должна иметь детерминированное состояние для генератора случайных чисел. Для получения дополнительной информации смотрите Повторные Случайные числа в циклах parfor.

Функция поддержки dollarAuctionStream берет четвертый аргумент, s. Эта функция поддержки использует заданный поток, чтобы произвести случайные числа. Чтобы просмотреть код, смотрите dollarAuctionStream.m.

Когда вы создаете поток, подпотоки того потока статистически независимы. Для получения дополнительной информации смотрите RandStream. Чтобы гарантировать, что ваш код производит то же распределение результатов каждый раз, создайте поток генератора случайных чисел в каждой итерации цикла, затем установите Substream свойство в индекс цикла. Замените dollarAuction с dollarAuctionStream, затем используйте s запускать dollarAuctionStream на рабочем.

Запишите время, потраченное, чтобы вычислить nTrials симуляции. Чтобы уменьшать статистический шум в прошедшее время, повторите этот процесс пять раз, затем займите минимальное прошедшее время.

t = zeros(1,5);
for j = 1:5
    tic
    parfor i = 1:nTrials
        s = RandStream('Threefry');
        s.Substream = i;
        bids = dollarAuctionStream(nPlayers,incr,dropoutRate,s);
        B(i) = bids.Bid(end);
    end
    t(j) = toc;
end
parforTime = min(t)
parforTime = 8.7355

Масштабируйте от рабочего стола до кластера

Можно масштабировать код от рабочего стола до кластера с большим количеством рабочих. Для получения дополнительной информации об увеличении масштаба от рабочего стола до кластера, смотрите, Масштабируют от Рабочего стола до Кластера.

Используйте delete закрывать существующий параллельный пул.

delete(p);

Вычислите функцию поддержки dollarAuctionStream в parfor- цикл. Запустите тот же parfor- цикл с различными количествами рабочих и запись прошедшие времена. Чтобы уменьшать статистический шум в прошедшее время, запустите parfor- цикл пять раз, затем займите минимальное прошедшее время. Запишите минимальные времена в массиве elapsedTimes. В следующем коде замените MyCluster с именем вашего кластерного профиля.

workers = [1 2 4 8 16 32];
elapsedTimes = zeros(1,numel(workers));

% Create a pool using the 'MyCluster' cluster profile
p = parpool('MyCluster', 32);
Starting parallel pool (parpool) using the 'MyCluster' profile ...
Connected to the parallel pool (number of workers: 32).
for k = 1:numel(workers)    
    t = zeros(1,5);
    for j = 1:5
        tic
        parfor (i = 1:nTrials, workers(k))
            s = RandStream('Threefry');
            s.Substream = i;
            bids = dollarAuctionStream(nPlayers,incr,dropoutRate,s);
            B(i) = bids.Bid(end);
        end
        t(j) = toc;
    end
    
    elapsedTimes(k) = min(t);
end
Analyzing and transferring files to the workers ...done.

Вычислите вычислительное ускорение путем деления elapsedTimes(1) ко временам в elapsedTimes. Исследуйте сильное масштабирование путем графического вывода ускорения против количества рабочих.

speedup = elapsedTimes(1) ./ elapsedTimes;
plot(workers,speedup)
xlabel('Number of workers')
ylabel('Computational speedup')

Вычислительное ускорение увеличивается с количеством рабочих.