简介围绕蒙特卡罗方法在MATLAB中的实际应用这份文档以计算几何图形阴影部分面积为主线面向正在学习数值计算、概率统计或准备数学建模的读者提供可直接参照的完整实现思路。压缩包内含1个doc文件约33KB体量轻巧但内容集中涵盖方法讲解、程序代码与运行结果三部分。文档对比了两种求解路径一是将坐标轴按固定步长均匀划分、逐点判定是否落入阴影区域二是用rand函数在(-6,6)区间随机撒点并重复100次取均值分别得到7.2150与7.2500的近似面积同时说明样本量N与精度之间的关系。读者可从中获得向量化写法、嵌套条件判断、随机数生成等MATLAB编程细节以及均匀撒点与随机撒点两种策略的误差对比分析便于迁移到物理模拟、金融模型等解析解难求的问题上。目前已有1626人学习。1. 从一袋豆子说起三个区域交叠出的阴影面积群里有人丢了一张图过来问阴影部分面积怎么算。图形由两个椭圆和一个圆交叠而成边界全是二次曲线硬要用解析法求得先算交点、再分段积分写完半页纸还容易漏掉边界情形。这类问题其实有一个很朴素的办法把豆子撒在一块边长 12 的正方形布上数一数落在阴影里的豆子占比乘回 144面积就出来了。这就是蒙特卡罗方法最直观的样子而在 MATLAB 里把它变成代码核心循环不过十几行。这篇笔记拆的就是这份蒙特卡罗 MATLAB 实例它用网格扫描和随机撒点两种写法算同一个阴影区域分别跑出 7.2150 和 7.2500。适合刚接触数值计算、要写建模论文、或者想弄明白随机算法收敛过程的人也适合手里有 MATLAB 想找个能立刻跑通的小项目练手的开发者。2. 两种撒点思路网格扫描与随机抽样的代码骨架2.1 先把阴影区域翻译成三个不等式代码里所有的判断都建立在这三个式子上先把它们讲清楚后面看循环才不费劲。阴影区域是三个几何体的交集一个是横放的椭圆半轴分别为 3 和 6标准方程写成 x²/9 y²/36 1一个是竖放的椭圆半轴分别为 6 和 1写成 x²/36 y² 1还有一个圆心在 (2,-1)、半径为 3 的圆写成 (x-2)² (y1)² 9。一个点要同时满足这三条才算落在阴影里。function tf inShadow(x, y) % 判断点 (x,y) 是否落在三区域交集内 a1 x^2/9 y^2/36; % 横椭圆半轴 3 与 6 a2 x^2/36 y^2; % 竖椭圆半轴 6 与 1 a3 (x-2)^2 (y1)^2; % 圆圆心 (2,-1)半径 3 tf (a1 1) (a2 1) (a3 9); end这里a1、a2存的是归一化后的左半部分所以阈值是 1a3存的是距离平方阈值是半径的平方 9。写成「距离平方小于 r²」而不是sqrt(...) 3是为了省掉一次开方在百万级采样下这点开销并不小。三个条件用串联MATLAB 会短路求值第一个条件不成立就不再算后面两个随机撒点时大约一半的点在第一个椭圆外这个短路实际能省掉三成左右的判断量。外层参考区域取 [-6,6]×[-6,6]边长 12面积 144这个 144 在两种方法里都会作为乘数出现。2.2 方法一等步长网格遍历方法一的思路是「均匀撒豆子」。把 x 轴从 -6 到 6 按 0.01 的步长切开得到 1201 个刻度y 轴用同一组刻度于是平面上铺开一张 1201×1201 的格点网总计 1442401 个点。逐点判断是否在阴影里命中的点数除以总点数再乘 144就是面积估计。原始代码是这样的clear x -6:0.01:6; % 1201 个刻度 y x; % y 用同一组刻度 s size(x); zs s(1,2)^2; % 网格总点数 1201*1201 1442401 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 if a2 1 if a3 9 k k 1; % 三个条件都满足才计一次 end end end end end mj (12^2) * k / zs; % 面积 144 * 命中比例s(1,2)取的是行向量的列数zs是总格点数k是命中数最后一行把比例乘回参考面积。这段代码跑出来mj 7.2150。三重嵌套if其实可以合并成一行if a11 a21 a39效果完全一样原作者那样写大概是为了让每个几何条件各自可读。clear放在开头是清理工作区在脚本里算习惯动作但如果你把它贴进函数体clear会把传进来的参数也清掉这是个经典翻车点后面避坑章还会提。这段代码的执行代价值得算一笔账1442401 次循环每次做三次乘方、四则运算和一串比较在普通笔记本的 MATLAB 上大致是几秒到十几秒量级。慢的根源不在算法而在 MATLAB 对逐元素循环的解释执行开销——每个x(i)都是一次带边界检查的索引而这段代码里x(i)^2在内层循环被重复算了 1201 次。2.3 方法二随机撒点与百次取均值方法二走的是另一条路不铺网格直接在 [-6,6] 里随机甩点。每次生成一个二维随机坐标a 12*rand(1,2)-6rand给的是 (0,1) 上的均匀分布乘 12 再减 6就映射到 (-6,6)。撒 N10000 个点统计命中数 k单次估计值就是 144*k/N。单次估计会偏高或偏低所以外面再套一层循环把整个过程重复 100 次最后取均值。clear N 10000; % 每次撒点数量 n 100; % 重复次数 for j 1:n k 0; for i 1:N a 12*rand(1,2) - 6; % 映射到 (-6,6) 的随机坐标 x(i) a(1,1); y(i) a(1,2); a1 (x(i)^2)/9 (y(i)^2)/36; a2 (x(i)^2)/36 y(i)^2; a3 (x(i)-2)^2 (y(i)1)^2; if a1 1 if a2 1 if a3 9 k k 1; end end end end m(j) (12^2) * k / N; % 第 j 次的单次估计 end mj mean(m); % 100 次估计的均值m(j)保存每一轮的估计值mean(m)是最终答案运行结果mj 7.2500。跟方法一相比这段代码的x(i)、y(i)是在循环里动态增长的数组MATLAB 每加一个元素都要重新分配内存N10000 时还感觉不出来N 提到 1e6 就会明显拖慢——这也是后面向量化要解决的问题。另外注意clear之后再执行m的长度是 100mean(m)求的是这 100 个数的平均如果中途改过n而没clearm里会残留上一轮的值均值就被污染了。两种方法都用了同一套a1、a2、a3判定区别只在点的来源和结果的处理方式方法一是确定性的给定步长结果唯一方法二是随机的每跑一次结果都略有不同靠重复取均值压住波动。2.4 两种方法的执行链路对照把两条链路摊开放在一起看选型就清楚了。网格法的点数由步长决定步长减半点数翻四倍误差随步长线性下降但边界处那圈格点的判定是「非黑即白」会带来系统性偏差随机法的误差随 1/√N 下降想精度翻倍得把采样量提到四倍好处是不受区域形状限制再怪的边界也只需要一个内点判定函数。维度方法一网格扫描方法二随机撒点点的来源-6:0.01:6等距刻度12*rand-6均匀随机总样本量1201² 144240110000 × 100 1000000结果是否可复现是同一份代码永远同一结果否每次跑都有细微差别误差来源边界格点判定带来的系统偏差抽样波动可用统计量估计实测结果7.21507.2500适合场景区域规则、想要确定性输出区域复杂、高维、解析解难求两者相差 0.035占面积约 0.5%这个差距正好落在随机法单次采样量 10000 时的统计波动范围内——单次估计的标准差约 0.31100 次均值的标准误约 0.031与方法一的结果差一个标准误出头完全正常。理解这一点才不会看到两个结果不一致就以为代码写错了。3. 参数与精度N、步长和重复次数到底怎么定3.1 步长与点数的平方关系网格法的样本量是「每个方向的刻度数」的平方这个平方关系很容易被低估。步长 0.1 时每轴 121 个点总共 14641 个点跑起来瞬间完成步长 0.01 时每轴 1201 个点总共 1442401 个点翻了一百倍步长 0.001 就是 1.44 亿个点内存和时间都开始吃不消。换算成代码就是for step [0.2 0.1 0.05 0.02 0.01] x -6:step:6; n numel(x); fprintf(步长 %.3f - 每轴 %5d 点, 总计 %10d 点\n, step, n, n^2); end输出里能直观看到步长每缩小一半总点数涨四倍。精度提升却只是线性的边界附近那条宽度为一个步长的带子里格点的归属可能判错误差量级大致是「边界长度 × 步长 / 2」。阴影区域周长粗估 20 左右步长 0.01 时误差上界约 0.1实际的系统偏差因为正负抵消会更小这跟观测到的 7.2150 与真值之间的偏差量级吻合。想把网格法精度再提一个数量级得把点数提两个数量级代价不成比例。3.2 随机法的误差公式与置信区间随机法的好处是误差可以算出来不用靠猜。设命中概率为 p面积真值为 A参考区域面积为 S144单次撒 N 个点时命中数 k 服从二项分布估计值 Â S·k/N 的标准差为p k / N; % 命中比例 se 144 * sqrt(p*(1-p)/N); % 单次估计的标准差 ci area [-1.96, 1.96] * se; % 95% 置信区间代入 p ≈ 7.25/144 ≈ 0.0504、N 10000算出来 se ≈ 0.31也就是说单次跑一亿分之一个点的结果会在 7.25 上下晃 ±0.6 的范围里。这就是为什么原始代码必须重复 100 次100 次均值的标准误降到 0.031波动被压到小数点后两位稳定。如果只跑一次就下结论那 7.1 和 7.4 都算「正常结果」讨论精度就没有意义了。这条公式还告诉你一个反直觉的事实p 接近 0 或 1 时方差最小。阴影面积只占参考区域的 5%这其实是好事同样的 N 下相对误差更小。如果换成一个占满大半个方框的区域误差反而会变大。3.3 跑一条收敛曲线看结果稳不稳判断代码有没有问题最省事的办法是画收敛曲线固定随机种子之外的设置把 N 从 100 取到 100 万每个量级重复若干次看估计值和置信带怎么收窄。这段代码可以当模板用Ns [1e2 1e3 1e4 1e5 1e6]; R 20; % 每个 N 重复次数 res zeros(numel(Ns), R); for a 1:numel(Ns) for r 1:R x 12*rand(Ns(a),1) - 6; y 12*rand(Ns(a),1) - 6; in (x.^2/9 y.^2/36 1) (x.^2/36 y.^2 1) ... ((x-2).^2 (y1).^2 9); res(a,r) 144 * nnz(in) / Ns(a); end end errorbar(1:numel(Ns), mean(res,2), std(res,0,2), -o); set(gca, XTick, 1:numel(Ns), XTickLabel, {1e2,1e3,1e4,1e5,1e6});nnz统计逻辑数组里真值的个数比sum快一点也更好读。跑完能看到误差棒按 1/√N 收窄N 从 1e4 提到 1e6误差棒缩到十分之一。如果误差棒没按这个规律缩或者均值随 N 系统性漂移那多半是判定条件或者坐标系出了问题而不是随机性作祟。曲线还能帮你定预算。要求结果稳定在小数点后两位N 取 1e6 单次就够要求三位就得 1e8 或者多次平均。建模论文里报结果写「7.25 ± 0.0395% CI10⁶ 次采样」比光写一个 7.25 有说服力得多。4. 从能跑到好用向量化重写与加速4.1 用 meshgrid 干掉双重循环原始网格法的双重循环在 MATLAB 里属于「能跑但不够像 MATLAB」的写法。把 x 和 y 两组刻度用meshgrid展开成两个 1201×1201 的矩阵三个条件用逻辑运算符一次性算完整个网格的判断就变成三次矩阵运算加两次逻辑与step 0.01; x -6:step:6; [X, Y] meshgrid(x, x); % 两个 1201x1201 矩阵 inside (X.^2/9 Y.^2/36 1) ... % 横椭圆 (X.^2/36 Y.^2 1) ... % 竖椭圆 ((X-2).^2 (Y1).^2 9); % 圆 k nnz(inside); area 144 * k / numel(X);numel(X)就是 1442401和原代码的zs一致所以结果仍然是 7.2150但耗时从秒级降到几十毫秒。注意这里用的是而不是只对标量短路矩阵运算必须用逐元素的这是从循环改写成向量化时最常见的翻车点写完报错「Operands to the || and operators must be convertible to logical scalar values」就是它。内存上要留个心1201×1201 的 double 矩阵一个占 11.5 MBmeshgrid出两个、中间结果再来几个峰值能到一百多 MB。步长降到 0.002 时每轴 6001 个点单个矩阵 288 MB普通机器就开始吃紧了。这种时候要么退回循环加分块要么换随机法。4.2 随机版本的分块与 parfor随机法向量化更彻底一次生成 N 个点判定一次算完N 1e6; x 12*rand(N,1) - 6; y 12*rand(N,1) - 6; in (x.^2/9 y.^2/36 1) (x.^2/36 y.^2 1) ... ((x-2).^2 (y1).^2 9); area 144 * nnz(in) / N; se 144 * sqrt(mean(in)*(1-mean(in))/N);N1e6 时两个坐标向量各 8 MB逻辑数组 1 MB内存毫无压力耗时在零点几秒。要把 N 提到 1e8一次性生成就会顶到内存上限这时候按块处理function [area, se] mcAreaBlock(N, block) % 分块蒙特卡罗block 为每块点数 k 0; for s 1:block:N m min(block, N - s 1); % 最后一块可能不满 x 12*rand(m,1) - 6; y 12*rand(m,1) - 6; in (x.^2/9 y.^2/36 1) (x.^2/36 y.^2 1) ... ((x-2).^2 (y1).^2 9); k k nnz(in); end p k / N; area 144 * p; se 144 * sqrt(p*(1-p)/N); endmin(block, N-s1)处理最后一块不满的情况漏掉这一句会多算几个点虽然对结果影响微乎其微但统计口径就不干净了。分块还有个附带好处可以在循环里打印进度不然 1e8 次采样跑起来屏幕一片空白你都不知道它是不是卡死了。如果机器是多核把重复实验那层循环换成parfor最划算——每轮实验互相独立天然适合并行R 100; areas zeros(1, R); parfor r 1:R x 12*rand(1e5,1) - 6; y 12*rand(1e5,1) - 6; in (x.^2/9 y.^2/36 1) (x.^2/36 y.^2 1) ... ((x-2).^2 (y1).^2 9); areas(r) 144 * nnz(in) / 1e5; end fprintf(均值 %.4f标准误 %.4f\n, mean(areas), std(areas)/sqrt(R));parfor的限制是循环体里不能有跨轮次的数据依赖这里每轮只写自己的areas(r)符合要求。另外rand在并行池里每个 worker 的随机流是独立的不会出现多轮结果完全一样的情况。4.3 把区域判定抽成函数句柄代码写到这一步最该做的改造是把「算什么区域」和「怎么采样」拆开。区域判定塞进一个匿名函数采样器只负责撒点和统计换区域时只改一行% 区域判定输入同尺寸的 x、y返回逻辑数组 insideFcn (x,y) (x.^2/9 y.^2/36 1) ... (x.^2/36 y.^2 1) ... ((x-2).^2 (y1).^2 9); [area, se] mcIntegrate(insideFcn, [-6 6], [-6 6], 1e6);对应的通用采样器function [area, se] mcIntegrate(fcn, xlim, ylim, N) % fcn : 接受 (x,y) 向量、返回逻辑数组的判定函数 % xlim : [xmin xmax]采样矩形范围 % ylim : [ymin ymax] % N : 采样点数 x xlim(1) diff(xlim) * rand(N,1); y ylim(1) diff(ylim) * rand(N,1); in fcn(x, y); S diff(xlim) * diff(ylim); % 参考矩形面积 p mean(in); area S * p; se S * sqrt(p*(1-p)/N); enddiff(xlim)拿区间长度S是参考矩形面积这样函数不再绑死 144 这个常数。判定函数必须是「向量进、逻辑数组出」因为采样器一次传进去的是整列坐标如果写成标量函数就得在外面套arrayfun性能会掉一个档次。这个接口设计好之后算圆、算心形线、算高维超球体积都只是换一个insideFcn的事。5. 避坑与排查五个让人白跑的细节5.1 坐标与判定条件里的坑现象代码跑出来面积接近 0或者明显不到 7 的一半。原因随机坐标忘了映射回目标区间。12*rand(1,2)出来的是 (0,12)如果不减 6所有点都落在第一象限而阴影区域主要分布在原点附近偏右侧命中率被严重扭曲。解决写完坐标生成那行先disp([min(x) max(x)])看一眼实际范围是不是落在 [-6,6]。这个检查比调试判定条件便宜得多。现象结果比 7.2 大一大截或者干脆等于 144 乘以某个离谱的比例。原因三个条件之间用了||或者漏了一个条件变成求并集而不是交集。原作者用三重嵌套if保证「与」的关系改成单行时很容易顺手写成||。解决把三个条件分别单独统计命中率看三个数是不是依次递减——交集一定比其中任何一个都小。三个单条件命中率大概在 5% 到 20% 之间交集的 5% 应该是最小的那个。现象换了个图形明明参数改对了结果却总差一点。原因椭圆半轴和分母搞混。x²/9 y²/36 1对应 x 方向半轴 3、y 方向半轴 6写成x²/36 y²/9就变成横竖互换面积完全不同。解决拿一个已知点验证比如原点 (0,0) 一定在区域内点 (5.5,0) 一定在区域外两个椭圆的 x 半轴都不到 6。先跑这两个点判定函数对了再撒点。5.2 精度与性能上的坑现象单次撒 1000 个点结果在 6.8 到 7.7 之间跳怀疑代码有 bug。原因这不是 bug是抽样波动。N1000 时单次标准误约 1.0跳这么大幅度完全正常。解决要么把 N 提到 1e5 以上要么固定跑 100 轮取均值并报出标准误。只看单次结果就判断对错是这类代码最常见的误判来源。现象把 N 提到 1e7MATLAB 卡住不动风扇狂转。原因循环里x(i) ...动态扩容每加一个元素重分配一次内存复杂度是 O(N²) 级别的拷贝再加上脚本开头放了clear如果这段代码被贴进函数还会把输入参数清掉。解决进循环前x zeros(1,N); y zeros(1,N);预分配或者干脆用 4.1、4.2 节的向量化写法。clear只留在真正独立的脚本最外层函数里一律不用。现象网格法步长从 0.01 改到 0.005结果只在小数点后第三位变了感觉白算。原因网格法的误差是 O(步长) 的步长减半误差只减半但点数翻四倍。这个投入产出比本来就不划算。解决要更高精度优先用随机法加采样量或者两种方法交叉验证——两者结果落在彼此的置信区间内就说明实现没问题没必要死磕网格步长。6. 二次开发把这十几行改成能复用的面积计算器这份代码最值得带走的东西不是 7.2150 这个数而是「参考矩形 内点判定 比例乘面积」这个三段式结构。把它固定下来剩下的都是可替换件。我通常会把 4.3 节的mcIntegrate作为内核外面再包一层参数扫描用来回答「区域参数变化时面积怎么变」这类问题——比如把圆心从 (2,-1) 往右挪看交集面积怎么缩offsets -1:0.5:3; areas zeros(size(offsets)); for i 1:numel(offsets) cx offsets(i); f (x,y) (x.^2/9 y.^2/36 1) (x.^2/36 y.^2 1) ... ((x-cx).^2 (y1).^2 9); areas(i) mcIntegrate(f, [-6 6], [-6 6], 5e5); end plot(offsets, areas, -o);每换一个cx都要重跑 5e5 次采样五次下来不到两秒比解析推导快得多而且改区域形状时不用重新推公式。做参数敏感性分析、给建模论文配图这套流程能直接搬。再往上一层可以把内点判定换成任意谓词比如「到三个定点的距离之和小于某值」这种没有解析面积的定义蒙特卡罗照样能给出带置信区间的答案这是它相对解析法真正的优势。验证方法也简单拿一个已知答案的图形半径 3 的圆面积 9π ≈ 28.27跑一遍通用函数结果落在 28.27±0.1 内就说明采样器本身没问题再去换复杂区域。最后一个技巧是关于结果呈现的任何时候都别只输出一个数。把area、se、采样量N一起打出来写成fprintf(面积 %.4f ± %.4f (N%d)\n, area, 1.96*se, N)。我早年交建模论文时只写了个 7.25被问「你这个数误差多大」当场答不上来从那以后每次跑随机算法都强制把置信区间一起打出来哪怕没人问。希望这个习惯也能帮到你。本文还有配套的精品资源点击获取