蒙特卡罗方法求面积：MATLAB随机采样与误差收敛实例

简介：一份面向MATLAB初学者与数值计算学习者的蒙特卡罗方法实例文档，围绕随机抽样估算几何图形阴影面积展开，重点演示均匀网格与随机撒点两种实现思路，并给出可直接运行的MATLAB程序与结果对照。压缩包内共1个doc文件，大小约33KB，以文字讲解配合代码片段的形式呈现，便于随查随用。已有1626人学习下载，说明该案例在入门实践中具有一定参考价值。文档从“豆子撒布”的直观比喻切入，先说明面积近似公式，再分别用固定步长网格和rand随机生成N个点统计落入阴影区域的点数，并讨论样本量、重复试验取均值对精度的影响；附带的程序涵盖向量化操作、循环结构与随机数生成，运行结果分别约为7.2150和7.2500，可与实际面积对照验证。读者可据此理解蒙特卡罗方法的核心流程，掌握两种精度控制策略，并迁移到面积估算、概率模拟等类似问题。

1. 撒豆子求面积：蒙特卡罗为什么能绕开解析解

一片由椭圆和圆交叠出来的阴影区域，解析几何算起来要联立好几条曲线方程，取交点、分段积分，稍不留神就出错。可有人往一个 12×12 的方形布上随手撒了一把豆子，数一数落在阴影里的豆子重量占总重量的比例，再乘布的面积，面积就出来了。这不是玄学，是蒙特卡罗方法最朴素的样子：当精确求解的路径太陡，就用大量独立采样把答案“逼”出来。搞懂它，你就明白为什么随机性反而能换来确定性——这正是蒙特卡罗能在数值积分、物理仿真、金融定价里反复出场的原因，也是我们这份 MATLAB 实例想让你亲手验证的底层机制。

上面那段撒豆子的描述，本质上是用物理过程模拟了一次面积测度。豆子撒得越均匀、数量越多，比例越接近真实面积比。换成代码，就是把豆子当成坐标点，把阴影边界当成判定条件，把比例乘回总面积。整件事没有任何高深数学，只有三个值得深挖的环节：采样怎么产生、命中怎么判断、误差怎么收敛。下面几章就沿着这条线，从概率论里那条期望公式一路拆到 MATLAB 循环里那几个if的嵌套。

2. 面积即概率：把几何测度翻译成命中率

2.1 从均匀分布到期望值的那一步

蒙特卡罗求面积的核心等式只有一行：对于一个完全包含目标区域 $S$ 的参考区域 $B$，在 $B$ 上均匀采样一个点，该点落入 $S$ 的概率恰好是两者面积之比，即 $P = \text{Area}(S)/\text{Area}(B)$。这不需要任何近似，是均匀分布的测度定义直接推出来的结论。于是面积估计变成 $\text{Area}(S) \approx \text{Area}(B) \cdot \hat{P}$，而 $\hat{P}$ 用频率 $k/N$ 去估。

写成一个更工程化的形式就是：

$$\hat{A} = \text{Area}(B) \cdot \frac{1}{N}\sum_{i=1}^{N} \mathbf{1}{(x_i,y_i) \in S}$$

那个指示函数 $\mathbf{1}{\cdot}$ 就是命中判定，取 1 表示点落在阴影内，取 0 表示在外。整段代码的复杂度几乎全花在计算这个指示函数上，占比统计只是一次求和。

值得强调的是，这个估计量是无偏的。因为每次采样独立同分布，$\mathbb{E}[\hat{P}] = P$，所以 $\mathbb{E}[\hat{A}] = \text{Area}(S)$。无偏不代表单次结果准，只代表你重复很多次取平均会收敛到真值——这正是后文方法二重复 100 次取均值的理论依据。

2.2 12×12 方块从哪来：参考区域的选取

实例里反复出现的 $12^2$ 不是随手写的。阴影区域由三个条件圈定：$\frac{x^2}{9}+\frac{y^2}{36}<1$、$\frac{x^2}{36}+y^2<1$、$(x-2)^2+(y+1)^2<9$。第一个是长半轴 6、短半轴 3 的椭圆，第二个是长半轴 6、短半轴 1 的椭圆，第三个是圆心 $(2,-1)$、半径 3 的圆。三者交集必然落在 $x \in [-6,6]$、$y \in [-6,6]$ 之内，所以边长为 12 的正方形是天然的、紧贴边界的参考区域。

