蒙特卡洛模拟从原理到实战：工期、库存与风险估算

上周某位做项目管理的朋友问了我一个问题：手上这个研发项目，合同里写死了六个月交付，可每个环节的工期都是估出来的，有的乐观有的悲观，最后按期完成的概率到底有多大？他手上有一张 Excel 表，每个任务的工期给了三个数：最短、最可能、最长。我说这张表其实已经够了，缺的只是把这三个数变成概率分布，然后让计算机替你把这套流程"跑"十万遍。这就是蒙特卡洛模拟最朴素的用法——当你没法用一个公式直接算出结果时，就让随机性替你把所有可能性走一遍，最后数一数有多少次落在了你想要的那个区间里。

这篇文章我想把蒙特卡洛模拟这件事从头到尾讲透：它到底在算什么、为什么能算对、随机数和概率分布这些零件该怎么选、四个从易到难的完整案例怎么落地、误差为什么会抖、以及在真机上跑的时候我踩过哪些坑。不管你是做项目管理、供应链、财务测算、工程设计，还是单纯想给自己写的模拟程序找一个靠谱的骨架，下面这些内容都能直接拿去用。代码全部用 Python，依赖只有 NumPy，复制粘贴就能跑。

1. 先搞清楚蒙特卡洛到底在算什么

1.1 用撒豆子测面积，比公式更直观

想象你在一张白纸上画了一个不规则的水洼，你想知道它的面积，但形状太怪，积分积不出来。这时候有个笨办法：画一个边长为 2 的正方形把水洼整个框住，然后闭着眼睛往这个正方形里随机撒一万颗豆子，数一数有多少颗落在了水洼里面。落进去的比例乘以正方形面积，就是水洼面积的估计值。撒的豆子越多，估计越准。

这个笨办法就是蒙特卡洛模拟的全部内核。它不关心水洼的边界曲线有多复杂，只关心一件事：随便丢一个点，它落在里面的概率是多少。而这个概率恰好等于面积之比。把"面积"换成"项目按期完成的概率"、"库存不够用的概率"、"投资组合亏损超过某个额度的概率"，思路完全一样——把问题转化成一个"随机试验中事件发生的频率"，然后用大量重复试验去逼近这个频率。

我特别喜欢这个类比的原因在于，它把蒙特卡洛最大的优势和最大的软肋同时暴露了出来。优势是它对问题形状完全不敏感，边界再奇怪也能算；软肋是它靠"数数"得到答案，精度提升得非常慢，你想让误差缩小一半，豆子得撒四倍。后面讲收敛性的时候我们会把这个"四倍"精确地算出来。

1.2 大数定律撑着，但撑得并不轻松

从数学上讲，蒙特卡洛模拟合法性的根基是大数定律：当独立重复试验的次数 n 趋于无穷时，事件发生的频率依概率收敛于真实概率。这句话翻译成人话就是——只要你撒得够多，频率总会靠拢概率。但"总会靠拢"这四个字里藏着一个让人不太舒服的事实：它没有告诉你多快。

实际工程里我们靠的是中心极限定理给出的近似：估计值的误差大致服从正态分布，标准差等于单次试验结果的标准差除以根号 n。这个公式是整篇文章后面所有优化手段的出发点。举个例子，假设单次试验的结果标准差是 0.5，你跑了一万次，那么估计值的标准误就是 0.5 / 100 = 0.005。想把这个 0.005 压到 0.0005，n 得从一万涨到一百万，计算量翻一百倍，精度只提升十倍。

注意：很多人第一次做模拟，看到结果小数点后第三位还在跳，第一反应是"程序有 bug"。其实不是，那就是蒙特卡洛的正常呼吸。判断有没有 bug 的方法不是看结果稳不稳，而是看结果有没有系统性地偏离预期，以及标准误是否和理论公式对得上。

1.3 什么场景该用它，什么场景别碰

蒙特卡洛不是万能锤子。我的判断标准很简单，满足下面任意一条就值得上：

• 问题里有多个不确定量在互相叠加，比如工期、成本、需求、汇率同时波动，解析解根本写不出来。
• 目标是一个概率或者分位数，而不是一个确定值。比如"95% 的情况下亏损不超过多少"这种问题，传统公式只能给你均值，给不了尾部。
• 系统里有非线性或者路径依赖，比如某个合同条款规定"超过预算 10% 的部分由乙方承担"，这种分段函数用解析方法处理起来极其痛苦。
• 维度很高但每一维都很容易采样，这一点是蒙特卡洛相对网格法的绝对优势，后面会专门讲。

