企业官网建设流程全解析

简介：围绕蒙特卡罗方法在 MATLAB 中的面积计算实例，这份 doc 文档面向正在学习数值计算、随机模拟与 MATLAB 编程的读者，重点解决如何用随机抽样思路近似求解不规则阴影区域面积的问题。压缩包内仅含 1 个 doc 文件，大小约 33KB，内容紧凑，适合作为课堂示例、建模入门或自学参考。文档给出均匀网格分割与随机投点两套实现方案，附有可运行程序、关键判断条件、循环与向量化写法，以及 7.2150 与 7.2500 两组运行结果，便于对照理解精度来源、样本量设置和多次试验取均值的处理方式。已有 1626 人学习，说明该案例在蒙特卡罗入门与 MATLAB 实践中具有一定参考价值。读者可从中获得完整的程序框架、两种方法的差异分析、随机数生成与结果验证思路，并据此迁移到面积估算、物理现象模拟等类似问题。

1. 7.25 这个数是怎么数出来的

7.25，是把一万个随机点撒进 12×12 的正方形、数出落在阴影里的那些点之后乘回来的面积。蒙特卡罗在 MATLAB 里的入门例子，十个里有八个是算圆周率，这份实例走的是另一条路——算一个由两个椭圆加一个圆叠出来的阴影区域。它值得单独拆一遍，原因很实在：这块区域的解析解要分段积分，写起来不短，而数点法三行代码就能出数，误差还能自己算出来。适合谁：手头有 MATLAB、需要快速估一个不规则区域面积，或者拿它验证某个积分近似值的人；也适合刚碰随机算法、想找一个能跑通、能改参数、能亲眼看到误差怎么变的练手项目。下面的内容不重复简介里的文件清单，直接从三个不等式怎么写、两种采样怎么选、哪里容易翻车讲起。

2. 三个不等式定生死：把阴影区域翻译成可判定的逻辑

拿到一个“求阴影面积”的图，第一件事不是打开 MATLAB，而是拿笔把边界写成不等式。写错了，后面采样再密也是白搭；写对了，剩下的只是数点。

2.1 三条不等式：从图形到布尔条件

这份实例的阴影区是三个区域同时成立的部分：一个竖放的椭圆、一个横放的椭圆、一个斜着偏在第四象限的圆。原始脚本里用a1、a2、a3三个变量表示，写法是标准形式：

判断一个点是否落在阴影里，就是把这三条同时取真。用集合语言说是交集，用代码说是逻辑与。这里有个细节值得停一下：a2里y^2的系数是 1，说明 y 半轴等于 1，不是 6；这个数写错一位，面积能差出好几倍，而且结果看上去还挺像那么回事，不会报错。我见过有人把a2抄成x^2/36 + y^2/6 < 1，跑出来 12 点多，还以为是自己 N 不够大。

还有一个隐含前提容易被忽略：这三个区域必须在参考正方形内有界。竖椭圆 y 方向到 ±6，横椭圆 x 方向到 ±6，正好顶到 [-6,6]² 的边；圆心 (2,-1) 半径 3，最远到 x=5、y=-4，也在框内。也就是说参考区域完全包住了阴影区，没有任何一部分漏在框外。如果哪天换成半径 5 的圆，最远到 x=7，超出边界，那 144 这个参考面积就不成立了，数出来的比例会系统性偏大。

2.2 判定写成一个函数：向量化以后省掉三层 if

原始脚本用的是三层嵌套if，每判断一个点要走进三层分支。逻辑没错，但把判定和循环缠在一起，改一个条件要动三处，也没法复用。我一般会把它抽成一个独立函数，输入输出都允许是数组：

function tf = inShadow(x, y) % 判断点 (x,y) 是否落在阴影区（三个区域取交集） % x、y 可以是标量、向量或同尺寸矩阵，返回同尺寸逻辑数组 tf = (x.^2/9 + y.^2/36 < 1) & ... % 竖椭圆 (x.^2/36 + y.^2 < 1) & ... % 横椭圆 ((x-2).^2 + (y+1).^2 < 9); % 圆心 (2,-1)、半径 3 end

逻辑说明部分：点运算符.^和./让表达式对数组逐元素生效，一个矩阵进来，一个逻辑矩阵出去，中间没有任何循环。三个条件之间用&而不是&&，因为&&只对标量做短路求值，喂矩阵进去会直接报错——这是新手改这段代码时最常见的报错来源。另外 MATLAB 里&对两个逻辑数组是逐元素与，不短路的代价是三个条件都会被完整算一遍，但换来的是整块矩阵运算，比逐点判断快得多。