参考区域选得越贴合目标，命中率越高，同样样本量下方差越小。如果图省事把参考区域取成 $[-100,100]^2$，命中率会掉到万分之一量级，$N=10000$ 个样本可能只命中个位数，估计值抖得没法看。这是一个容易被忽略的设计选择：蒙特卡罗的精度不只取决于 $N$，还取决于参考区域与目标区域的面积比。常见的实现思路是先用一次粗采样或几何估算把包围盒收紧，再正式撒点。

提示：参考区域必须完整包住目标区域，否则估计值会系统性偏小，且这种偏差不会随 $N$ 增大而消失。

2.3 大数定律：为什么 $N$ 一大结果就稳

估计量的方差是 $\text{Var}(\hat{A}) = \text{Area}(B)^2 \cdot P(1-P)/N$，标准差按 $1/\sqrt{N}$ 衰减。这个衰减速度是整个方法最关键的工程事实：想让误差减半，样本量要翻四倍；想多一位有效数字，样本量得涨一百倍。所以蒙特卡罗从来不是“免费的高精度”，它用可预测的收敛速率换取对复杂边界的免疫能力。

也正因为标准差与 $\sqrt{N}$ 成反比、与 $P(1-P)$ 相关，当 $P$ 接近 0 或 1 时方差会变小。这解释了一个反直觉现象：估算一个几乎占满参考区域的图形，反而比估算一个中等占比的图形更稳。实例里的阴影区域面积约为 7.2，参考面积 144，$P \approx 0.05$，属于低命中率场景，所以 $N=10000$ 时单次估计的波动肉眼可见，必须靠重复取均值压下去。

3. 命中判定：三个不等式如何拼出阴影边界

3.1 嵌套 if 的布尔合取

实例代码里判定一个点是否属于阴影，用的是三层嵌套：

a1 = (x(i)^2)/9 + (y(j)^2)/36; % 第一个椭圆：长半轴6、短半轴3 a2 = (x(i)^2)/36 + y(j)^2; % 第二个椭圆：长半轴6、短半轴1 a3 = (x(i)-2)^2 + (y(j)+1)^2; % 圆：圆心(2,-1)、半径3 if a1 < 1 if a2 < 1 if a3 < 9 k = k + 1; % 三个条件同时满足才计为命中 end end end

这三个不等式分别对应三条边界曲线，嵌套if实现的是逻辑与：只有 $a1<1$ 且 $a2<1$ 且 $a3<9$ 同时成立，点才在阴影内。参数含义很直接：9、36是半轴平方，9是半径平方。把a3 < 9写成开平方形式sqrt((x-2)^2+(y+1)^2) < 3结果一致，但会多一次开方运算，在百万级循环里这种差别会被放大。

这段代码有个值得称道的细节：三个判定按“最容易失败”的顺序排列。第一个椭圆和第二个椭圆的交集已经把大部分点排除了，第三个圆条件放在最内层，实际执行次数最少。这属于分支预测层面的优化，MATLAB 的解释器虽然不如编译型语言敏感，但在 $N$ 上千万时仍然能测出差异。

3.2 从嵌套到向量化：同一个判定的两种写法

嵌套if可读性好，但 MATLAB 的精髓在矩阵运算。等价的向量化写法是把整个采样点集当成数组一次性判定：

N = 1e6; pts = 12 * rand(N, 2) - 6; % N×2 的采样点，范围[-6,6] x = pts(:,1); y = pts(:,2); in1 = (x.^2)/9 + (y.^2)/36 < 1; % 逻辑数组，逐元素比较 in2 = (x.^2)/36 + y.^2 < 1; in3 = (x-2).^2 + (y+1).^2 < 9; hit = in1 & in2 & in3; % 三个逻辑数组按位与 area = 144 * sum(hit) / N; % sum 对 true 计数

&对逻辑数组做逐元素与，sum把true累加成命中数。这种写法省掉了显式循环，执行时间通常比双重for低一到两个数量级，因为循环被下沉到了高度优化的 BLAS 层。代价是内存：$N=10^7$ 时，pts本身就占 160 MB，再加上几个临时逻辑数组，普通机器会开始吃紧。我的经验是单次采样量控制在 $10^6$ 到 $10^7$ 之间，需要更大样本就分批累加，别一次性申请。

注意：向量化写法里x.^2的点号不能漏。漏掉会变成矩阵乘法，维度对不上直接报错，或者更糟——维度凑巧对上时算出完全错误的结果。

