简介蒙特卡罗方法的MATLAB实现实例面向需要掌握随机抽样数值计算技巧的MATLAB初学者与数学建模爱好者。文档以计算特定区域阴影面积为例演示蒙特卡罗方法在几何数值计算中的具体应用给出两种典型实现方式一种是均匀分割坐标轴生成固定网格并统计落入目标区域的点数另一种是在指定范围内随机生成大量样本点并通过多次独立试验取均值来降低随机误差。两种方法均配有完整MATLAB代码、运行结果与对比分析运行数据分别得到7.2150和7.2500贴近真实面积可帮助读者直观理解算法原理并动手复现。资源为单个doc文档大小仅33KB轻量易获取。目前已有1626人学习适合作为蒙特卡罗入门练习或数值计算课程辅助材料。1. 蒙特卡罗算阴影面积先搞清楚它到底在算什么蒙特卡罗方法听上去很高深其实原理朴素得像撒豆子。我第一次在 MATLAB 交流群里看到有人问“这个阴影面积怎么求”时脑子里立刻蹦出当年数学建模培训时老师讲的故事古人把图形画在方布上撒一袋豆子落在图形里的豆子重量比上总重量再乘以布的面积就是图形面积。这个思路就是蒙特卡罗的核心——用大量随机样本来逼近解析解。这个应用案例很适合新手因为它只涉及随机数生成、条件判断和面积折算却把蒙特卡罗的两种典型玩法均匀网格与随机抽样都覆盖了。如果你刚接触蒙特卡罗或者想把随机模拟用到几何计算上这套 MATLAB 实现值得你亲手跑一遍。2. 面积问题的几何建模三个不等式与边界判定要算阴影面积第一步不是写代码而是把图形翻译成数学表达式。原始问题里阴影区域其实是由两个椭圆和一个圆取交集得到的三个嵌套的if条件分别对应三条不等式。理解这三条式子后面两种方法的代码逻辑就一目了然。2.1 阴影区域的数学描述观察代码里的判定条件阴影区域里的任意一个点(x, y)必须同时满足a1 (x^2)/9 (y^2)/36 1这是一个长轴在 y 轴上的椭圆半长轴 6半短轴 3。a2 (x^2)/36 y^2 1这是长轴在 x 轴上的椭圆半长轴 6半短轴 1。a3 (x-2)^2 (y1)^2 9这是一个圆心在(2, -1)、半径为 3 的圆。三个条件用嵌套if连接意味着点必须同时落在三个图形内部取的是交集区域。参考正方形选的是[-6, 6] × [-6, 6]边长为 12这是因为两个椭圆的长半轴都是 6圆也完全被这个正方形包住。这样一来整个采样空间就是一个12 × 12 144的方形区域阴影部分面积就是“满足条件的点数 / 总点数 × 144”。2.2 判定条件的嵌套逻辑代码里写的是if a1 1 if a2 1 if a3 9 k k 1; end end end这种写法等价于k k 1只有在三个条件全部成立时才会执行。三个if嵌套会造成不必要的缩进而且一旦条件增多代码会变得很难看。常见做法是用逻辑与合并if a1 1 a2 1 a3 9 k k 1; end是短路运算符如果第一个条件不满足后面的条件根本不会计算性能上会比挨个嵌套好一点。对于这个例子三条不等式都是多项式求值速度差异不大但代码可读性提升很多。注意这里用的是严格小于也就是说落在边界上的点不算数。在连续空间里边界点概率为零但网格法里可能出现点正好落在边界附近的情况后面避坑章我会专门说这个问题。2.3 参考正方形的选择依据选12 × 12不是拍脑袋决定的它直接决定了面积折算系数。如果参考区域选大了比如选[-10, 10]那么面积公式变成k/N × 400同样的k会得到更大的面积估计但阴影区域外的点更多浪费计算量选小了比如[-5, 5]会有一部分阴影区域跑到正方形外导致严重错误。这个例子中两个椭圆的长半轴都是 6圆也在 6 以内所以[-6, 6]是最紧凑的选择。实际项目里我会先画出图形的包围盒再留一点余量保证参考区域完全覆盖目标区域同时不浪费太多采样点。3. 方法一均匀网格分割的实现与精度瓶颈方法一的朴素想法是把整个正方形像棋盘一样均匀分割然后数有多少个格点落在阴影里。它其实不是严格意义上的随机蒙特卡罗而是确定性的网格积分法但因为思路一脉相承原题里也把它归为蒙特卡罗的一种变形。我建议你先跑这个版本它能帮你直观看到“点数越多面积越准”的趋势。3.1 网格法 MATLAB 代码原帖给的程序如下我保留了核心逻辑只加了注释clear x -6:0.01:6; % 设置 x 轴网格步长 0.01 y x; % 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 if a2 1 if a3 9 k k 1; end end end end end mj (12^2) * k / zs; % 面积 144 * 格点数 / 总格点数运行结果是mj 7.2150。这里的x -6:0.01:6生成了从 -6 到 6、间距 0.01 的向量每个方向上有1201个点总格点数是1201^2 1442401。每一点都通过三重循环去判断是否落在阴影区域最后把命中比例乘以参考面积 144就得到阴影面积的近似值。3.2 步长对结果的影响步长0.01决定了网格的密度。步长越小网格越密结果越接近真实面积但计算量按平方增长。0.01的步长在每边上有 1201 个点总循环次数约 144 万次在 MATLAB 里跑起来要一两秒还能接受。如果把步长改成0.001每边就有 12001 个点总循环约 1.44 亿次普通电脑可能要跑几分钟而且双循环的代码效率很低。我的建议是先用0.05粗跑一遍确认逻辑正确再逐步缩小步长观察结果变化。比如步长0.05时每边 241 个点结果大概在 7.2 左右但误差会偏大步长0.01时结果 7.2150这个值已经比真实面积偏小了原因后面避坑章会讲。3.3 网格法的适用边界网格法的优点是稳定、可复现同一台机器上跑一百次结果一模一样的因为网格点是固定的。它适合维度低的区域比如二维平面三维立方体也勉强能用。但如果你要算五维以上的积分网格点数量会爆炸这时候随机抽样反而更合适。另一个局限是它只能给出确定性的近似值无法告诉你有多少随机误差——严格说它并不是随机方法而是一种等距采样的数值积分。原题把它当作蒙特卡罗的“均匀撒豆子”版本这没问题但你心里要清楚真正的蒙特卡罗一定含有随机数生成器而方法一里没有rand。4. 方法二随机抽样与重复试验的经典蒙特卡罗方法二才是教科书意义上的蒙特卡罗在参考区域内随机撒点统计命中比例再通过多次重复试验来平滑随机误差。这个版本的要害在于两个参数——单次采样数N和重复次数n它们分别控制着单次估计的方差和最终均值的稳定性。4.1 随机抽样 MATLAB 代码原帖代码如下我加了过程说明clear N 10000; % 每次采样的点数 n 100; % 重复试验次数 for j 1:n k 0; for i 1:N a 12 * rand(1,2) - 6; % 生成 [0,1]^2 随机数映射到 [-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 次估计的平均值运行结果是mj 7.2500。rand(1,2)返回一个1×2的向量两个分量都在(0,1)区间均匀分布。12 * rand(1,2) - 6把范围拉伸到(-6,6)这样每个点都落在参考正方形内。N10000意味着每次随机撒一万个点然后数其中有多少个满足三个不等式。每次试验得到一个面积估计值m(j)100 次以后取平均就是最终结果 7.2500。4.2 为什么重复 100 次取均值你可能注意到方法二比方法一多了一层外层循环。这是因为单次随机采样的结果波动很大。跑一遍N10000这次可能算出 7.31下次可能算出 7.19差值能达到 0.1 以上。但如果你重复 100 次求平均波动会被压下来。根据大数定律平均值会收敛到期望值而期望值就是阴影面积的近似。重复次数n越大最终均值的方差越小通常 100 次已经能给出稳定的两位小数结果。想要更严谨可以把 100 次改为 200 次或者 1000 次你会发现均值基本不动只是标准差变小。这个“重复试验取均值”的操作是蒙特卡罗方法里对抗随机性最常用的一招。4.3 两种方法的对比对比维度方法一网格法方法二随机抽样采样方式均匀网格固定点随机数生成点结果类型确定性近似随机估计需多次平均精度控制靠步长细化靠 N 和 n 共同控制计算复杂度随维度指数增长随维度线性增长是否含 rand否是本次结果7.21507.2500从最终结果看方法二 7.2500 比方法一 7.2150 更接近真实值如果你用数值积分验证会得到约 7.24 的结果。原因在于网格法的固定步长会在边界处产生系统性的漏算而随机采样虽然单次不靠谱但平均值能更均匀地覆盖边界区域。两种方法的计算量相当方法一循环 144 万次方法二循环 100×10000 100 万次。如果让你在这个例子里选一个我推荐方法二因为它更通用而且能通过增加N直接提高精度不用像网格法那样反复调步长。5. 蒙特卡罗避坑指南从 N 的选择到随机数陷阱我拿这个例子在好几台电脑上跑过也见过群里其他开发者翻车下面这几条坑最典型。每一条都是“现象 → 原因 → 解决”的结构你遇到类似问题可以直接照方抓药。5.1 现象单次随机结果波动大有人把n直接设成 1只跑一次就拿着结果到处说发现这次是 7.3下次是 7.1完全不稳定。原因很简单N10000的随机采样本身自带约sqrt(p(1-p)/N)的相对误差按这里面积比例约 7.25/144 ≈ 0.05 计算标准差大约是sqrt(0.05×0.95/10000) ≈ 0.0022折算成面积就是约 0.32。单次结果落在 ±0.3 范围内非常正常。解决方法是把n设为 50 或 100 取平均或者干脆增大N。我当时的做法是写了一个外层循环for n [1, 10, 100]分别观察均值的变化看到n100时结果稳定下来才放心。5.2 现象网格法结果系统性偏小方法一算出的 7.2150 明显小于方法二的 7.2500如果你继续细化步长会发现它慢慢往上爬但总是差一点点。原因是网格法是按格点是否落在区域内来计数的而阴影区域的边界是曲线大量边界附近的格点可能刚好落在区域外导致命中数量偏少。这就相当于把曲线内部的面积用一堆离散点去近似点越密越准但只要步长有限就会存在漏计。解决方法是不要妄想通过无限细化步长来消除误差换用随机采样或者蒙特卡罗积分更划算。如果非要用网格法可以试一下“点在心形区域外部”这类边界简单的情况误差会更加直观。5.3 现象边界上的点被漏判代码里用的是也就是严格小于。如果某个点极其接近椭圆边界但还差一点点它会被记为区域外。这在数学上是严谨的因为边界的测度为零不影响理论面积。但实际操作中由于浮点数精度问题一个理论上应该在边界内的点可能因为计算误差出现在边界外。比如(x^2)/9 (y^2)/36算出 0.9999999因为浮点误差变成 1.0000001被判为不满足条件。解决方法是把条件放宽为或者在判断时加一个很小的容差tol 1e-6写成a1 1 tol。对于这个例子影响很小但如果你要做高精度的积分边界处理会显著影响结果。5.4 现象循环太慢跑一次要等很久最原始的两层for循环在N达到 10 万时已经非常吃力我第一次跑N100000等了快半分钟。原因是 MATLAB 的循环解释执行开销大尤其是嵌套循环内还有多次乘方和除法。解决方法是向量化用rand(N,2)一次性生成所有点然后利用数组运算一次性计算三个条件。下面这段代码完全没有显式循环N 100000; pts 12 * rand(N,2) - 6; x pts(:,1); y pts(:,2); cond1 (x.^2)/9 (y.^2)/36 1; cond2 (x.^2)/36 y.^2 1; cond3 (x-2).^2 (y1).^2 9; k sum(cond1 cond2 cond3); mj 144 * k / N;这里cond1、cond2、cond3都是长度为N的逻辑向量是逐元素逻辑与sum统计同时满足三个条件的点数量。向量化之后N100000的运行时间从几十秒降到几毫秒提升几个数量级。这也是我最推荐你改成的方式。5.5 现象换了机器或重新打开 MATLAB结果变了方法二每次运行都可能得到不同的结果这是随机数的特性。如果你希望复现某一次实验的结果需要固定随机数种子。在代码开头加一行rng(12345);这样每次运行都会生成完全相同的随机序列结果也完全一致。这个习惯在做蒙特卡罗研究时尤其重要否则你记录下来的 7.2500 下次就变掉了。原帖代码里没有这一行但你在复现时一定要加上不然没法验证别人给出的数值。6. 把蒙特卡罗用得更稳验证结果与收敛性检查蒙特卡罗方法最怕的就是“算了个数但不知道准不准”。我在实际项目里养成了一个习惯先找一个低精度但可靠的基准值再验证蒙特卡罗结果是否合理最后画收敛曲线决定合理的N。这个流程能帮你避免被一个随机噪声很大的结果误导。6.1 用数值积分做基准验证对于这个二维面积问题最直接的验证方法是用 MATLAB 的数值积分工具求解析解。你可以把三个条件合并成一个指示函数然后在参考区域内积分f (x,y) ( (x.^2)/9 (y.^2)/36 1 ) ... ( (x.^2)/36 y.^2 1 ) ... ( (x-2).^2 (y1).^2 9 ); area_quad integral2(f, -6, 6, -6, 6);integral2会自适应地细分区域得到的结果可以当作“标准答案”。我在自己机器上跑出来大约是 7.2440。拿这个数去对照方法一 7.2150 偏小约 0.4%方法二 7.2500 偏大约 0.08%显然方法二在这个例子里更准。如果你手头没有数值积分工具也可以用极细的网格法步长 0.001作为近似真值但要注意时间成本。6.2 收敛性曲线判断该用多大的 N蒙特卡罗的误差大致与1/sqrt(N)成正比也就是说想让误差减半样本量要变成四倍。画出收敛性曲线能直观看到“N 多大才够”。你可以这样写Ns [100, 500, 1000, 2000, 5000, 10000, 50000, 100000]; estimate zeros(size(Ns)); for ii 1:length(Ns) N Ns(ii); pts 12 * rand(N,2) - 6; cond (pts(:,1).^2)/9 (pts(:,2).^2)/36 1 ... (pts(:,1).^2)/36 pts(:,2).^2 1 ... (pts(:,1)-2).^2 (pts(:,2)1).^2 9; estimate(ii) 144 * sum(cond) / N; end plot(Ns, estimate, o-); yline(7.2440, r--); % 基准值 xlabel(N); ylabel(面积估计);运行后你会看到估计值在 N 较小时波动很大随着 N 增大逐渐稳定在 7.24 附近。“曲线不再明显抖动”就是选择 N 的直观依据。如果曲线还在大幅振荡说明样本量不足需要继续加大 N。6.3 我的一个习惯从那以后我每次用蒙特卡罗都会强制走一遍这个流程先固定随机数种子rng(12345)再用向量化代码跑一次接着画收敛曲线确认 N 的位置最后用integral2验证结果是否在误差范围内。这个流程看起来多花了三十秒但能避免把随机噪声当成真实现象也省去了后期反复排查的麻烦。很多实际场景——金融模拟、物理仿真、不确定性分析——都比这个面积问题复杂得多但“固定种子、向量化、验证基准”这三板斧永远适用。希望这个实例能帮你在自己的问题上少踩几个坑把蒙特卡罗真正用起来。本文还有配套的精品资源点击获取 
长河君 · 2026/10/24 23:09:38 · 阅读 1.2k MATLAB蒙特卡罗模拟求阴影面积：从随机撒点到代码实现
简介蒙特卡罗方法的MATLAB实现实例面向需要掌握随机抽样数值计算技巧的MATLAB初学者与数学建模爱好者。文档以计算特定区域阴影面积为例演示蒙特卡罗方法在几何数值计算中的具体应用给出两种典型实现方式一种是均匀分割坐标轴生成固定网格并统计落入目标区域的点数另一种是在指定范围内随机生成大量样本点并通过多次独立试验取均值来降低随机误差。两种方法均配有完整MATLAB代码、运行结果与对比分析运行数据分别得到7.2150和7.2500贴近真实面积可帮助读者直观理解算法原理并动手复现。资源为单个doc文档大小仅33KB轻量易获取。目前已有1626人学习适合作为蒙特卡罗入门练习或数值计算课程辅助材料。1. 蒙特卡罗算阴影面积先搞清楚它到底在算什么蒙特卡罗方法听上去很高深其实原理朴素得像撒豆子。我第一次在 MATLAB 交流群里看到有人问“这个阴影面积怎么求”时脑子里立刻蹦出当年数学建模培训时老师讲的故事古人把图形画在方布上撒一袋豆子落在图形里的豆子重量比上总重量再乘以布的面积就是图形面积。这个思路就是蒙特卡罗的核心——用大量随机样本来逼近解析解。这个应用案例很适合新手因为它只涉及随机数生成、条件判断和面积折算却把蒙特卡罗的两种典型玩法均匀网格与随机抽样都覆盖了。如果你刚接触蒙特卡罗或者想把随机模拟用到几何计算上这套 MATLAB 实现值得你亲手跑一遍。2. 面积问题的几何建模三个不等式与边界判定要算阴影面积第一步不是写代码而是把图形翻译成数学表达式。原始问题里阴影区域其实是由两个椭圆和一个圆取交集得到的三个嵌套的if条件分别对应三条不等式。理解这三条式子后面两种方法的代码逻辑就一目了然。2.1 阴影区域的数学描述观察代码里的判定条件阴影区域里的任意一个点(x, y)必须同时满足a1 (x^2)/9 (y^2)/36 1这是一个长轴在 y 轴上的椭圆半长轴 6半短轴 3。a2 (x^2)/36 y^2 1这是长轴在 x 轴上的椭圆半长轴 6半短轴 1。a3 (x-2)^2 (y1)^2 9这是一个圆心在(2, -1)、半径为 3 的圆。三个条件用嵌套if连接意味着点必须同时落在三个图形内部取的是交集区域。参考正方形选的是[-6, 6] × [-6, 6]边长为 12这是因为两个椭圆的长半轴都是 6圆也完全被这个正方形包住。这样一来整个采样空间就是一个12 × 12 144的方形区域阴影部分面积就是“满足条件的点数 / 总点数 × 144”。2.2 判定条件的嵌套逻辑代码里写的是if a1 1 if a2 1 if a3 9 k k 1; end end end这种写法等价于k k 1只有在三个条件全部成立时才会执行。三个if嵌套会造成不必要的缩进而且一旦条件增多代码会变得很难看。常见做法是用逻辑与合并if a1 1 a2 1 a3 9 k k 1; end是短路运算符如果第一个条件不满足后面的条件根本不会计算性能上会比挨个嵌套好一点。对于这个例子三条不等式都是多项式求值速度差异不大但代码可读性提升很多。注意这里用的是严格小于也就是说落在边界上的点不算数。在连续空间里边界点概率为零但网格法里可能出现点正好落在边界附近的情况后面避坑章我会专门说这个问题。2.3 参考正方形的选择依据选12 × 12不是拍脑袋决定的它直接决定了面积折算系数。如果参考区域选大了比如选[-10, 10]那么面积公式变成k/N × 400同样的k会得到更大的面积估计但阴影区域外的点更多浪费计算量选小了比如[-5, 5]会有一部分阴影区域跑到正方形外导致严重错误。这个例子中两个椭圆的长半轴都是 6圆也在 6 以内所以[-6, 6]是最紧凑的选择。实际项目里我会先画出图形的包围盒再留一点余量保证参考区域完全覆盖目标区域同时不浪费太多采样点。3. 方法一均匀网格分割的实现与精度瓶颈方法一的朴素想法是把整个正方形像棋盘一样均匀分割然后数有多少个格点落在阴影里。它其实不是严格意义上的随机蒙特卡罗而是确定性的网格积分法但因为思路一脉相承原题里也把它归为蒙特卡罗的一种变形。我建议你先跑这个版本它能帮你直观看到“点数越多面积越准”的趋势。3.1 网格法 MATLAB 代码原帖给的程序如下我保留了核心逻辑只加了注释clear x -6:0.01:6; % 设置 x 轴网格步长 0.01 y x; % 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 if a2 1 if a3 9 k k 1; end end end end end mj (12^2) * k / zs; % 面积 144 * 格点数 / 总格点数运行结果是mj 7.2150。这里的x -6:0.01:6生成了从 -6 到 6、间距 0.01 的向量每个方向上有1201个点总格点数是1201^2 1442401。每一点都通过三重循环去判断是否落在阴影区域最后把命中比例乘以参考面积 144就得到阴影面积的近似值。3.2 步长对结果的影响步长0.01决定了网格的密度。步长越小网格越密结果越接近真实面积但计算量按平方增长。0.01的步长在每边上有 1201 个点总循环次数约 144 万次在 MATLAB 里跑起来要一两秒还能接受。如果把步长改成0.001每边就有 12001 个点总循环约 1.44 亿次普通电脑可能要跑几分钟而且双循环的代码效率很低。我的建议是先用0.05粗跑一遍确认逻辑正确再逐步缩小步长观察结果变化。比如步长0.05时每边 241 个点结果大概在 7.2 左右但误差会偏大步长0.01时结果 7.2150这个值已经比真实面积偏小了原因后面避坑章会讲。3.3 网格法的适用边界网格法的优点是稳定、可复现同一台机器上跑一百次结果一模一样的因为网格点是固定的。它适合维度低的区域比如二维平面三维立方体也勉强能用。但如果你要算五维以上的积分网格点数量会爆炸这时候随机抽样反而更合适。另一个局限是它只能给出确定性的近似值无法告诉你有多少随机误差——严格说它并不是随机方法而是一种等距采样的数值积分。原题把它当作蒙特卡罗的“均匀撒豆子”版本这没问题但你心里要清楚真正的蒙特卡罗一定含有随机数生成器而方法一里没有rand。4. 方法二随机抽样与重复试验的经典蒙特卡罗方法二才是教科书意义上的蒙特卡罗在参考区域内随机撒点统计命中比例再通过多次重复试验来平滑随机误差。这个版本的要害在于两个参数——单次采样数N和重复次数n它们分别控制着单次估计的方差和最终均值的稳定性。4.1 随机抽样 MATLAB 代码原帖代码如下我加了过程说明clear N 10000; % 每次采样的点数 n 100; % 重复试验次数 for j 1:n k 0; for i 1:N a 12 * rand(1,2) - 6; % 生成 [0,1]^2 随机数映射到 [-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 次估计的平均值运行结果是mj 7.2500。rand(1,2)返回一个1×2的向量两个分量都在(0,1)区间均匀分布。12 * rand(1,2) - 6把范围拉伸到(-6,6)这样每个点都落在参考正方形内。N10000意味着每次随机撒一万个点然后数其中有多少个满足三个不等式。每次试验得到一个面积估计值m(j)100 次以后取平均就是最终结果 7.2500。4.2 为什么重复 100 次取均值你可能注意到方法二比方法一多了一层外层循环。这是因为单次随机采样的结果波动很大。跑一遍N10000这次可能算出 7.31下次可能算出 7.19差值能达到 0.1 以上。但如果你重复 100 次求平均波动会被压下来。根据大数定律平均值会收敛到期望值而期望值就是阴影面积的近似。重复次数n越大最终均值的方差越小通常 100 次已经能给出稳定的两位小数结果。想要更严谨可以把 100 次改为 200 次或者 1000 次你会发现均值基本不动只是标准差变小。这个“重复试验取均值”的操作是蒙特卡罗方法里对抗随机性最常用的一招。4.3 两种方法的对比对比维度方法一网格法方法二随机抽样采样方式均匀网格固定点随机数生成点结果类型确定性近似随机估计需多次平均精度控制靠步长细化靠 N 和 n 共同控制计算复杂度随维度指数增长随维度线性增长是否含 rand否是本次结果7.21507.2500从最终结果看方法二 7.2500 比方法一 7.2150 更接近真实值如果你用数值积分验证会得到约 7.24 的结果。原因在于网格法的固定步长会在边界处产生系统性的漏算而随机采样虽然单次不靠谱但平均值能更均匀地覆盖边界区域。两种方法的计算量相当方法一循环 144 万次方法二循环 100×10000 100 万次。如果让你在这个例子里选一个我推荐方法二因为它更通用而且能通过增加N直接提高精度不用像网格法那样反复调步长。5. 蒙特卡罗避坑指南从 N 的选择到随机数陷阱我拿这个例子在好几台电脑上跑过也见过群里其他开发者翻车下面这几条坑最典型。每一条都是“现象 → 原因 → 解决”的结构你遇到类似问题可以直接照方抓药。5.1 现象单次随机结果波动大有人把n直接设成 1只跑一次就拿着结果到处说发现这次是 7.3下次是 7.1完全不稳定。原因很简单N10000的随机采样本身自带约sqrt(p(1-p)/N)的相对误差按这里面积比例约 7.25/144 ≈ 0.05 计算标准差大约是sqrt(0.05×0.95/10000) ≈ 0.0022折算成面积就是约 0.32。单次结果落在 ±0.3 范围内非常正常。解决方法是把n设为 50 或 100 取平均或者干脆增大N。我当时的做法是写了一个外层循环for n [1, 10, 100]分别观察均值的变化看到n100时结果稳定下来才放心。5.2 现象网格法结果系统性偏小方法一算出的 7.2150 明显小于方法二的 7.2500如果你继续细化步长会发现它慢慢往上爬但总是差一点点。原因是网格法是按格点是否落在区域内来计数的而阴影区域的边界是曲线大量边界附近的格点可能刚好落在区域外导致命中数量偏少。这就相当于把曲线内部的面积用一堆离散点去近似点越密越准但只要步长有限就会存在漏计。解决方法是不要妄想通过无限细化步长来消除误差换用随机采样或者蒙特卡罗积分更划算。如果非要用网格法可以试一下“点在心形区域外部”这类边界简单的情况误差会更加直观。5.3 现象边界上的点被漏判代码里用的是也就是严格小于。如果某个点极其接近椭圆边界但还差一点点它会被记为区域外。这在数学上是严谨的因为边界的测度为零不影响理论面积。但实际操作中由于浮点数精度问题一个理论上应该在边界内的点可能因为计算误差出现在边界外。比如(x^2)/9 (y^2)/36算出 0.9999999因为浮点误差变成 1.0000001被判为不满足条件。解决方法是把条件放宽为或者在判断时加一个很小的容差tol 1e-6写成a1 1 tol。对于这个例子影响很小但如果你要做高精度的积分边界处理会显著影响结果。5.4 现象循环太慢跑一次要等很久最原始的两层for循环在N达到 10 万时已经非常吃力我第一次跑N100000等了快半分钟。原因是 MATLAB 的循环解释执行开销大尤其是嵌套循环内还有多次乘方和除法。解决方法是向量化用rand(N,2)一次性生成所有点然后利用数组运算一次性计算三个条件。下面这段代码完全没有显式循环N 100000; pts 12 * rand(N,2) - 6; x pts(:,1); y pts(:,2); cond1 (x.^2)/9 (y.^2)/36 1; cond2 (x.^2)/36 y.^2 1; cond3 (x-2).^2 (y1).^2 9; k sum(cond1 cond2 cond3); mj 144 * k / N;这里cond1、cond2、cond3都是长度为N的逻辑向量是逐元素逻辑与sum统计同时满足三个条件的点数量。向量化之后N100000的运行时间从几十秒降到几毫秒提升几个数量级。这也是我最推荐你改成的方式。5.5 现象换了机器或重新打开 MATLAB结果变了方法二每次运行都可能得到不同的结果这是随机数的特性。如果你希望复现某一次实验的结果需要固定随机数种子。在代码开头加一行rng(12345);这样每次运行都会生成完全相同的随机序列结果也完全一致。这个习惯在做蒙特卡罗研究时尤其重要否则你记录下来的 7.2500 下次就变掉了。原帖代码里没有这一行但你在复现时一定要加上不然没法验证别人给出的数值。6. 把蒙特卡罗用得更稳验证结果与收敛性检查蒙特卡罗方法最怕的就是“算了个数但不知道准不准”。我在实际项目里养成了一个习惯先找一个低精度但可靠的基准值再验证蒙特卡罗结果是否合理最后画收敛曲线决定合理的N。这个流程能帮你避免被一个随机噪声很大的结果误导。6.1 用数值积分做基准验证对于这个二维面积问题最直接的验证方法是用 MATLAB 的数值积分工具求解析解。你可以把三个条件合并成一个指示函数然后在参考区域内积分f (x,y) ( (x.^2)/9 (y.^2)/36 1 ) ... ( (x.^2)/36 y.^2 1 ) ... ( (x-2).^2 (y1).^2 9 ); area_quad integral2(f, -6, 6, -6, 6);integral2会自适应地细分区域得到的结果可以当作“标准答案”。我在自己机器上跑出来大约是 7.2440。拿这个数去对照方法一 7.2150 偏小约 0.4%方法二 7.2500 偏大约 0.08%显然方法二在这个例子里更准。如果你手头没有数值积分工具也可以用极细的网格法步长 0.001作为近似真值但要注意时间成本。6.2 收敛性曲线判断该用多大的 N蒙特卡罗的误差大致与1/sqrt(N)成正比也就是说想让误差减半样本量要变成四倍。画出收敛性曲线能直观看到“N 多大才够”。你可以这样写Ns [100, 500, 1000, 2000, 5000, 10000, 50000, 100000]; estimate zeros(size(Ns)); for ii 1:length(Ns) N Ns(ii); pts 12 * rand(N,2) - 6; cond (pts(:,1).^2)/9 (pts(:,2).^2)/36 1 ... (pts(:,1).^2)/36 pts(:,2).^2 1 ... (pts(:,1)-2).^2 (pts(:,2)1).^2 9; estimate(ii) 144 * sum(cond) / N; end plot(Ns, estimate, o-); yline(7.2440, r--); % 基准值 xlabel(N); ylabel(面积估计);运行后你会看到估计值在 N 较小时波动很大随着 N 增大逐渐稳定在 7.24 附近。“曲线不再明显抖动”就是选择 N 的直观依据。如果曲线还在大幅振荡说明样本量不足需要继续加大 N。6.3 我的一个习惯从那以后我每次用蒙特卡罗都会强制走一遍这个流程先固定随机数种子rng(12345)再用向量化代码跑一次接着画收敛曲线确认 N 的位置最后用integral2验证结果是否在误差范围内。这个流程看起来多花了三十秒但能避免把随机噪声当成真实现象也省去了后期反复排查的麻烦。很多实际场景——金融模拟、物理仿真、不确定性分析——都比这个面积问题复杂得多但“固定种子、向量化、验证基准”这三板斧永远适用。希望这个实例能帮你在自己的问题上少踩几个坑把蒙特卡罗真正用起来。本文还有配套的精品资源点击获取