反过来，如果问题本身有解析解，或者维度只有一两维、函数又很光滑，那老老实实用数值积分或者查表，没必要杀鸡用牛刀。我见过有人用十万次模拟去算一个标准正态分布的 95% 分位数，算出来 1.6448，和查表得到的 1.645 差 0.0002，白白烧了 CPU。工具选型的第一步永远是问一句：这个问题真的需要模拟吗？

2. 搭一个能跑的最小系统：核心零件拆解

2.1 随机数发生器：整个链条的质量源头

所有蒙特卡洛模拟的地基是随机数。但计算机生成的是伪随机数，由确定性算法产生，只要种子相同，序列就完全相同。这不是缺陷，反而是工程上极其重要的特性——它保证了结果可复现。你调参调到半夜发现结果对不上，第一件事就是回头检查种子有没有固定住。

Python 里现在推荐用numpy.random.default_rng()，底层是 PCG64 算法，周期长、统计性质好、支持并行时的独立流。老代码里常见的np.random.seed()走的是 Mersenne Twister，虽然也能用，但在并行场景下容易出问题，新项目我基本不再用了。

import numpy as np # 推荐写法：显式构造生成器，种子固定 rng = np.random.default_rng(20240607) # 生成一万个 [0,1) 均匀随机数 u = rng.random(10_000) # 生成一万个标准正态随机数 z = rng.standard_normal(10_000) print(u.mean(), u.std()) # 期望接近 0.5 和 0.2887 print(z.mean(), z.std()) # 期望接近 0 和 1

这段代码看着平淡，但它是你每次动手前的"体检"。跑一下，看看均值和标准差是不是在合理范围内。如果u.mean()出来是 0.9，那后面所有结果都是废的，原因可能是你误用了别的函数或者种子被覆盖了。这个习惯帮我省过好几次时间。

还有一个细节容易被忽略：不同维度的随机数要用同一次调用生成。比如你要模拟一百个任务各自的工期，正确的做法是一次性rng.standard_normal((n_sim, 100))生成一个矩阵，而不是在循环里每次调用一次。前者速度快两个数量级，后者在 Python 里循环一百万次光函数调用开销就够你喝一壶。

2.2 概率分布的建模与参数选择

蒙特卡洛的输出质量，九成取决于输入分布的合理性。这是整个流程里最需要业务判断、最没法自动化的一步。我的经验是按"你手上有多少信息"来选分布：

三角分布是我在项目估算里用得最多的。原因很实在：专家不愿意给你一个精确的均值和方差，但愿意说"最快两个月，最慢五个月，最可能三个月"。这三个数直接对应三角分布的三个参数，中间那个"最可能"就是众数。PERT 分布是三角分布的改良版，它在众数附近更平滑，尾部也不会像三角分布那样突然截断，做风险分析时我更倾向用它。

这里有个坑必须提醒：三角分布的均值不等于"最可能值"。三个参数是 (2, 3, 5) 时，均值是 (2+3+5)/3 = 3.33，而不是 3。如果你拿"最可能值"当均值去汇报，整体结果会系统性偏低。这个偏差在任务链很长的时候会被放大，我自己就吃过一次亏，最后给出去的完工日期比实际乐观了将近两周。

对数正态的参数换算也是个高频错误点。它的两个参数是对数域的均值和标准差，不是原始尺度的。假设你观察到某个成本项的均值是 100 万、标准差是 30 万，想用对数正态拟合，得先换算：

mu_orig, sigma_orig = 100.0, 30.0 sigma_ln = np.sqrt(np.log(1 + (sigma_orig / mu_orig) ** 2)) mu_ln = np.log(mu_orig) - 0.5 * sigma_ln ** 2 samples = rng.lognormal(mean=mu_ln, sigma=sigma_ln, size=100_000) print(samples.mean(), samples.std()) # 应接近 100 和 30

直接把 100 和 30 塞进lognormal是个经典错误，出来的分布完全不是你想的那样，均值会跑到几百上去。写完分布拟合，一定养成回头验证一下样本均值和标准差的习惯，花五秒钟，省半天排查。

2.3 结果统计：别只盯着均值

跑完十万次模拟，很多人的第一反应是取个平均值就交差了。这浪费了模拟能提供的绝大部分信息。真正的价值在于分布的形状和尾部。

我通常固定输出这几样东西：