3.3 边界点算不算命中：一个不影响结果但影响心态的细节

代码用的是严格小于<，意味着恰好落在边界上的点被判为不命中。理论上边界是零测集，连续分布下恰好命中的概率为零，所以<和<=在数学上等价。但实际运行时rand产生的是离散浮点数，边界点确实可能出现，只不过概率低到可以忽略。

真正需要警惕的是另一类边界问题：当目标区域本身由多个不相交子区域组成，或者边界条件写错方向（比如把<写成>），结果会系统性错误而不是随机波动。排查这类问题的办法是固定随机种子跑一个小样本，把命中点画出来和预期阴影叠在一起看。散点图不会骗人，比盯着数字猜快得多。

4. 两种撒点策略：均匀网格与随机重复的取舍

4.1 方法一：规则网格是一种“去随机化”的蒙特卡罗

方法一的思路是把 $[-6,6]$ 按步长 0.01 切成 1201 个格点，两两组合成约 144 万个采样点，逐个判定后统计命中数：

x = -6:0.01:6; y = x; s = size(x); zs = s(1,2)^2; % 总点数为 1201^2 k = 0; for i = 1:s(1,2) for j = 1:s(1,2) a1 = (x(i)^2)/9 + (y(j)^2)/36; a2 = (x(i)^2)/36 + y(j)^2; a3 = (x(i)-2)^2 + (y(j)+1)^2; if a1 < 1 && a2 < 1 && a3 < 9 k = k + 1; end end end mj = (12^2) * k / zs; % 运行结果 7.2150

这里的“随机性”被彻底去掉了，采样点是确定性的网格。它严格来说不是蒙特卡罗，而是黎曼和的一种变体——用规则网格上的示性函数均值去近似积分。好处是结果完全可复现，没有随机波动；坏处是精度被网格步长锁死，而且存在系统性偏差：网格点与边界曲线的相对位置固定，边界附近的取舍总是偏同一方向。

4.2 方法二：重复实验把随机性“平均掉”

方法二走标准蒙特卡罗路线，每次生成 $N=10000$ 个随机点，重复 100 轮，最后对 100 个面积估计取均值：

N = 10000; n = 100; for j = 1:n k = 0; for i = 1:N a = 12*rand(1,2) - 6; % 生成[-6,6]区间内的随机点 x = a(1,1); y = a(1,2); a1 = (x^2)/9 + (y^2)/36; a2 = (x^2)/36 + y^2; a3 = (x-2)^2 + (y+1)^2; if a1 < 1 && a2 < 1 && a3 < 9 k = k + 1; end end m(j) = (12^2) * k / N; % 单轮估计 end mj = mean(m); % 100 轮均值，运行结果 7.2500

12*rand(1,2)-6把 $[0,1)$ 上的均匀随机数线性映射到 $[-6,6)$，这是rand最常用的变换。单轮 $N=10000$ 的标准差约为 $144 \times \sqrt{0.05 \times 0.95 / 10000} \approx 0.31$，单次结果在 7.2 附近上下跳是正常的。取 100 轮均值后，标准差降到约 0.031，稳定性提升一个数量级。

4.3 两种策略的对比与适用边界

二维问题里规则网格其实很有竞争力，因为 $O(h)$ 在低维下衰减快。但维度一高，网格点数量指数爆炸，$d$ 维、每维 $m$ 个点就是 $m^d$，而随机采样的误差衰减速率 $1/\sqrt{N}$ 与维度无关——这才是蒙特卡罗真正的杀手锏，也是它在高维积分里不可替代的原因。实例只做二维演示，但理解了这一点，你就明白为什么它值得被单独拿出来讲。

提示：方法二里rand的随机种子默认随启动变化，想复现某次结果需要在开头调用rng(固定整数)。做对比实验时忘了这一步，两次运行结果对不上，很容易误判成代码有 bug。

5. 常见问题排查：五个让结果对不上的坑

5.1 现象：估计值稳定地偏小约 10%

原因通常是参考区域没有完整包住目标区域，或者边界条件的方向写反。比如把a3 < 9误写成a3 > 9，命中集合直接变成圆外区域，结果会莫名其妙地偏大且形状完全不对。解决方式是把三个条件分别单独画图验证：对每个条件生成一批点，只画满足该条件的散点，确认单条边界正确后再组合。分而治之比一次性调试四个不等式快得多。