参数说明：<用的是严格小于。对连续均匀分布来说，边界上取不取等号不影响期望值，因为落在曲线上的概率是零。但对规则网格来说就有影响了，这一点放到第 4 章细说。

写完先别急着跑大循环，拿几个已知点验一下：

>> inShadow(0, 0) ans = logical 1 >> inShadow(5, 0) ans = logical 0 >> inShadow(2, -1) ans = logical 1

(0,0) 三个条件分别给出 0<1、0<1、5<9，全真；(5,0) 第一条 25/9≈2.78>1，直接出局；(2,-1) 是圆心，必然在内。这三行能挡掉大部分抄错系数的低级错误。

2.3 参考面积 144 什么时候会骗你

原始脚本里mj=(12^2)*k/zs，12^2就是参考正方形面积，(12)² = 144。这个数能直接写死，前提是采样范围严格等于 [-6,6]×[-6,6]。方法一里x=-6:0.01:6，步长 0.01 能整除 12，点数 1201，左右端点都取到，覆盖面积正好是 144；zs=s(1,2)^2用点数的平方当分母，等于假设 x 和 y 方向的网格完全一致。

这个假设在步长不整除的时候会崩。比如把步长改成 0.013，-6:0.013:6的最后一个点落在 5.994，网格实际覆盖的是 [-6, 5.994]，面积 11.994² = 143.86，不是 144。误差 0.1%，在精度要求不高的时候看不出来，一旦拿来做收敛性分析，这条系统性偏差会一直挂在那里，怎么加密度都消不掉。稳妥的写法是让分母跟着实际覆盖范围走：

g = -6:h:6; [X, Y] = meshgrid(g, g); k = nnz(inShadow(X, Y)); covered = (g(end)-g(1)) * (g(end)-g(1)); % 实际覆盖面积，不写死 144 A = covered * k / numel(X);

nnz数非零元素，逻辑数组里就是数 true 的个数，比自己写sum(sum(M))干净。numel(X)是全网格点数，比s(1,2)^2更稳——它不关心你是方形网格还是矩形网格，哪怕以后把 x、y 步长设成不一样，这条式子照样成立。这几行改动看着小，是把一段“只能跑这一个例子”的脚本改成“换个区域也能跑”的第一步。

3. 网格法与随机法：两种撒豆子的方式，账本不一样

原始脚本给了两条路，跑出来一个 7.2150，一个 7.2500。差 0.035，约 0.5%。这个差距不是谁写错了，是两种采样策略各自的误差结构不同，理解这一点比记住结果重要得多。

3.1 网格法：步长 0.01 到底意味着多少个点

方法一的思路是把正方形切成规则格子，逐格判断。x=-6:0.01:6生成 1201 个坐标，meshgrid铺开就是 1201×1201 = 1,442,401 个点。原始脚本用两层for逐个走这 144 万个点，在现在的机器上大概能跑完，但纯属浪费——每一轮循环里做的都是同样的三个表达式，完全可以一次性算完。

h = 0.01; % 网格步长，直接决定点数：n = 12/h + 1 g = -6:h:6; [X, Y] = meshgrid(g, g); % X、Y 都是 1201×1201 M = inShadow(X, Y); % 一次算完 144 万个点的判定 k = nnz(M); A_grid = 144 * k / numel(M); % 结果约 7.2150

逻辑说明：meshgrid把两个一维坐标铺成二维网格，X 的每一行都是g，Y 的每一列都是g，这样inShadow里的每个表达式都作用在整张网格上。nnz(M)数出命中的格点数，除以总格点数得到比例，再乘参考面积。整个过程只有三次矩阵运算，没有显式循环。

参数说明：h是这段代码唯一需要调的旋钮。它和点数的关系是点数 = (12/h + 1)²，减半步长等于点数翻四倍。h=0.01 出 144 万个点，h=0.005 就是 577 万个点，内存占用跟着涨——两个 1201×1201 的 double 矩阵大概 23 MB，到 h=0.0025 就要 370 MB 上下，普通笔记本还能扛，再细就得考虑分块了。

网格法的误差来源很明确：边界处的格子要么整个算进去，要么整个不算，实际边界穿过的那些格子被粗暴地二值化了。步长越小，被误判的格子占总数的比例越低，误差随 h 线性下降。这个特性决定了它适合低维、边界相对规整、需要确定性结果的场合——同一份代码跑一百遍，结果一位不差。