• 均值和中位数：两者差距大，说明分布偏斜严重，这时候用均值做决策是有风险的。
• 分位数：P10、P50、P90、P95，这是给管理层看的最有用的数字。
• 超过某个阈值的概率：比如"超预算的概率"、"延期超过一个月的概率"。
• 标准误：这个数字告诉你结果本身的抖动幅度，不给标准误的模拟结果是不严谨的。

def summarize(samples, threshold=None): s = np.asarray(samples) n = s.size out = { "n": n, "mean": s.mean(), "median": np.median(s), "std": s.std(ddof=1), "se": s.std(ddof=1) / np.sqrt(n), "P10": np.percentile(s, 10), "P50": np.percentile(s, 50), "P90": np.percentile(s, 90), } if threshold is not None: out["P_exceed"] = (s > threshold).mean() # 比例的标准误 p = out["P_exceed"] out["P_exceed_se"] = np.sqrt(p * (1 - p) / n) return out

注意P_exceed那一行的标准误算法。估计一个概率时，它的标准误是 √(p(1-p)/n)，和我们估计均值用的公式不一样。如果你要汇报"延期概率是 12%"，顺手带上"±0.3%"，可信度立刻不一样。这个细节我在给决策层做汇报时反复强调，因为他们很容易把 12% 和 13% 当成实质性差异，其实完全在噪声范围内。

3. 四个从易到难的完整实操案例

3.1 案例一：用撒点法估圆周率，验证框架是否正确

这是所有蒙特卡洛教程的第一课，但我不建议跳过它。它的作用是给你一个已知答案的沙盒，用来确认你的随机数、循环、统计这一整套流程没问题。在正方形 [0,1)² 里均匀撒点，落在单位圆四分之一内的比例乘以 4 就是 π 的估计。

def estimate_pi(n, seed=42): rng = np.random.default_rng(seed) x = rng.random(n) y = rng.random(n) inside = (x * x + y * y) <= 1.0 return 4.0 * inside.mean(), inside est, mask = estimate_pi(1_000_000) print(f"估计值: {est:.6f}") # 真值 3.141593

跑一百万次，结果大概落在 3.140 到 3.143 之间。现在算一下理论误差：单次试验的结果是一个 0/1 变量，取 1 的概率是 π/4 ≈ 0.7854，所以方差是 0.7854 × 0.2146 ≈ 0.1685，标准差约 0.4105。乘以 4 之后标准差是 1.642。标准误 = 1.642 / √1000000 = 0.00164。也就是说 95% 置信区间大约是 ±0.0032。这个数字和实测的偏差完全对得上。

这里可以顺手做一个样本量的反向推算，这是实操中真正有用的技能。假设你要求 π 的估计误差不超过 0.001（95% 置信），需要多少样本？

n = (1.96 × 1.642 / 0.001)² ≈ 1.06 × 10⁷

一千万次。你看，为了把误差从 0.0032 压到 0.001，样本量涨了十倍多。这个成本意识一定要建立起来，否则很容易写出跑了三个小时还没收敛的脚本。实际做项目的时候，我一般先跑一万次看看数量级，再根据目标精度倒推需要多少，而不是一上来就设个一千万的无脑参数。

3.2 案例二：项目工期的完工概率模拟

回到开头那位朋友的问题。假设一个项目拆成五个串联的关键任务，每个任务的工期用三角分布描述，单位是天。同时还有两个任务存在"返工风险"——有 20% 的概率需要额外增加一段工期。这种带条件分支的结构，用公式几乎没法算，但用模拟就是几行代码。

def simulate_project(n_sim=200_000, seed=7): rng = np.random.default_rng(seed) # 五个串联任务：(最短, 最可能, 最长) tasks = [ (20, 30, 55), (15, 22, 40), (25, 35, 62), (10, 16, 30), (30, 45, 80), ] total = np.zeros(n_sim) for low, mode, high in tasks: total += rng.triangular(low, mode, high, size=n_sim) # 两个独立返工事件：20% 概率，额外工期用指数分布 for p_rework, mean_extra in [(0.20, 12.0), (0.15, 8.0)]: hit = rng.random(n_sim) < p_rework extra = rng.exponential(mean_extra, size=n_sim) total += np.where(hit, extra, 0.0) return total days = simulate_project() stats = summarize(days, threshold=180) for k, v in stats.items(): print(f"{k}: {v:.4f}")

跑完你会拿到一条完整的工期分布曲线。我在实际项目里最常被问到的三个问题，这套模拟能一次全答上：按期（180 天）完成的概率是多少、如果要保证 90% 的把握，合同该签多少天、最坏情况下（P95）要准备多少缓冲。

