MATLAB蒙特卡罗仿真：从随机抽样到可复现的置信区间

简介面向MATLAB初学者与数值计算学习者的蒙特卡罗实例文档围绕随机抽样估算复杂区域面积这一典型问题展开适合用来理解概率统计方法如何落到可运行代码中。压缩包仅含1个doc文件大小约33KB内容紧凑便于直接查阅与课堂演示。相关页面显示已有1627人学习说明该案例在入门练习中具有一定参考价值。文档以计算阴影部分面积为主线给出两种实现思路一是按固定步长均匀划分坐标网格并逐点判定二是用rand在(-6,6)范围内生成大量随机点、重复100次取均值。其中同时呈现嵌套循环、条件判断与向量化写法并对比7.2150与7.2500两种结果及精度差异可帮助读者理解样本量、随机误差与计算效率之间的取舍也能为课程作业和小型数值实验提供参考。1. 蒙特卡罗仿真在 MATLAB 里到底适合算什么从抛硬币到可靠性评估“蒙特卡罗 MATLAB 实现”最容易翻车的地方不是写不出rand而是把随机实验当成了“多跑几次取平均”。一个典型场景是你要评估某零件失效概率解析公式推不动于是用 MATLAB 抽样上万次最后却说不清结果误差有多大。蒙特卡罗适合处理高维积分、随机服务过程、可靠性抽样这类解析解不友好的问题不适合确定性极强、单次仿真成本又高的场景。读者如果是做课程设计、算法验证或工程风险评估先把随机变量、估计量和收敛判据立住再谈代码实例才可复现。2. 把随机问题写成可复现实验蒙特卡罗 MATLAB 的三个建模要素一个能复现的蒙特卡罗实验至少要有三样东西随机变量怎么来、估计量怎么算、误差怎么报。很多人一上来就堆rand(n,1)代码写得很快结果却经不起换一个种子再跑一遍。比较稳的做法是先把问题写成“输入随机变量 X、函数 g(X)、目标量 θE[g(X)]”的形式再决定样本量和置信区间。下面的写法不追求花哨但每一步都能对应到一个可检查的量。2.1 先定随机变量和估计量再碰 rand以估计 π 为例常见做法是取 X~U(-1,1)令 g(X)sqrt(1-X^2)那么 E[g(X)]π/4估计量就是 4 倍样本均值。这个例子小但结构完整随机变量、变换、估计量、标准误都能落地。MATLAB 里先把函数写成一个独立文件后面换参数、换种子、做重复实验都方便。function est mcPiBasic(n, seed) % 估计 pi 的最小蒙特卡罗函数 rng(seed, twister); % 固定随机流便于复现 u rand(n, 1); % U(0,1) x 2*u - 1; % X ~ U(-1,1) g sqrt(1 - x.^2); % 被积函数 est.mean 4 * mean(g); % pi 的估计 est.se 4 * std(g) / sqrt(n); % 标准误 est.ci est.mean [-1, 1] * 1.96 * est.se; end代码里n是样本量seed决定随机流mean(g)是样本均值std(g)/sqrt(n)是标准误的常用估计1.96 对应大样本下的 95% 正态近似置信区间。这里不要急着把n拉到 1e8先用小样本确认估计量方向对不对mcPiBasic(1e4, 1)的结果应该落在 3.14 附近标准误在 0.01 量级。若结果明显偏到 2 或 4先查随机变量范围再查倍数是否写错。标准误不是装饰品它决定了你还要不要继续加样本。蒙特卡罗误差大致按 1/sqrt(n) 下降样本量翻 4 倍标准误才减半。这个规律很朴素但实际项目里经常被忽略有人把n从 1e5 加到 2e5发现置信区间只窄了一点点就以为代码有问题。其实这是正常现象先看目标误差再决定预算。2.2 用固定种子和分层实验骨架做重复实验单次运行只能给出一个估计值不能说明这次估计靠不靠谱。更稳的验证方式是做多次独立重复每次用不同种子记录估计值、标准误和置信区间宽度再看“不同种子之间的波动”和“单次报告的标准误”是否接近。两者接近说明估计量没有明显的实现错误两者差很多就要查随机流、相关性或统计口径。seeds 1:30; N 20000; res zeros(numel(seeds), 3); for k 1:numel(seeds) est mcPiBasic(N, seeds(k)); res(k,:) [est.mean, est.se, est.ci(2)-est.ci(1)]; end mean(res(:,1)) % 30 次重复的均值 std(res(:,1)) % 不同种子之间的波动 mean(res(:,2)) % 单次运行报告的标准误这里seeds是 30 个独立随机流入口N是单次样本量。std(res(:,1))衡量重复实验之间的实际波动mean(res(:,2))是单次标准误的平均。两者在同一量级才正常。若实际波动远大于报告的标准误常见原因是样本不独立、循环里复用了同一段随机数或者把同一个种子反复传给每次实验。这个检查成本很低但能提前挡掉不少“看起来能跑”的假结果。做重复实验时不要把 30 次结果只画一条折线就结束。更实用的做法是同时看均值和区间均值是否稳定、区间是否随 N 增大而收窄、有没有异常种子导致离群。若某一次结果明显偏离先别删掉查清它对应的输入和随机流往往能发现边界条件。2.3 向量化之前先算内存账分块累加的最小写法MATLAB 的向量化很舒服但大样本下内存会先翻车。n1e8时一个 double 列向量就是 800 MBrand、x、g、g.^2同时存在普通机器很容易 Out of memory。比较稳的方式是分块累加每次只生成blockSize个样本累加一阶和二阶和最后统一算均值和方差。这样代码稍长但内存可控结果也可复现。function est mcPiBlock(n, blockSize, seed) % 分块计算 pi避免一次性生成超大数组 rng(seed, twister); sumG 0; sumG2 0; done 0; while done n m min(blockSize, n - done); u rand(m, 1); x 2*u - 1; g sqrt(1 - x.^2); sumG sumG sum(g); sumG2 sumG2 sum(g.^2); done done m; end meanG sumG / n; varG (sumG2 - n * meanG^2) / (n - 1); est.mean 4 * meanG; est.se 4 * sqrt(varG / n); est.ci est.mean [-1, 1] * 1.96 * est.se; endblockSize一般取 1e5 到 1e7看机器内存和单块计算量sumG、sumG2用 double 累加避免单精度误差积累varG用一阶、二阶和反推样本方差。这个骨架可以直接套到其他估计量上只要把g换成你的目标函数把4*meanG换成对应倍数即可。分块写法还有一个好处块大小固定后结果不随机器内存变化而漂移复现时更安心。3. 用 MATLAB 跑通三类最小实例积分、排队与可靠性蒙特卡罗的实例很多但真正能覆盖大多数落地场景的无非三类积分估计、随机过程仿真、失效概率评估。下面三个例子都尽量小代码可以直接复制运行参数也能按需改。每个例子都保留标准误或置信区间避免只给一个孤零零的均值。3.1 蒙特卡罗积分估计 π 与定积分的可复现脚本先看一个比 π 更通用的定积分例子估计 ∫₀¹ exp(-x²) dx。这个积分没有初等原函数适合用蒙特卡罗做教学和验证。取 X~U(0,1)估计量是(b-a)*mean(f(X))标准误也按同样倍数缩放。f (x) exp(-x.^2); a 0; b 1; n 1e6; seed 7; rng(seed, twister); x a (b - a) * rand(n, 1); y f(x); Ihat (b - a) * mean(y); se (b - a) * std(y) / sqrt(n); ci Ihat [-1, 1] * 1.96 * se; fprintf(Ihat%.6f, se%.3e, CI[%.6f, %.6f]\n, Ihat, se, ci(1), ci(2));a、b是积分区间n是样本量x是从区间均匀抽样得到的点y是对应函数值Ihat是积分估计se是标准误。这个积分真值约 0.746824运行结果应落在附近。若区间换成[0,2]只需改a、b(b-a)会同步缩放不要手动再乘错。对于高维积分思路一样只是rand的列数变成维度内存和样本量要重新算账。3.2 离散事件排队M/M/1 等待时间的蒙特卡罗近似排队仿真是蒙特卡罗很典型的落地场景。以 M/M/1 为例到达间隔服从指数分布服务时间也服从指数分布利用率 ρλ/μ 必须小于 1。代码里用 Lindley 递推算每个顾客的开始服务时间和等待时间再去掉前一段 burn-in降低初始状态影响。lambda 0.8; mu 1.0; n 2e5; burn 2e4; seed 12; rng(seed, twister); A exprnd(1/lambda, n, 1); % 到达间隔均值 1/lambda S exprnd(1/mu, n, 1); % 服务时间均值 1/mu arrival cumsum(A); start zeros(n,1); start(1) arrival(1); for i 2:n start(i) max(arrival(i), start(i-1) S(i-1)); end Wq start - arrival; Wq Wq(burn1:end); WqMean mean(Wq); WqSE std(Wq) / sqrt(numel(Wq)); theory lambda / (mu * (mu - lambda)); fprintf(WqMean%.4f, naiveSE%.4f, theory%.4f\n, WqMean, WqSE, theory);exprnd的参数是均值不是速率这一点很容易混lambda0.8表示每小时平均到达 0.8 个平均到达间隔是 1/0.81.25 小时。burn是预热样本数去掉初始空队列带来的偏差theory是 M/M/1 的平均排队等待时间理论值。这里报告的是 naive 标准误实际排队输出有自相关标准误会偏窄后续要用批均值法修正。先把均值和理论值对上再谈区间。3.3 可靠性评估失效抽样与置信区间可靠性里常见的问题是抗力 R 和荷载 S 都有随机性失效事件是 RS求失效概率 Pf。解析方法推不动时直接抽样是最直观的入口。下面用两个正态分布做演示fail是 0/1 指示变量Pf 的估计就是失效率标准误按伯努利公式算。n 1e6; seed 21; rng(seed, twister); R 10 1.5 * randn(n, 1); % 抗力 S 6 1.0 * randn(n, 1); % 荷载 fail R S; pf mean(fail); se sqrt(pf * (1 - pf) / n); ci pf [-1, 1] * 1.96 * se; fprintf(pf%.5f, se%.3e, CI[%.5f, %.5f]\n, pf, se, ci(1), ci(2));R、S的均值和标准差按你的工程背景改fail是逻辑向量mean(fail)直接给出失效比例se是二项比例的标准误。这里有一个边界如果 Pf 很小比如 1e-6而 n1e6可能只抽到个位数失效样本置信区间会很不稳。此时要么加样本量要么用重要性抽样、分层抽样等方差缩减手段。先把小例子跑通再决定要不要上这些技巧。4. 参数怎么设才不玄学样本量、方差缩减与收敛判据参数设置是蒙特卡罗 MATLAB 实现里最容易“凭感觉”的部分。样本量拍脑袋、方差缩减随便套、收敛判据只看一条曲线最后结果可能很漂亮但不可信。比较稳的顺序是先用试点样本估变异系数再按目标相对误差算样本量能用对偶变量、控制变量、分层抽样降低方差就先用最后用批均值或独立重复检查收敛。4.1 样本量 N 的三种估算口径标准误、相对误差与预算样本量不是越大越好而是要和目标误差匹配。若目标是相对标准误 εr常用近似是 N≈(z·CV/εr)²其中 CVσ/μ 是变异系数z 取 1.96 对应 95% 置信。这个公式要求均值不为零且试点样本能代表总体。下面给一张速查表方便先估预算。目标相对误差CV0.2CV0.5CV1.010%16973855%6238515371%1537960438416表里的数只是量级参考实际还要考虑单次仿真成本和内存。若单次抽样很贵可以先跑 1000 到 5000 个试点样本估计 CV 后再决定总量。下面这个函数把试点估计和样本量计算串起来。function n sampleSizePilot(pilot, errRel, alpha) % 用试点样本估计所需样本量errRel 是相对标准误目标 z norminv(1 - alpha/2); cv std(pilot) / mean(pilot); n ceil((z * cv / errRel)^2); endpilot是试点样本的输出值不是输入随机数errRel是目标相对误差比如 0.05 表示 5%alpha一般取 0.05。这个函数假设试点样本独立同分布且均值远离零。如果mean(pilot)接近零相对误差会爆炸这时改用绝对误差或直接报告标准误。样本量算出来后再乘一个安全系数常见做法是 1.2 到 1.5防止试点 CV 偏小。4.2 方差缩减对偶变量、控制变量和分层抽样的适用边界方差缩减不是万能药但对偶变量在单调函数上很好用。以 π 为例取 u 和 1-u 生成一对负相关样本再对配对结果取平均。这样有效样本量还是 2n但估计量方差可能明显下降。代码里用pairAvg做统计单位标准误按配对数量算不要按总样本数算。nPairs 5e5; seed 31; rng(seed, twister); u rand(nPairs, 1); x1 2*u - 1; x2 -x1; g1 sqrt(1 - x1.^2); g2 sqrt(1 - x2.^2); pairAvg 0.5 * (g1 g2); theta 4 * mean(pairAvg); se 4 * std(pairAvg) / sqrt(nPairs); ci theta [-1, 1] * 1.96 * se; fprintf(theta%.6f, se%.3e\n, theta, se);nPairs是配对数量总抽样次数是 2*nPairsx2-x1构造对偶样本pairAvg是每对的平均se用std(pairAvg)/sqrt(nPairs)。对偶变量适合 g 单调或近似单调的场景若函数在区间内剧烈振荡负相关可能变弱甚至消失。控制变量则适合找一个与 g 相关、均值已知的辅助变量分层抽样适合输入分布可以自然分层的场景。选哪种先看问题结构不要为了“高级”硬套。4.3 收敛判据与批均值法别只看一条曲线收敛判据不能只看估计值曲线“看起来平了”。更实用的做法是看标准误是否稳定、批均值是否漂移、不同独立重复是否给出一致区间。对排队这类有自相关的输出批均值法很常用把长序列切成若干批先算批均值再用批均值算标准误。批数一般取 20 到 50批太长会导致批数太少批太短又去不掉自相关。function [mb, seb] batchMeans(x, k) % 把相关样本切成 k 批用批均值估计均值和标准误 x x(:); n numel(x); b floor(n / k); x x(1:b*k); X reshape(x, b, k); mb mean(X, 1); seb std(mb) / sqrt(k); endx是相关输出序列k是批数b是每批长度mb是总均值seb是批均值标准误。用的时候先对Wq去掉 burn-in再调用batchMeans(Wq, 50)。如果批均值随批次有明显趋势说明 burn-in 不够或系统还没进入稳态如果seb随k变化很大说明批长选择还不稳需要加长序列再试。这个检查比盯一条曲线靠谱得多。5. 常见踩坑排查随机种子、向量化内存与统计口径的 5 个血泪教训蒙特卡罗代码跑起来不难难的是结果可信。下面五条都是实际写 MATLAB 时容易遇到的坑每条按“现象、原因、解决”拆开。它们不针对某个特定项目换一个仿真场景也大概率会撞上。5.1 随机种子设了但结果还变并行与函数作用域现象主脚本里写了rng(42)换成parfor后每次汇总结果仍有差异。原因并行 worker 有各自独立的随机流主进程的种子不会自动同步到每个 worker函数内部的rng也可能影响调用方的随机状态。解决给每个 worker 显式分配独立子流不要依赖全局种子。nWorkers 4; out zeros(1, nWorkers); streams RandStream.create(mrg32k3a, NumStreams, nWorkers, ... StreamIndices, 1:nWorkers, Seed, 42); parfor i 1:nWorkers rng(streams(i)); % 每个 worker 使用独立子流 u rand(1000, 1); out(i) mean(u); endRandStream.create生成互不重叠的子流StreamIndices指定每个 worker 用哪一路rng(streams(i))在 worker 内设置随机源。这样并行结果可复现也不会因为 worker 调度顺序不同而漂移。若不用并行也建议把rng放在函数入口而不是散落在脚本中间。5.2 向量化后内存炸掉分块与单精度陷阱现象n1e8时 MATLAB 卡死或报 Out of memory。原因rand、x、g、g.^2同时存在多个 double 数组叠加内存迅速翻倍。解决分块处理临时数组可以用 single 降内存但累加器保持 double。function [m, se] piBlockSingle(n, blockSize, seed) rng(seed, twister); sumG 0; sumG2 0; done 0; while done n m min(blockSize, n - done); u rand(m, 1, single); x 2*u - 1; g sqrt(1 - x.^2); gd double(g); % 累加用 double避免误差累积 sumG sumG sum(gd); sumG2 sumG2 sum(gd.^2); done done m; end meanG sumG / n; varG (sumG2 - n*meanG^2) / (n - 1); m 4*meanG; se 4*sqrt(varG/n); endblockSize取 1e5 到 1e6 比较稳rand(...,single)只降低临时数组内存sumG、sumG2仍用 doublegd显式转换避免单精度累加导致尾数误差。分块代码比纯向量化长但内存可预测尤其适合在普通工作站上跑大样本。5.3 置信区间算窄了独立同分布假设被破坏现象排队等待时间的 95% 置信区间很窄但换一个种子重跑均值偏移明显超过区间宽度。原因Wq序列有自相关 naive 标准误把样本当成独立同分布低估了真实波动。解决用批均值法或独立重复实验别直接对整段序列套std/sqrt(n)。[mb, seb] batchMeans(Wq, 50); batchCI mb [-1, 1] * 1.96 * seb; naiveCI mean(Wq) [-1, 1] * 1.96 * std(Wq) / sqrt(numel(Wq)); fprintf(batch CI[%.4f, %.4f], naive CI[%.4f, %.4f]\n, ... batchCI(1), batchCI(2), naiveCI(1), naiveCI(2));batchMeans用前面的函数k50是常用起点batchCI通常比naiveCI宽这才更接近真实不确定性。若批均值仍有趋势先加 burn-in若批均值方差随批数变化很大增加总样本量或调整批长。置信区间不是越窄越好窄而不可信比宽更危险。5.4 结果“看起来对”但量纲错rate 与 mean 混用现象M/M/1 的平均等待时间比理论值小很多或者算出来的利用率大于 1。原因exprnd的参数是均值不是速率λ0.8 每小时表示平均到达间隔 1.25 小时不是 0.8 小时。解决把速率和均值分开命名先检查 ρλ/μ 是否小于 1。lambda 0.8; mu 1.0; % 每小时 meanInterarrival 1/lambda; meanService 1/mu; A exprnd(meanInterarrival, n, 1); S exprnd(meanService, n, 1); rho lambda / mu; assert(rho 1, M/M/1 需要 rho1);meanInterarrival、meanService用显式变量名减少把速率直接塞进exprnd的概率assert在参数不合理时直接停不要等跑完才发现利用率大于 1。单位也要统一小时、分钟、秒不能混着进同一个模型。5.5 图表复现不了保存随机流与输入现象昨天跑出来的图很平滑今天重跑形状变了但代码看起来没改。原因只保存了图片没有保存 seed、样本量、分块大小、分布参数和随机流类型。解决每次关键运行保存一个运行清单至少包含 seed、n、参数和输出摘要。runInfo struct(seed, 42, n, 1e6, blockSize, 1e5, ... lambda, 0.8, mu, 1.0, rngType, twister); save(mc_run_info.mat, runInfo);runInfo不用很大但要把能影响结果的输入都记下来save存成 MAT 文件和结果图放在同一目录。后面复现时先加载清单再按清单里的参数重跑。这个习惯不显眼但能省掉很多“黑匣子”式的排查时间。6. 从能跑到可信用对偶变量与置信区间做交叉验证的进阶技巧代码能跑之后下一步不是立刻加样本量而是做一次低成本交叉验证。我一般会用同一预算对比朴素抽样和对偶变量如果两者估计接近但对偶变量的标准误明显更小说明方差缩减起了作用如果两者差得很远先查抽样逻辑而不是急着相信区间更窄的那个。n 1e6; seed 99; rng(seed, twister); uNaive rand(n,1); gNaive sqrt(1 - (2*uNaive-1).^2); thetaNaive 4*mean(gNaive); seNaive 4*std(gNaive)/sqrt(n); uPair rand(n/2,1); x 2*uPair - 1; g1 sqrt(1 - x.^2); g2 sqrt(1 - (-x).^2); pairAvg 0.5*(g1 g2); thetaAnti 4*mean(pairAvg); seAnti 4*std(pairAvg)/sqrt(n/2); fprintf(naive: %.6f (se %.3e)\n, thetaNaive, seNaive); fprintf(anti : %.6f (se %.3e)\n, thetaAnti, seAnti);n是总抽样预算朴素法用 n 个独立样本对偶法用 n/2 对样本总抽样次数相同seNaive和seAnti分别对应两种估计的标准误。看结果时先看两个均值是否本文还有配套的精品资源点击获取