3.2 随机法：为什么单次结果不能信，要取 100 次均值

方法二才是教科书意义上的蒙特卡罗。早年间在一次建模培训上，某导师讲过一个故事：古人要量一块不规则图形的面积，把图样画在方布上，抓一把豆子随手撒上去，数落在图形内的豆子重量，比上全部豆子的重量，再乘布的面积，就是要求的面积。这个故事我记了很多年，因为它把“随机撒点 + 数比例”讲得比任何公式都清楚。

代码落地是这样的：

rng(20240501); % 固定随机种子，保证可复现 N = 10000; % 单次撒点数 nRep = 100; % 重复次数 est = zeros(nRep, 1); for r = 1:nRep X = 12*rand(N,1) - 6; % 均匀撒在 [-6,6) Y = 12*rand(N,1) - 6; est(r) = 144 * nnz(inShadow(X, Y)) / N; end A_mc = mean(est); % 约 7.2500 A_se = std(est) / sqrt(nRep); % 均值的标准误

逻辑说明：rand(N,1)产生 [0,1) 上的均匀随机数，乘 12 减 6 就映射到 [-6,6)。注意右端点取不到，但连续分布下这不影响结果。每次撒 N 个点，算出一个面积估计存进est，重复 100 次后取均值。取均值这一步不是可有可无的装饰——单次估计的涨落能到 ±0.3 左右，直接报单次结果等于报了一个随机数。

参数说明：N控制单次精度，nRep控制均值的稳定性。这两者的作用不能互相替代：把 N 从 1 万加到 100 万，单次标准差降到原来的十分之一；而把重复次数从 100 加到 10000，均值的标准误也降到十分之一，但总采样点数是一样的。区别在于，前者同时改善了每一次估计，后者只改善最终均值。如果只需要一个数，两者等价；如果要看分布，重复次数不能太少。

rng(20240501)这行是血泪经验。原始脚本没写种子，rand会用系统时间初始化，同一份代码今天跑 7.2500，明天跑 7.1832，把结果贴给别人，对方复现不出来，沟通成本陡增。种子值本身随便填，写进脚本注释里就行。

3.3 两种采样方式的账本：精度、耗时、可复现性

把两条路的特性摆在一张表里，选型就不用拍脑袋了。

这张表里最关键的是倒数第二行。三维空间里步长减半，点数涨 8 倍；十维空间里涨 1024 倍，直接跑不动。蒙特卡罗的真正价值不在二维面积——二维用网格法又快又准——而在于维度一高，它是唯一还能线性扩展的办法。这份实例用二维面积做载体，练的是“把问题写成指示函数 + 采样 + 求均值”这套动作，这套动作搬到十维积分上一样成立。

还有一点：两者的结果不能混着比。7.2150 和 7.2500 差 0.035，看着不大，但网格法那边是确定性偏差，随机法那边是统计涨落，性质不同。真要判断谁更准，得先估出随机法的标准误——按 p≈7.25/144≈0.0503 算，单次 N=1e4 的标准差约 0.315，100 次均值的标准误约 0.0315。也就是说 7.2500 这个数本身的波动范围就有 ±0.06（95% 置信），和 7.2150 差了不到一个多标准误，两者是一致的，不存在谁对谁错。

4. 避坑与排查：五条踩过的坑，按现象对号入座

这段代码短，坑却不少。下面五条是我自己改这份脚本时真踩过的，按“现象→原因→解决”写，遇到问题直接对号入座。

4.1 同一份代码两次结果不一样

现象：把脚本发给同事，他跑出 7.1832，你这边是 7.2500，代码一字不差，怀疑是不是 MATLAB 版本问题。

原因：rand默认用系统时间做种子，每次启动都是新的随机序列。没固定种子，脚本就不是确定性的程序，而是一个随机实验。

解决：在脚本最前面加rng(某个固定整数)。种子值随手填一个六位数即可，重要的是写进注释并随代码一起交付。如果要看结果的分布而不是单点，就别固定种子，改成固定nRep并把est整个数组保存下来，附上均值和标准误一起报。

4.2 把 N 加大十倍，精度没跟着涨十倍

现象：N 从 10000 加到 100000，以为误差能降一个数量级，实际只降了约 3 倍，感觉代码有 bug。