第二个问题尤其值得说。做法是把样本排序，取 90% 分位数，那个数字就是"有九成把握能完成"的工期。这比拍脑袋加 20% 缓冲科学得多，因为它是从不确定性结构里推出来的，而不是拍出来的。我给某制造企业做过一次类似的测算，他们原来习惯在所有任务上加 15% 缓冲，模拟出来的结果显示，这个统一比例在不同任务链上的实际覆盖水平从 62% 到 94% 不等，等于有的地方浪费、有的地方裸奔。

还有一个实操技巧：把敏感度分析一起做掉。做法很简单，固定其他任务的参数，只把某一个任务的"最长工期"从 55 调到 70，看最终分布的 P90 变化了多少。变化最大的那个任务，就是你应该重点去盯、去压缩不确定性的对象。这个信息在资源有限的时候极其值钱，比平均分配管理精力有效得多。

3.3 案例三：安全库存水平和缺货概率的确定

这个案例换个领域，讲供应链。某仓库面对的日需求服从正态分布，均值 120 件，标准差 30 件；补货提前期也不是固定的，均值 7 天，标准差 2 天。问题是：订货点设在多少，才能把缺货概率控制在 5% 以内？

简化公式是 ROP = μ_d × L + z × σ_d × √L，但注意这个公式假设提前期是固定的。一旦提前期本身也随机，公式就变成近似，误差在参数波动大时相当可观。这时候蒙特卡洛的优势就出来了。

def simulate_rop(rop, n_sim=300_000, seed=11): rng = np.random.default_rng(seed) # 每个模拟周期：提前期内的总需求 lead_time = rng.normal(7.0, 2.0, size=n_sim) lead_time = np.clip(lead_time, 1.0, None) # 提前期不能为负 # 提前期内需求 = 逐日需求累加，用泊松过程近似更贴近现实 daily_mean = 120.0 # 逐日累加：先生成足够大的矩阵再掩码求和 max_days = int(np.ceil(lead_time.max())) + 1 daily = rng.normal(120.0, 30.0, size=(n_sim, max_days)) daily = np.clip(daily, 0, None) day_idx = np.arange(max_days)[None, :] mask = day_idx < lead_time[:, None].astype(int) demand = (daily * mask).sum(axis=1) stockout_prob = (demand > rop).mean() return stockout_prob, demand for rop in [900, 1000, 1100, 1200, 1300]: p, _ = simulate_rop(rop) print(f"ROP={rop} 缺货概率={p:.4f}")

这个脚本里有两个细节值得单拎出来说。第一是np.clip(lead_time, 1.0, None)，正态分布会生成负的提前期，虽然概率很小，但在三十万次模拟里总会碰到几个，一旦出现负值，下面的掩码求和就会算错。这类"物理上不可能但数学上会生成"的值，是模拟代码里最隐蔽的 bug 来源。第二是逐日累加的写法用了矩阵掩码，而不是在 Python 里写循环，速度差了一个数量级。

拿这个结果去和简化的解析公式对比，你会发现当提前期波动较大的时候，蒙特卡洛给出的缺货概率明显高于公式值。原因是公式把提前期的波动"吸收"进了需求的标准差里，但那个近似在提前期变异系数超过 0.3 之后就开始失真。这也是我在实际项目里更愿意用模拟的原因——它不需要你相信任何一个近似假设。

3.4 案例四：风险敞口的尾部估算

第四个案例往金融方向走一步，讲风险度量。假设某组合的日收益率近似服从一个均值 0.03%、标准差 1.2% 的分布，但真实收益率的尾部比正态更肥。如果我们直接用正态假设去算 99% 分位数的损失，会系统性低估风险。用蒙特卡洛配合 t 分布采样，可以把肥尾的特征带进来。

def simulate_var(n_sim=500_000, seed=2024, df=4): rng = np.random.default_rng(seed) mu, sigma = 0.0003, 0.012 # t 分布标准化后缩放到目标均值方差 raw = rng.standard_t(df, size=n_sim) raw = (raw - 0) / np.sqrt(df / (df - 2)) # 方差归一化到 1 ret = mu + sigma * raw loss = -ret var_99 = np.percentile(loss, 99) es_99 = loss[loss >= var_99].mean() # 期望尾部损失 return var_99, es_99, ret var99, es99, ret = simulate_var() print(f"99% VaR = {var99:.4%}, 99% ES = {es99:.4%}")

