2. 基于蒙特卡洛的输入不确定性传播

问题设定：一个含三个随机输入的悬臂梁端部挠度模型

本节选用一个可手算校核的物理模型：等截面悬臂梁在端部集中载荷下的最大挠度。其解析式为

δ = F·L³ / (3·E·I)，  I = b·h³ / 12

把 F（端部载荷）、E（弹性模量）、h（截面高度）视为随机输入，L、b 视为确定性常数。这样输出 δ 的不确定性完全由三个输入传播而来，便于观察蒙特卡洛估计的收敛行为。

注意单位一致性：计算时把 E 换算为 MPa（1 GPa = 1000 MPa），长度统一用 mm，则 δ 的单位为 mm。对数正态参数由中位数 m 与变异系数 CV 反推：σ_ln = sqrt(ln(1+CV²))，μ_ln = ln(m)。

蒙特卡洛传播的实现

核心只有三步：按各自分布独立抽样、逐样本代入模型、对输出样本做统计。下面给出可直接运行的 Python 实现，同时输出不同样本量下的统计量，用于观察收敛性。

import numpy as np

L, b = 1000.0, 20.0          # mm
mu_F, sd_F = 500.0, 50.0     # N
med_E, cv_E = 200.0, 0.08    # GPa
h_lo, h_hi = 18.0, 22.0      # mm

# 对数正态参数反推（E 以 GPa 抽样，再换算为 MPa）
sig_ln = np.sqrt(np.log(1 + cv_E**2))
mu_ln = np.log(med_E)

def deflection(F, E_GPa, h):
    E = E_GPa * 1000.0                 # GPa -> MPa (N/mm^2)
    I = b * h**3 / 12.0                # mm^4
    return F * L**3 / (3.0 * E * I)    # mm

def mc_run(n, seed=0):
    rng = np.random.default_rng(seed)
    F = rng.normal(mu_F, sd_F, n)
    E = rng.lognormal(mu_ln, sig_ln, n)
    h = rng.uniform(h_lo, h_hi, n)
    return deflection(F, E, h)

for n in [100, 1000, 10000, 100000]:
    d = mc_run(n, seed=42)
    print(f"n={n:>6d}  mean={d.mean():.4f}  std={d.std(ddof=1):.4f}  "
          f"P5={np.percentile(d,5):.4f}  P50={np.percentile(d,50):.4f}  "
          f"P95={np.percentile(d,95):.4f}")

固定随机种子 seed=42 是为了让不同样本量的结果可比较；实际工程中应改用独立种子并做多次重复，以评估估计量本身的波动。

收敛性观察：样本量如何影响估计

蒙特卡洛估计量的标准误随样本量按 1/√n 衰减。对均值而言，SE(mean) ≈ σ_δ/√n；对分位数而言，其标准误还与分布在该分位处的密度有关。因此样本量从 100 增到 10000（100 倍），均值标准误约缩小 10 倍，而分位数估计的改善通常略慢一些。

下面的图示把“输入分布 → 抽样 → 模型传播 → 输出统计 → 收敛判读”这条链与三个随机输入的分布形态放在一起，便于对照理解每个输入如何影响输出尾部。

结果判读与常见陷阱

• 均值与中位数不等：由于 δ 对 E、h 呈非线性（尤其 h³），输出分布右偏，均值通常大于中位数。报告时应同时给出 P50 与 P95，而不是只给均值。
• 样本量不足的典型表现：P95 在不同种子间跳动明显，而均值已相对稳定。这说明尾部估计比中心矩更“吃”样本。
• 单位错误是最常见的假结果来源：若忘记把 E 从 GPa 换成 MPa，δ 会整体偏小 1000 倍，且这种错误不会因增大样本量而消失。
• 独立性假设要写明：本节假设三个输入独立。若实际中 E 与 h 存在相关性，需改用联合分布或 Copula 抽样，否则方差估计会失真。

可复用的执行清单