原因：蒙特卡罗的误差按 1/√N 收敛，这是方法本身的数学性质，不是实现问题。N 涨 10 倍，标准差只降到 1/√10 ≈ 0.316 倍。

解决：认清这条曲线，别跟它较劲。要把误差降一个数量级，N 得涨 100 倍。如果算力吃不消，就换降方差的手段（分层抽样、控制变量），或者干脆承认这个精度够用——工程上估个面积，1% 的相对误差通常远远够。

4.3 网格法结果总是比随机法小一点

现象：网格法跑出 7.2150，随机法 7.2500，反复调步长，网格法的结果总在 7.21 附近晃，上不去。

原因：两个因素叠加。一是严格小于号把边界上的点排除了；二是规则网格的采样点位置固定，边界穿过的格子被整格判为内或外，偏差不会正负抵消，而是单向累积。

解决：把采样点从格点改成格心，g = -6+h/2 : h : 6-h/2，每个格子用中心点代表，边界附近的判定更接近真实占比。想验证是不是这个原因，把 h 从 0.01 减到 0.005 再跑一次，如果结果单调往 7.25 靠，就确认是量化偏差而非逻辑错误。

4.4 逐点循环在大 N 下跑不动

现象：把原始脚本的嵌套for直接用在 N=1e6 上，进度条半天不动，风扇狂转。

原因：MATLAB 的循环虽然这些年优化过，但每个点都要走一遍解析器和三层if分支，开销远大于矩阵运算。144 万个点还好，1 亿个点就是灾难。

解决：把判定写成接受数组的函数（见 2.2），采样点一次性生成成列向量，整块传进去。向量化以后，同样的 N 只受内存带宽限制，不受解释器开销限制。注意监控内存：一次生成 1e7 个 double 列向量是 80 MB，同时有好几个临时变量就得几百 MB，机器小的分批处理。

4.5 分母写死 144 导致结果对不上

现象：换了采样范围或者改了步长，结果和手算的比例对不上，差个千分之几。

原因：mj=(12^2)*k/zs里 144 和zs都是写死的，隐含“采样范围正好是 [-6,6]²、网格正好是正方形”两个假设。任一假设被打破，分母就不对了。

解决：参考面积用实际覆盖范围算，点数用numel或实际长度算，见 2.3 的改法。这条在只跑原始例子时不会暴露，一旦想把这套代码用到别的图上，就是第一个绊脚石。

5. 从算一个面积到算一类积分：把脚本改成能复用的函数

原始脚本是为这一张图写的。真要反复用，得把“撒点、判定、求均值、给误差”这套动作封起来，让换一个区域只需要换一个判定函数。

5.1 抽象成估计期望：面积只是被积函数取 0/1 的特例

蒙特卡罗求面积的数学形式是 A = V · E[f(U)]，其中 V 是参考区域体积（这里是 144），U 在参考区域上均匀分布，f 是示性函数——点落在阴影里取 1，否则取 0。E[f(U)] 就是命中概率 p，用样本均值 k/N 去估。

这个形式一写出来，适用范围立刻宽了：把 f 从 0/1 换成任意函数，同一套代码就能算 ∫f dV。算物理里的碰撞概率、金融里的期望收益、可靠性里的失效概率，骨架完全一样。区别只在于有没有解析解、方差多大、要不要降方差。

我一般建议初学时先用示性函数跑通，因为结果能和几何面积对照，错了看得出来。等这套动作熟了，再换成连续被积函数，那时候没有几何直觉兜底，只能靠收敛曲线和标准误判断对错。

5.2 封装成 mcArea：接口设计与调用示例

把重复采样的骨架抽出来，判定逻辑作为函数句柄传进去：

function [A, se, raw] = mcArea(indicator, N, nRep, seed) % indicator: 函数句柄，接收 X、Y 两个列向量，返回同长度逻辑向量 % N: 单次采样点数；nRep: 重复次数；seed: 随机种子 % 返回：面积估计 A、均值的标准误 se、每次的原始估计 raw if nargin < 4 || isempty(seed) seed = 1; end rng(seed); raw = zeros(nRep, 1); for r = 1:nRep X = 12*rand(N,1) - 6; Y = 12*rand(N,1) - 6; raw(r) = 144 * mean(indicator(X, Y)); end A = mean(raw); se = std(raw) / sqrt(nRep); end