对比一下同一组参数下正态假设给出的 VaR：1.96 × 1.2% - 0.03% ≈ 2.32%。而 t 分布在自由度 4 时给出的结果通常明显更高，因为尾部更厚。这个差距在平静市场里看不出来，在极端行情下就是生死之别。ES（期望尾部损失）这个指标比 VaR 更有参考价值，因为它回答的是"一旦突破，平均亏多少"，而不是"跌破这条线的概率是 1%"。

提示：做尾部估计时，样本量必须成倍放大。因为 99% 分位数意味着你只用了最极端的 1% 样本，五十万次里也就五千个点。要估 99.9% 分位数，样本量不上千万根本站不住。这是尾部模拟最容易翻车的地方，很多人拿一万次模拟去报 99.9% VaR，那个数字基本是随机数。

4. 收敛太慢、结果抖动？误差控制与加速手段

4.1 误差量级的数学来源，以及为什么高维反而占优

前面反复提到 1/√n，这里把它的来源说清楚。设单次模拟的结果是随机变量 X，其真实均值为 μ，标准差为 σ。跑 n 次得到的样本均值 X̄，其标准差就是 σ/√n。这个式子告诉我们两件事：误差随 n 的平方根倒数下降，以及误差和问题的维度完全无关。

第二点经常被误解。很多人听说"维度灾难"，就以为蒙特卡洛在高维下也会崩。事实恰恰相反：如果要在 d 维空间做确定性网格积分，每个维度切 m 个点，总点数就是 m^d，d 稍微大一点就爆炸；而蒙特卡洛只关心采样次数 n，不管是 2 维还是 200 维，1/√n 的收敛速度一模一样。这就是它在高维积分、多因子风险模型里不可替代的原因。当然代价是它收敛得比低维的确定性方法慢——一维积分用辛普森法则误差是 O(n⁻⁴)，蒙特卡洛完全打不过。

理解了这个权衡，工具选型就有依据了：低维、光滑、有解析结构的问题，用数值积分；高维、不规则、带跳变的问题，用蒙特卡洛。

4.2 方差缩减：不动样本量也能提精度

既然误差是 σ/√n，那除了加大 n，还有一条路——把 σ 压小。这就是方差缩减技术的全部出发点。我常用的有四种，按上手难度排序：

对偶变量法最容易实现，代码改造只要两行：