5.2 现象：重复跑结果忽大忽小，差距超过预期

这是单轮样本量太小、又没做重复平均的典型症状。$N=1000$ 时单轮标准差约 0.98，结果在 6 到 8 之间跳完全正常。解决办法有两个：把单轮 $N$ 提到 $10^4$ 以上，或者保持小 $N$ 但重复几百轮取均值。两者总计算量相当，但后者还能顺便估计出标准差，给结果配一个置信区间，工程上更有说服力。

5.3 现象：向量化版本报错“矩阵维度必须一致”

多半是.^写成了^，或者两个数组的形状不匹配。rand(N,2)返回 $N \times 2$，pts(:,1)是 $N \times 1$ 列向量，两个列向量做.^2没问题；但如果某处误用了行向量，MATLAB 会尝试矩阵乘法并给出维度错误。排查时用whos看每个变量的尺寸，比读报错信息定位更快。

5.4 现象：大样本运行时内存耗尽

一次性生成 $10^8$ 个点需要约 1.6 GB 仅存坐标，加上中间逻辑数组很容易撑爆。解决办法是分批：外层循环每次处理 $10^6$ 个点，累加命中数，最后统一除以总样本量。分批不改变估计的无偏性，只是把内存压力摊平。写成分批版本后，普通笔记本也能跑到 $10^9$ 量级样本。

5.5 现象：结果和解析解对不上，但代码看起来没错

先怀疑参考面积乘错了。实例里参考区域是边长 12 的正方形，面积 144，代码里写的是12^2。如果误写成12*12之外的数，或者把半轴当成长轴，结果会整体缩放。第二个常见原因是把k/N里的 $N$ 写成了循环变量而不是总点数。建议在计算面积那一行之前加一句断言，检查k不超过N、N等于预期值，这两条能拦下大部分低级错误。

6. 精度、维度与设计取舍：什么时候该信蒙特卡罗

蒙特卡罗的诱惑在于实现简单，陷阱在于误差不可控地累积。理解它的设计取舍，比会写那几行循环重要得多。第一层取舍是精度与成本的兑换关系：误差按 $1/\sqrt{N}$ 衰减，意味着从两位有效数字提升到三位，计算量要涨一百倍。所以但凡存在解析解或者高精度数值积分可用，就不该无脑上蒙特卡罗——它是给“解析路径走不通”的问题准备的备用通道，不是通用加速器。

第二层取舍是维度。二维下规则网格的 $O(h)$ 收敛其实优于随机采样的 $O(1/\sqrt{N})$，实例里方法一和方法二精度接近，正是这个原因。但维度升到十几维，网格法彻底失效，随机采样的收敛速率却与维度无关，这时候蒙特卡罗的劣势瞬间变成优势。判断标准很清晰：低维光滑边界优先网格或求积公式，高维或边界病态才轮到随机采样登场。

第三层取舍在方差缩减。朴素蒙特卡罗的方差只靠堆样本压，而重要性采样、分层采样、对偶变量这些技巧能在同样样本量下把方差砍掉一大截。实例里“重复 100 次取均值”其实就是最朴素的一种方差缩减——它等价于把总样本量提到 $10^6$，只是顺手给出了误差估计。再往深走一步，可以在采样时按区域面积比例分配样本，让命中率从 5% 提到 30% 以上，同样计算量下标准差能降一半。

验证方法也很直接：把估计值和解析解对比只是第一步，更重要的是看估计值随 $N$ 增长的收敛曲线。健康的结果应该在对数坐标下呈现斜率约 $-1/2$ 的直线，如果斜率明显偏离，说明采样有偏或者边界判定有系统性错误。我一般会跑三组不同 $N$ 的试验，画出这条线，比单看一个数字可靠得多。

注意：固定随机种子只保证复现，不保证正确。种子固定后结果依然可能系统性偏离真值，验证收敛性才是判断对错的唯一手段。

从那以后我每次写蒙特卡罗的代码，都强制先跑一个小样本把命中点画出来看一眼，再放大样本量做正式计算。这一步多花两分钟，能省掉后面半小时盯着数字怀疑人生的时间。这份 MATLAB 实例把两种撒点策略和完整代码都摆在一起，正好适合拿来跑通这条从采样到收敛的验证链路，动手改改参数、画画散点，比读十篇原理文章都管用。希望帮到你。

本文还有配套的精品资源，点击获取