逻辑说明：把参考正方形、采样、求均值这三件事固定下来，变化的部分只有indicator。mean作用在逻辑向量上会自动转成 0/1 求平均，等价于nnz/N，少一次转换。返回值多给一个se，是为了让调用方知道这个数有多可信，而不是把它当成精确值。

参数说明：N默认给 1e4，nRep默认给 100，这两个是原始脚本验证过的组合。seed留空时用 1，保证默认行为可复现。调用时这样写：

[A, se] = mcArea(@inShadow, 1e4, 100, 20240501); fprintf('面积 %.4f，标准误 %.4f\n', A, se);

换成别的区域，只写一个新的indicator函数，其余不动。要是参考区域不是正方形，把里面写死的12、-6、144提成参数就行，结构不用改。

5.3 降方差三板斧：分层抽样有用，对偶变量对 0/1 无效

误差按 1/√N 收敛这件事改变不了，但可以改前面的系数。

分层抽样是把参考区域切成 m×m 个子块，每个子块内独立撒点。子块内命中率的波动比整块小，合起来整体方差就降了。这是三板斧里对示性函数最实在的一招：

function A = mcStratified(indicator, m, perCell, seed) % m: 每边分层数；perCell: 每个子块内撒点数 rng(seed); h = 12/m; xs = -6 + h*(0:m-1)'; % 每个子块的左边界 acc = 0; for i = 1:m for j = 1:m X = xs(i) + h*rand(perCell,1); Y = xs(j) + h*rand(perCell,1); acc = acc + nnz(indicator(X, Y)); end end A = 144 * acc / (m*m*perCell); end

参数说明：m每边分几层，perCell每块撒几点，总采样数 = m²·perCell。同样撒 100 万个点，分层版把 100 万撒在一个大方块里，还是撒在 400 个小方块里各 2500 个，方差的差别在边界附近最明显——子块让边界穿过的区域被更均匀地覆盖了。m取 10 到 50 之间通常够，太大每块点数太少，单块内的估计反而不稳。

对偶变量（用 U 和 1-U 配对）对光滑被积函数很有效，但对 0/1 示性函数基本没用——f(U)和f(1-U)都是 0 或 1，要么同时命中要么不同时命中，协方差捞不到多少好处。这条我试过，白折腾半天，后来想明白是因为被积函数不连续，降方差技巧的收益和函数光滑性直接挂钩。控制变量同理，得有和 f 高度相关、期望已知的辅助函数才行，示性函数场景下不好构造。

6. 先算 N 再动手：误差换算与收敛曲线验证

写完代码别急着把 N 拉到最大，先估一下要多大的 N 才够。

6.1 一张换算表：N 和误差的对应关系

按 p ≈ 7.25/144 ≈ 0.0503、单次估计标准误 = 144·√(p(1-p)/N) 算：

这张表反过来用最有价值：要 1% 的相对误差，单次 N 得 1e5 左右；要 0.1%，得上 1e7。每次看到有人开口就要 0.01% 的精度，我都先请他算一下 N 需要多大——1e9 个点，光是生成随机数就够喝一壶，这时候该考虑的是降方差或者换个方法，不是加机器。

表里的数字是按这个特定区域的命中率算的，换个区域 p 变了，数值会变，但量级规律不变。我习惯在脚本里把se打印出来，跟结果一起报，看的人心里有数。

6.2 收敛曲线：一条比单个数更可信的证据

单次结果对不上预期时，最省事的诊断是画收敛曲线：把 N 取一组等比的值，每个 N 重复若干次取均值，画 log-log 图，看它是不是沿着斜率 -1/2 的直线往下走。

Ns = [1e3 3e3 1e4 3e4 1e5 3e5 1e6]; err = zeros(size(Ns)); for t = 1:numel(Ns) [A, se] = mcArea(@inShadow, Ns(t), 30, 100+t); err(t) = se; end loglog(Ns, err, 'o-'); hold on; loglog(Ns, err(end)*sqrt(Ns(end)./Ns), 'k--'); % 参考斜率 -1/2

如果实测曲线和虚线基本平行，说明采样逻辑没问题，误差就是纯粹的统计涨落；如果实测曲线明显高于虚线，或者在某一段突然抬起来，那多半是判定函数写错了、采样范围没覆盖全、或者随机数相关性出了问题——这三种情况都不是加 N 能救的。

从那以后我每次改这类采样脚本，都会先跑一遍这条曲线，再谈结果准不准。7.25 这个数字本身没什么稀奇，稀奇的是它能被一条斜率 -1/2 的线解释清楚。希望帮到你。

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