def antithetic_normal(n, rng): half = rng.standard_normal(n // 2) return np.concatenate([half, -half])

把生成的正态随机数取一半，另一半直接取负，然后用拼接起来的两倍长度去做模拟。因为 -Z 和 Z 的均值都是 0、方差都是 1，但两者负相关，合并之后估计量的方差会显著下降。前提是被积函数得是单调的，如果函数有对称性或者振荡，这招可能完全失效甚至帮倒忙。

控制变量法威力更大但需要动脑子。思路是找一个和你的目标量高度相关、且均值已知的量 Y，然后用 X - β(Y - E[Y]) 去估计，β 取最优系数 Cov(X,Y)/Var(Y)。举个实际例子：模拟期权价格的时候，用标的资产本身做控制变量，因为它的期望是已知的。我在做工期模拟时也用过这招，用"所有任务都取均值时的总工期"作为控制变量，因为那个数字是确定的，结果方差降了七八成，同样的精度只要原来五分之一的样本。

分层抽样对一维问题效果极好。把均匀分布切成十层，每层里采一个点，等于强制让样本在 [0,1) 上均匀铺开，消除了聚集造成的波动。但维度一高就没法用了，因为层数会指数爆炸，这时候得换成拉丁超立方抽样，它保证每一维上的投影都是分层的，是个折中方案。

4.3 并行化与工程落地时的几个坑

样本量上到千万级，单进程就跑得慢了。蒙特卡洛天然适合并行，因为每次试验互相独立。但并行有两个坑。

第一个是随机数流的独立性。如果你在多个进程里都用同一个种子，那每个进程生成的序列完全一样，等于白跑。正确做法是用SeedSequence派生独立的子流：

from numpy.random import SeedSequence, default_rng ss = SeedSequence(20240607) child_seeds = ss.spawn(8) streams = [default_rng(s) for s in child_seeds]

这样八个进程拿到的随机序列在统计上互相独立，结果合并之后才是有效的。我见过有人用"主种子 + 进程号"的方式凑合，短样本下看着没事，样本一多就能测出不希望的相关性。

第二个坑是结果合并时忘了加权。如果不同进程跑的样本量不一样（比如按批次动态分配），合并的时候不能简单取平均，得按样本量加权。更隐蔽的是估计概率时的合并：正确做法是把所有进程的"命中次数"加总，再除以总样本数，而不是把各个进程的概率取平均。后者在样本量不均衡时会产生偏差。

第三个是内存。有人喜欢一次性生成一个 (10_000_000, 50) 的矩阵，那是五亿个浮点数，四个 G 起步，直接爆内存。正确的做法是分批生成、增量累加。写模拟脚本时养成"能增量就不要全存"的习惯，尤其是并行环境下每个进程都要吃一份内存。

5. 我踩过的坑与问题排查速查表

5.1 六类高频问题与定位方法

做了这么多年模拟，出问题的地方其实高度集中。整理成表，出问题时按顺序排查：

第一条我要多啰嗦两句。判断样本量够不够，有个很实用的自检方法：跑三组不同种子，看三个结果的离散程度。如果三组结果的差异和你算出来的标准误量级一致，说明没问题；如果差异远大于标准误，那肯定有别的问题，比如某处有隐藏的相关性。这个方法我在没有统计工具的环境下也常用，属于手边就能做的体检。

第三条也值得一提。"概率估计总是 0"这个现象，八成是因为你的阈值设在分布的支撑集之外了。比如你用三角分布模拟工期，最短 20 天，然后问"10 天内完成的概率"，那必然是 0，不是 bug。但如果阈值明明在合理范围内，概率却是 0，那就得回头检查单位——天数还是小时，元还是万元，这种错误在多个数据源合并时极其常见。

5.2 几条没人写在文档里的实操心得

第一条：先跑小样本验证逻辑，再上大样本。我的固定流程是先用 1000 次跑通，打印中间变量，确认每一步的分布形状、量纲、边界都合理，再切到十万次以上。直接上大样本的坏处是，跑十分钟出来一个错的结果，你还得再花十分钟验证它错在哪。小样本跑一遍只要几毫秒，改起来毫无心理负担。

第二条：永远把你的分布假设写进输出。我在交付任何模拟结果时，都会附上一段"本次模拟基于以下假设"，把每个输入分布的类型和参数列清楚。这不是形式主义。同一个问题，换一组假设，结论可能完全反过来。把假设亮出来，别人才知道该质疑哪里，你也能在事后追溯到底哪一步的假设出了问题。我自己就有过一次教训，某个参数的分布选错了，但因为当初没记录，排查花了整整两天。

第三条：对极端值保持警惕，但不要急着删。模拟出来的样本里偶尔会出现离谱的值，比如工期模拟里冒出个 3000 天。第一反应不该是截断，而是问一句：这个值出现的概率有多大，如果真实世界真发生了会怎样。如果用np.clip截断，等于人为制造了一个概率堆积点，会扭曲尾部分位数。我现在的做法是保留原始样本，在统计时单独报告"超过某个阈值的样本占比"，让决策者自己判断要不要关注。真要截断，也得明确写出来，不能悄悄处理。

第四条：结果的可视化比数字更有说服力。一张直方图配上几条分位数竖线，比一堆表格更容易让人理解"风险在哪里"。我做汇报时基本都用这个组合：左边直方图，右边累积分布曲线，关键分位数用竖线标出来。累积分布曲线特别好用，因为它能直接读出"低于某个值的概率"，而这正是绝大多数决策者真正想知道的东西。相比一张满是均值、标准差、偏度的表格，这条曲线的沟通效率高出一个量级。

第五条：把模拟脚本当成可复现的实验来管理。种子固定、参数外置、输出带时间戳和参数快照。这听起来像工程规范，但实际收益很直接——当你两周后需要回答"当时那版模型的假设是什么"的时候，翻一下输出文件就够了，不用去考古自己的代码。我现在的每个模拟脚本都会自动把参数 dump 成一份 JSON 存在结果旁边，这个习惯养成之后，返工率明显下降。

最后再分享一个我常用的组合拳。面对一个原来没做过模拟的新问题，我的套路是：第一步，找一个有解析解的简化版本，用模拟跑一遍，确认两者对得上，这验证了框架；第二步，逐步把简化的假设放开，每放开一条就记一次结果变化，这让你知道哪个假设最关键；第三步，在全复杂版本上做敏感度分析，找出影响最大的输入。这三步走完，结论基本就站得住了，而且你对每个结论的来龙去脉都心里有数。相比一上来就堆一个几百行的大模型，这套流程慢一点，但省下来的返工时间远不止那点。