蒙特卡洛模拟方法

蒙特卡洛模拟（或方法）是一种将随机性与大数定律应用于解决各类科学与工程问题的技术。其最早有文献记载的应用之一，是二战后 Stanislaw Ulam 和 John von Neumann 将其用于核武器模拟。

在本文中，我想举例说明蒙特卡洛模拟的一些简单用途。我们将从通过投掷飞镖计算 \pi 值的经典示例开始。

估计 pi

假设我们取一块正方形板，并在其中内接一个四分之一圆。然后我们向板上投掷飞镖，并记录每支飞镖是落在四分之一圆内还是圆外。在投掷了许多这样的飞镖之后，我们计算落在圆内的飞镖数量与总投掷数量的比值。

假设我们的飞镖在正方形上均匀分布，根据大数定律，这个比值应当趋近于四分之一圆的面积 A_{circ} 与整个正方形面积 A_{square} 之比。

当正方形边长为 1 时，我们有：

\[\frac{A_{circ}}{A_{square}}=\frac{\pi/4}{1}=\frac{\pi}{4}\]

因此：

\[\pi\approx 4\frac{\text{hits inside}}{\text{total throws}}\]

下面是一个可视化：

单位正方形内半径为 1 的四分之一圆。

更改样本数量（飞镖投掷次数）并点击“Run”以重新生成。代码非常简单——下面是一个稍作整理的版本：

const total = Number(samples.value);
let inside = 0;
for (let i = 0; i < total; i++) {
  const x = Math.random();
  const y = Math.random();
  const inCircle = x * x + y * y <= 1;
  if (inCircle) inside++;
}
estimateValue = 4 * inside / total;

你会注意到这个估计值相当差——即使有 1000 个样本——如果你多次点击“运行”，有些数字会偏离得很远。虽然这种方法确实能估计 \pi，但它并不是一个特别好的估计。我发现需要运行约 100 亿个样本，才能以合理的可靠性将其估计到小数点后 4 位。

一般来说，对于像这样的独立试验，典型的抽样误差与 1/\sqrt{N} 成比例减小，其中 N 是试验次数。这意味着要将误差减半，需要四倍的样本量。

虽然 \pi 的估计看起来有些异想天开，但它是一个重要问题类别的示例，而蒙特卡洛模拟正是应用于此类问题：数值积分。我们的模拟估计了四分之一圆曲线下的面积，这是一个定积分。

许多积分很难解析求解，数值分析领域已经做了大量研究来开发计算积分的方法。蒙特卡洛方法对于高维积分特别有用，在其他数值方法可能变得过于昂贵的情况下尤为如此。

组合模拟——SET 游戏

蒙特卡洛方法的一个常见用途是估计复杂的组合计算。这些计算通常没有解析解，而且由于涉及的数字规模，枚举所有选项是难以处理的。作为一个例子，让我们考虑 SET 游戏)。每张 SET 卡牌有四个属性：

1. 形状数量（1、2 或 3）
2. 颜色（红色、绿色或紫色）
3. 形状类型（椭圆形、菱形或波浪形）
4. 填充方式（空白、条纹或实心）

目标是找出一个“套组”——三张牌，且对每个属性分别而言，要么全部不同，要么全部相同。举个例子，下面这手牌中只有一个套组；看看你能否找出来：

含有一个套组的 SET 手牌

而下一手牌中没有任何套组：

不含套组的 SET 手牌

这里有一个问题：给定一副刚洗好的 SET 牌堆，前 12 张牌中没有任何套组的概率是多少？不使用计算机的话，这个问题很难回答。

计算从 81 张牌的牌堆中发出 12 张手牌的方式数很容易：

\[\binom{81}{12}=70,724,320,184,700\]

但这些手牌中有多少不含套组呢？枚举 70 万亿手 SET 牌并逐一检查会花费相当长的时间，而且没有直接的计数公式可以回答这个问题。可以采用一些巧妙的方法，利用 SET 的对称性和其他数学性质来大幅缩减搜索空间。Donald Knuth 本人就研究过这个问题，并编写了一个简洁的程序（他的程序页面上的 setset-all），找出了 2,284,535,476,080 手这样的牌。因此，我们问题的答案是：

\[\frac{2,284,535,476,080}{70,724,320,184,700}\approx 0.0323\]

从一副完整的 SET 牌堆中随机发出的 12 张手牌中不含任何套组的概率为 3.23%。

我们来看看如何用蒙特卡洛模拟以相对较小的代价回答这个问题，而无需深入了解 SET 的数学性质——那些性质使得 Knuth 式地缩减搜索空间成为可能。我们可以使用以下伪代码：

C = 0
run N times:
  draw a random 12-card hand from a fresh deck
  count sets in the hand
  if no sets:
    C += 1

Estimated probability = C / N

在运行了 1000 万次模拟抽样后，我得到了 0.0323 的结果，这与真实答案非常接近。

蒙特卡洛方法让我们能够以非常简单的方式解决相当复杂的问题。假设我们想对一手 15 张牌回答同样的问题；这会使搜索空间大幅膨胀——选择 15 张牌的组合方式大约比选择 12 张牌的组合方式多 100 倍。但对于蒙特卡洛模拟来说，我们只需调整一个小参数，就能得到一个非常可靠的结果（顺便说一句，大约是 0.00037）。

退休规划

蒙特卡洛模拟无处不在的一个领域是退休投资组合的预测。假设某人准备退休时投资组合总额为 100 万美元；他们希望每年能从投资组合中取出 3 万美元用于生活开支。这样可行吗？

分析这个问题时需要考虑大量因素，但为简单起见，我们只关注两个：投资组合回报率和通货膨胀。我们可以做一个朴素的估算，假设取平均值：假设投资组合年均回报率为 4%，年均通货膨胀率为 2%。

我们将投资组合回报率记为 r=0.04，通货膨胀率记为 q=0.02。那么每年的实际回报率为：

\[r^{\mathrm{real}}=\frac{1+r}{1+q}-1\approx0.0196\]

从 $1,000,000 起步，到年底我们将拥有 $1,019,600，然后取出 $30,000 用于生活开支，最终剩下 $989,600。如果继续这样下去，资金将在约 55 年后耗尽，这意味着一个 65 岁退休的人应该相当安全，对吧？

但这非常简化；仅假设平均回报率是有风险的，因为它们无法很好地代表现实，而且许多因素都存在不确定性。例如，回报的顺序非常重要；一个糟糕的年份（-10%） followed by 一个出色的年份（+18%）仍然会被算作“平均 4%”，但产生的资金显著少于连续两个 +4% 的年份。通货膨胀也是不可预测的，有时与投资组合回报相关；可能会出现高通胀和低/波动回报的糟糕年份，这会对投资组合造成严重破坏。

随着我们添加因素（年度取款的方差、股票、债券、房地产的混合投资组合、预期寿命、意外事件、变化的税法等），依赖单一的平均估计变得越来越危险。这就是为什么蒙特卡洛模拟在该领域非常受欢迎：通过从基于历史数据的合理分布中抽取，蒙特卡洛模拟可以轻松运行一百万个不同的场景并提供估计：例如，在死亡之前资金耗尽的可能性有多大。

这是我运行的一个模拟中的有用图表：

在上图中：

• 虚线显示了之前提到的恒定假设：即年回报率始终为 4%、通货膨胀率始终为 2% 时的情况。

…