使用 parfor 加速蒙特卡罗代码

美元拍卖

美元拍卖是一种非零和游戏，由马丁·舒比克于 1971 年首次提出。在游戏中，玩家竞标一张一美元的钞票。当一个玩家出价后，其他每个玩家都可以选择做出比之前出价者更高的出价。当没有其他玩家决定出价时，拍卖结束。最高出价者将获得一美元钞票，然而，与典型的拍卖不同的是，最高出价者和第二高出价者都会向拍卖师出价。

随机模型

您可以使用随机模型来仿真类似于美元拍卖的游戏。可以使用马尔可夫进程对状态（当前出价和活跃参与者的数量）进行建模，因此可以预期结果（市场价值）具有明确的统计数据。结果是根据条件分布得出的，因此美元拍卖非常适合蒙特卡罗分析。市场价值受以下因素影响：

在这个示例中，以下算法根据状态确定玩家采取什么行动（出价或退出）。

支持函数 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 的钞票的市场价值。在这里，您生成一个蒙特卡罗模型并比较有和没有 Parallel Computing Toolbox 的速度。设置用于随机抽样结果的试验次数 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 轻松加速您的蒙特卡罗代码。首先，使用 'Processes' 配置文件创建一个具有四个工作单元的并行池。

p = parpool('Processes',4);

Starting parallel pool (parpool) using the 'Processes' profile ...
Connected to parallel pool with 4 workers.

用 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

结果表明，在四个工作单元的情况下，使用 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')

计算速度随着工作单元数量的增加而增加。