14. 故障树不确定性#

除了Anaconda中已有的库外，本讲座还需要以下库：

!pip install quantecon tabulate

Hide code cell output

 
Requirement already satisfied: quantecon in /home/runner/miniconda3/envs/quantecon/lib/python3.13/site-packages (0.13.1)
Requirement already satisfied: tabulate in /home/runner/miniconda3/envs/quantecon/lib/python3.13/site-packages (0.10.0)
Requirement already satisfied: numba>=0.56.0 in /home/runner/miniconda3/envs/quantecon/lib/python3.13/site-packages (from quantecon) (0.65.1)
Requirement already satisfied: numpy>=1.17.0 in /home/runner/miniconda3/envs/quantecon/lib/python3.13/site-packages (from quantecon) (2.4.6)
Requirement already satisfied: scipy>=1.5.0 in /home/runner/miniconda3/envs/quantecon/lib/python3.13/site-packages (from quantecon) (1.18.0)
Requirement already satisfied: llvmlite<0.48,>=0.47.0dev0 in /home/runner/miniconda3/envs/quantecon/lib/python3.13/site-packages (from numba>=0.56.0->quantecon) (0.47.0)

14.1. 概述#

本讲将运用基本工具来近似计算由多个关键部件组成的系统的年度故障率的概率分布。

我们将使用对数正态分布来近似关键部件的概率分布。

为了近似描述系统总故障率（表示为 \(n\) 个对数正态随机变量之和）的概率分布，我们计算这些分布的卷积。

我们将使用以下概念和工具：

  • 对数正态分布

  • 描述独立随机变量之和的概率分布的卷积定理

  • 用于近似多组件系统故障率的故障树分析

  • 用于描述不确定概率的层次概率模型

  • 傅里叶变换和傅里叶逆变换作为计算序列卷积的高效方法

参见

关于傅里叶变换的更多信息，请参见 循环矩阵 以及 协方差平稳过程 和 谱估计。

El-Shanawany et al. [2018] 和 Greenfield and Sargent [1993] 应用了这些方法来近似核设施安全系统的故障概率。

这些技术响应了 Apostolakis [1990] 提出的关于量化安全系统可靠性不确定性的建议。

本讲座将使用以下导入和设置：

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
FONTPATH = "fonts/SourceHanSerifSC-SemiBold.otf"
mpl.font_manager.fontManager.addfont(FONTPATH)
plt.rcParams['font.family'] = ['Source Han Serif SC']

from scipy.signal import fftconvolve
from tabulate import tabulate
import quantecon as qe

14.2. 对数正态分布#

如果随机变量 \(x\) 服从均值为 \(\mu\)、方差为 \(\sigma^2\) 的正态分布，那么 \(y = \exp(x)\) 服从参数为 \(\mu, \sigma^2\) 的对数正态分布。

备注

我们将 \(\mu\) 和 \(\sigma^2\) 称为参数而不是均值和方差，因为：

  • \(\mu\) 和 \(\sigma^2\) 是 \(x = \log(y)\) 的均值和方差

  • 它们不是 \(y\) 的均值和方差

  • \(y\) 的均值是 \(\exp(\mu + \frac{1}{2}\sigma^2)\)，方差是 \((e^{\sigma^2} - 1) e^{2\mu + \sigma^2}\)

对数正态随机变量 \(y\) 始终是非负的。

\(y\) 的概率密度函数是

(14.1)#\[f(y) = \frac{1}{y \sigma \sqrt{2 \pi}} \exp \left( \frac{- (\log y - \mu)^2 }{2 \sigma^2} \right), \quad y \geq 0\]

对数正态随机变量的重要特性是：

(14.2)#\[\begin{split}\begin{aligned} \text{均值:} & \quad e ^{\mu + \frac{1}{2} \sigma^2} \\ \text{方差:} & \quad (e^{\sigma^2} - 1) e^{2 \mu + \sigma^2} \\ \text{中位数:} & \quad e^\mu \\ \text{众数:} & \quad e^{\mu - \sigma^2} \\ \text{0.95 分位数:} & \quad e^{\mu + 1.645 \sigma} \\ \text{0.95/0.05 分位数比:} & \quad e^{3.29 \sigma} \end{aligned}\end{split}\]

14.2.1. 稳定性性质#

回顾独立正态分布随机变量具有以下稳定性性质：

如果 \(x_1 \sim N(\mu_1, \sigma_1^2)\) 和 \(x_2 \sim N(\mu_2, \sigma_2^2)\) 是独立的，那么 \(x_1 + x_2 \sim N(\mu_1 + \mu_2, \sigma_1^2 + \sigma_2^2)\)。

独立的对数正态分布具有不同的稳定性性质：独立对数正态随机变量的乘积也是对数正态分布。

具体来说，如果 \(y_1\) 是参数为 \((\mu_1, \sigma_1^2)\) 的对数正态分布，且 \(y_2\) 是参数为 \((\mu_2, \sigma_2^2)\) 的对数正态分布，那么 \(y_1 y_2\) 是参数为 \((\mu_1 + \mu_2, \sigma_1^2 + \sigma_2^2)\) 的对数正态分布。

警告

虽然两个对数正态分布的乘积是对数正态分布，但两个对数正态分布的和却不是对数正态分布。

这个观察结果引出了本讲座的核心挑战：近似独立对数正态随机变量之和的概率分布。

14.3. 卷积定理#

设 \(x\) 和 \(y\) 是概率密度分别为 \(f(x)\) 和 \(g(y)\) 的独立随机变量，其中 \(x, y \in \mathbb{R}\)。

设 \(z = x + y\)。

那么 \(z\) 的概率密度为

(14.3)#\[h(z) = (f * g)(z) \equiv \int_{-\infty}^\infty f(\tau) g(z - \tau) d\tau\]

其中 \((f*g)\) 表示 \(f\) 和 \(g\) 的卷积。

对于非负随机变量，这可以特化为

(14.4)#\[h(z) = (f * g)(z) \equiv \int_{0}^z f(\tau) g(z - \tau) d\tau\]

14.3.1. 离散卷积#

我们将使用卷积公式的离散化版本。

我们将 \(f\) 和 \(g\) 都替换为离散化的对应形式，并归一化使其和为 1。

离散卷积公式为

(14.5)#\[h_n = (f*g)_n = \sum_{m=0}^n f_m g_{n-m}, \quad n \geq 0\]

这计算了两个离散随机变量之和的概率质量函数。

14.3.2. 示例：离散分布#

考虑两个概率质量函数：

\[ f_j = \mathbb{P}\{X = j\}, \quad j = 0, 1 \]

和

\[ g_j = \mathbb{P}\{Y = j\}, \quad j = 0, 1, 2, 3 \]

\(Z = X + Y\) 的分布由卷积 \(h = f * g\) 给出。

# 定义概率质量函数
f = [0.75, 0.25]
g = [0.0, 0.6, 0.0, 0.4]

# 使用两种方法计算卷积
h = np.convolve(f, g)
hf = fftconvolve(f, g)

print(f"f = {f}, sum = {np.sum(f):.3f}")
print(f"g = {g}, sum = {np.sum(g):.3f}")
print(f"h = {h}, sum = {np.sum(h):.3f}")
print(f"hf = {hf}, sum = {np.sum(hf):.3f}")

f = [0.75, 0.25], sum = 1.000
g = [0.0, 0.6, 0.0, 0.4], sum = 1.000
h = [0.   0.45 0.15 0.3  0.1 ], sum = 1.000
hf = [0.   0.45 0.15 0.3  0.1 ], sum = 1.000

numpy.convolve 和 scipy.signal.fftconvolve 都得到相同的结果，但对于长序列，fftconvolve 要快得多。

为了提高效率，本讲座将始终使用 fftconvolve。

14.4. 近似连续分布#

现在我们验证离散化分布能否准确近似来自底层连续分布的样本。

我们从三个独立的对数正态随机变量中生成25,000个样本，并计算它们的两两之和与三者之和。

然后我们将样本的直方图与离散化分布的直方图进行比较。

# 设置对数正态分布的参数
μ, σ = 5.0, 1.0
n_samples = 25000

# 生成样本
rng = np.random.default_rng(1234)
s1 = rng.lognormal(μ, σ, n_samples)
s2 = rng.lognormal(μ, σ, n_samples)
s3 = rng.lognormal(μ, σ, n_samples)

# 计算和
ssum2 = s1 + s2
ssum3 = s1 + s2 + s3

# 绘制 s1 的直方图
fig, ax = plt.subplots()
ax.hist(s1, 1000, density=True, alpha=0.6)
ax.set_xlabel('数值')
ax.set_ylabel('密度')
plt.show()

findfont: Failed to find font weight normal, now using 600.
_images/dee0f367991220e202c1995bae17de3dda8e44799af4cbb376d7142078703aed.png

图 14.1 单个对数正态分布的样本直方图#

# 绘制两个对数正态分布之和的直方图
fig, ax = plt.subplots()
ax.hist(ssum2, 1000, density=True, alpha=0.6)
ax.set_xlabel('数值')
ax.set_ylabel('密度')
plt.show()
_images/9d694bbf6b28dc334b2b1a26bc8a28e4c5dff8938dd8959e143c210ff117eb61.png

图 14.2 两个对数正态分布之和的直方图#

# 绘制三个对数正态分布之和的直方图
fig, ax = plt.subplots()
ax.hist(ssum3, 1000, density=True, alpha=0.6)
ax.set_xlabel('数值')
ax.set_ylabel('密度')
plt.show()
_images/cb8dcb503d1754c863cb2c084903b38e70be150662923dfd61b77f5ff4444b9c.png

图 14.3 三个对数正态分布之和的直方图#

让我们验证样本均值是否与理论均值相匹配：

samp_mean = np.mean(s2)
theoretical_mean = np.exp(μ + σ**2 / 2)

print(f"理论均值: {theoretical_mean:.3f}")
print(f"样本均值: {samp_mean:.3f}")

理论均值: 244.692
样本均值: 243.197

14.5. 离散化对数正态分布#

我们定义辅助函数来创建对数正态概率密度函数的离散化版本。

我们手动写出该密度函数，以便与公式 (14.1) 对照；scipy.stats.lognorm(s=σ, scale=np.exp(μ)).pdf(x) 计算的是同样的内容。

def lognormal_pdf(x, μ, σ):
    """
    计算对数正态概率密度函数。
    """
    p = 1 / (σ * x * np.sqrt(2 * np.pi)) \
            * np.exp(-0.5 * ((np.log(x) - μ) / σ)**2)
    return p


def discretize_lognormal(μ, σ, I, m):
    """
    在网格 0, m, 2m, ..., 直到 I 上离散化对数正态分布。

    参数
    ----------
    μ, σ : 对数正态分布的参数
    I    : 网格的上限，用于截断右尾
    m    : 网格点之间的间距，决定分辨率

    返回
    -------
    p_array      : 在网格上计算得到的密度
    p_array_norm : 隐含的概率质量函数，其和为1
    x            : 网格本身，共有 I / m 个点
    """
    x = np.arange(1e-7, I, m)
    p_array = lognormal_pdf(x, μ, σ)
    p_array_norm = p_array / np.sum(p_array)
    return p_array, p_array_norm, x

有两个独立的选择决定了这一近似的质量，最好将它们分清楚。

  • \(I\) 决定网格的终止点，因此它控制我们舍弃了多少右尾部分

  • \(m\) 决定网格点之间的间距，因此它控制分辨率

网格共有 \(I/m\) 个点，因此在 \(m\) 固定时增大 \(I\) 会扩大覆盖范围，而在 \(I\) 固定时减小 \(m\) 会提高精度。

一旦 \(I\) 足够大，使得几乎没有概率质量落在其之外，再进一步增大它就不会改变结果，此时只有 \(m\) 起作用。

练习 14.1 要求你验证这一点。

备注

scipy.signal.fftconvolve 会在内部自动将输入填充到合适的长度，因此无需将 \(I/m\) 选为 2 的幂。

# 设置网格参数
p = 15
I = 2**p  # 网格的终止点：截断右尾
m = 0.1   # 网格点之间的间距：决定分辨率

让我们直观地看一下离散化分布对连续对数正态分布的近似效果：

# 计算离散化的概率密度函数
pdf, pdf_norm, x = discretize_lognormal(μ, σ, I, m)

# 绘制离散化的概率密度函数与直方图的对比
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(x, pdf, 'r-', lw=2, label='离散化概率密度函数')
ax.hist(s1, 1000, density=True, alpha=0.6, label='样本直方图')
ax.set_xlim(0, 2500)
ax.set_xlabel('数值')
ax.set_ylabel('密度')
ax.legend()
plt.show()
_images/039c1b0a24fdf6a264566c8842b037fbad999411883e39a4be3f72e780f3b898.png

图 14.4 离散化密度与样本的对比#

现在让我们验证离散化分布是否具有正确的均值：

# 从离散化的概率密度函数计算均值
mean_discrete = np.sum(x * pdf_norm)
mean_theory = np.exp(μ + 0.5 * σ**2)

print(f"理论均值: {mean_theory:.3f}")
print(f"离散化均值: {mean_discrete:.3f}")

理论均值: 244.692
离散化均值: 244.691

14.6. 概率质量函数的卷积#

现在我们使用卷积定理来计算上面参数化的两个对数正态随机变量之和的概率分布。

我们还将计算上面构造的三个对数正态分布之和的概率分布。

对于长序列，scipy.signal.fftconvolve 比 numpy.convolve 快得多，因为它使用了快速傅里叶变换。

让我们先定义傅里叶变换和傅里叶逆变换

14.6.1. 快速傅里叶变换#

序列 \(\{x_t\}_{t=0}^{T-1}\) 的傅里叶变换是

(14.6)#\[x(\omega_j) = \sum_{t=0}^{T-1} x_t \exp(-i \omega_j t)\]

其中 \(\omega_j = \frac{2\pi j}{T}\)，\(j = 0, 1, \ldots, T-1\)。

序列 \(\{x(\omega_j)\}_{j=0}^{T-1}\) 的傅里叶逆变换是

(14.7)#\[x_t = T^{-1} \sum_{j=0}^{T-1} x(\omega_j) \exp(i \omega_j t)\]

序列 \(\{x_t\}_{t=0}^{T-1}\) 和 \(\{x(\omega_j)\}_{j=0}^{T-1}\) 包含相同的信息。

方程对 (14.6) 和 (14.7) 说明了如何从一个序列恢复其傅里叶对应序列。

程序 scipy.signal.fftconvolve 利用了两个序列 \(\{f_k\}\)、\(\{g_k\}\) 的卷积可以通过以下方式计算的定理：

  • 计算序列 \(\{f_k\}\) 和 \(\{g_k\}\) 的傅里叶变换 \(F(\omega)\)、\(G(\omega)\)

  • 形成乘积 \(H (\omega) = F(\omega) G (\omega)\)

  • 卷积 \(f * g\) 是 \(H(\omega)\) 的傅里叶逆变换

快速傅里叶变换和相关的快速傅里叶逆变换能够非常快速地执行这些计算。

这就是 fftconvolve 使用的算法。

让我们做一个预热计算，比较 numpy.convolve 和 scipy.signal.fftconvolve 所需的时间

我们的三个分量是同分布的，因此单次离散化就足以适用于所有分量。

# 离散化对数正态分布；三个分量是独立同分布的
_, pmf1, x = discretize_lognormal(μ, σ, I, m)
pmf2 = pmf3 = pmf1

# 直接卷积的成本为 O(N²)，因此我们只对一个较短的前缀部分计时
short = pmf1[:20_000]

with qe.Timer() as timer_numpy:
    np.convolve(short, short)
time_numpy = timer_numpy.elapsed

with qe.Timer() as timer_fft:
    fftconvolve(short, short)
time_fft = timer_fft.elapsed

print(f"在 {len(short):,} 个点上：")
print(f"  np.convolve: {time_numpy:.4f} 秒")
print(f"  fftconvolve: {time_fft:.4f} 秒")
print(f"  加速倍数:     {time_numpy / time_fft:.0f}x")

0.0424 seconds elapsed
0.0016 seconds elapsed
在 20,000 个点上：
  np.convolve: 0.0424 秒
  fftconvolve: 0.0016 秒
  加速倍数:     26x

随着序列长度的增加，这一差距会迅速扩大，因为直接卷积的计算成本为 \(O(N^2)\)，而 FFT 方法的计算成本为 \(O(N \log N)\)。

在下面使用的完整网格上，直接方法的速度会更慢。

# 使用快速方法完成完整计算
conv_fft = fftconvolve(fftconvolve(pmf1, pmf2), pmf3)
print(f"每个分量的网格点数: {len(pmf1):,}")

每个分量的网格点数: 327,680

现在让我们将计算得到的两个对数正态随机变量之和的概率质量函数近似值与我们上面形成的样本直方图进行对比绘制

# 计算两个分布的卷积以进行比较
conv2 = fftconvolve(pmf1, pmf2)

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(x, conv2[:len(x)] / m, 'r-', lw=2, label='卷积 (FFT)')
ax.hist(ssum2, 1000, density=True, alpha=0.6, label='样本直方图')
ax.set_xlim(0, 5000)
ax.set_xlabel('数值')
ax.set_ylabel('密度')
ax.legend()
plt.show()
_images/9b31f532aee4c6687fb16801d73ac947aa2ae2f35807e230c2c5b1f2d9a57d98.png

图 14.5 卷积与样本的对比，两个分量#

现在我们展示三个对数正态随机变量之和的图：

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(x, conv_fft[:len(x)] / m, 'r-', lw=2, label='卷积 (FFT)')
ax.hist(ssum3, 1000, density=True, alpha=0.6, label='样本直方图')
ax.set_xlim(0, 5000)
ax.set_xlabel('数值')
ax.set_ylabel('密度')
ax.legend()
plt.show()
_images/d09b1740605dabf080a5b9d29107b89e2f2a009c09a7dd131c2f1e36e6ba9395.png

图 14.6 卷积与样本的对比，三个分量#

让我们验证均值是否正确

# 两个分布之和的均值
mean_conv2 = np.sum(x * conv2[:len(x)])
mean_theory2 = 2 * np.exp(μ + 0.5 * σ**2)

print(f"两个分布之和:")
print(f"  理论均值: {mean_theory2:.3f}")
print(f"  计算均值: {mean_conv2:.3f}")

两个分布之和:
  理论均值: 489.384
  计算均值: 489.381

# 三个分布之和的均值
mean_conv3 = np.sum(x * conv_fft[:len(x)])
mean_theory3 = 3 * np.exp(μ + 0.5 * σ**2)

print(f"三个分布之和:")
print(f"  理论均值: {mean_theory3:.3f}")
print(f"  计算均值: {mean_conv3:.3f}")

三个分布之和:
  理论均值: 734.076
  计算均值: 734.071

14.7. 故障树分析#

我们即将应用卷积定理来计算故障树分析中顶事件的概率。

在应用卷积定理之前，我们首先描述将组成事件与我们要量化其故障率的顶事件连接起来的模型。

正如 El-Shanawany et al. [2018] 所描述的，故障树分析是一种广泛使用的评估系统可靠性的技术。

为了构建统计模型，我们反复使用所谓的稀有事件近似。

14.7.1. 稀有事件近似#

我们想要计算事件 \(A \cup B\) 的概率。

对于事件 \(A\) 和 \(B\)，并集的概率为

\[ P(A \cup B) = P(A) + P(B) - P(A \cap B) \]

其中 \(A \cup B\) 是事件 \(A\) 或 \(B\) 发生的情况，\(A \cap B\) 是事件 \(A\) 和 \(B\) 都发生的情况。

如果 \(A\) 和 \(B\) 是独立的，那么 \(P(A \cap B) = P(A) P(B)\)。

当 \(P(A)\) 和 \(P(B)\) 都很小时，\(P(A) P(B)\) 就更小。

稀有事件近似为

\[ P(A \cup B) \approx P(A) + P(B) \]

这种近似方法在系统故障分析中被广泛使用。

14.7.2. 系统故障概率#

考虑一个具有 \(n\) 个关键组件的系统，当任何一个组件发生故障时，系统就会发生故障。

我们假设：

  • 每个组件 \(A_i\) 的故障概率 \(P(A_i)\) 都很小

  • 组件故障在统计上是独立的

我们反复应用稀有事件近似，得到系统故障概率问题的以下公式：

\[ P(F) \approx P(A_1) + P (A_2) + \cdots + P (A_n) \]

或

(14.8)#\[P(F) \approx \sum_{i=1}^n P(A_i)\]

其中 \(P(F)\) 是系统故障概率。

每个事件的概率以每年故障率的形式记录。

备注

严格来说，每年的故障率与一年内的故障概率是不同的概念。

对于稀有事件，两者几乎相等，因为当 \(\lambda\) 很小时，\(1 - e^{-\lambda} \approx \lambda\)。

正是同样让我们能够跨组件相加概率的近似方法，也让我们能够在故障率与概率之间相互转换，因此我们遵循可靠性文献的惯例，在此交替使用这两个术语。

14.8. 未知的故障率#

现在我们来讨论真正感兴趣的问题，遵循 El-Shanawany et al. [2018] 和 Greenfield and Sargent [1993] 的方法，秉承 Apostolakis [1990] 的精神。

组件故障率 \(P(A_i)\) 并非精确已知，需要进行估计。

我们通过指定概率的概率来解决这个问题，这体现了不了解作为故障树分析输入的构成概率的一种概念。

因此，我们假设系统分析师对系统组件的故障率 \(P(A_i), i =1, \ldots, n\) 存在不确定性。

分析师通过将系统的故障概率 \(P(F)\) 和每个组件概率 \(P(A_i)\) 视为随机变量来应对这种情况。

  • \(P(A_i)\) 概率分布的离散程度表征了分析师对故障概率 \(P(A_i)\) 的不确定性

  • \(P(F)\) 的隐含概率分布的离散程度表征了他对系统故障概率的不确定性

这就是所谓的层次化模型，其中分析师对概率 \(P(A_i)\) 本身也有概率估计。

备注

该模型中出现了两种截然不同的随机性，值得加以区分。

偶然性（Aleatory）不确定性是指某个组件在给定年份内是否发生故障的随机性，它由故障率 \(P(A_i)\) 来描述。

认知性（Epistemic）不确定性是指分析师对该故障率取值的无知，它由分析师赋予 \(P(A_i)\) 的对数正态分布来描述。

我们下面计算的分布是一个认知性对象：它描述的是分析师对某个故障率的了解程度，而不是系统实际发生故障的频率。

将两者区分开来是 Apostolakis [1990] 的核心建议。

分析师通过以下假设来形式化他的不确定性：

  • 故障概率 \(P(A_i)\) 本身是一个对数正态随机变量，其参数为 \((\mu_i, \sigma_i)\)。

  • 对于所有 \(i \neq j\) 的配对，故障率 \(P(A_i)\) 和 \(P(A_j)\) 在统计上是相互独立的。

分析师通过阅读工程论文中的可靠性研究来校准故障事件 \(i = 1, \ldots, n\) 的参数 \((\mu_i, \sigma_i)\)，这些研究考察了与所研究系统中使用的组件尽可能相似的组件的历史故障率。

分析师假设，这些关于年度故障率或故障时间的观测分散性的信息，可以帮助他预测零件在其系统中的性能表现。

分析师假设随机变量 \(P(A_i)\) 在统计上是相互独立的。

警告

独立性是一个很强的假设，也是可靠性分析师最为担心的一点。

设计缺陷、共用的电源、共同的维护团队，或单一的环境冲击，都可能同时使多个组件趋向故障。

这类共因故障会使 \(P(F)\) 分布的尾部远比独立性假设下计算出的结果更为肥厚，而这恰恰是安全监管者最为关心的区域。

练习 14.5 对这种差异究竟有多大进行了量化。

分析师想要近似系统的故障概率 \(P(F)\) 的概率质量函数和累积分布函数。

  • 我们说概率质量函数是因为我们对每个随机变量进行了离散化，正如前文描述的那样。

分析师通过重复应用卷积定理来计算顶事件 \(F\)（即系统故障）的概率质量函数，以计算独立对数正态随机变量之和的概率分布，如方程 (14.8) 所述。

14.9. 应用：废物提升机失效率#

现在我们分析一个具有 \(n = 14\) 个组件的真实案例。

该应用估计了核废料设施中一个关键提升机的年度故障率。

监管机构要求系统的设计能够使顶事件的故障率以高概率保持在较小值。

14.9.1. 模型设定#

这个例子是 Greenfield and Sargent [1993] 第27页表10中描述的设计方案B-2（案例I）。

该表描述了十四个对数正态随机变量的参数 \(\mu_i, \sigma_i\)，这些随机变量由七对独立同分布的随机变量组成。

  • 在每一对内，参数 \(\mu_i, \sigma_i\) 是相同的

  • 如 Greenfield and Sargent [1993] 第27页表10所述，七个唯一概率 \(P(A_i)\) 的对数正态分布参数已被校准为以下Python代码中的值：

# 组件故障率参数
# (参见 Greenfield & Sargent 1993 表10)
params = [
    (4.28, 1.1947),   # 组件类型 1
    (3.39, 1.1947),   # 组件类型 2
    (2.795, 1.1947),  # 组件类型 3
    (2.717, 1.1947),  # 组件类型 4
    (2.717, 1.1947),  # 组件类型 5
    (1.444, 1.4632),  # 组件类型 6
    (-0.040, 1.4632), # 组件类型 7 (出现8次)
]

备注

由于故障率都很小，这些对数正态分布实际上描述的是 \(P(A_i) \times 10^{-9}\)。

所以我们将在概率质量函数和相关累积分布函数的 \(x\) 轴上标注的概率应该乘以 \(10^{-09}\)

我们定义一个辅助函数来查找数组索引：

def find_nearest(array, value):
    """
    数组中最接近给定值的元素的索引。

    应用于累积分布函数时，这会返回累积概率最接近目标值的网格点，
    对于足够精细离散化的分布而言，这与将分位数定义为满足
    CDF(x) >= q 的最小 x 值是没有区别的。
    """
    array = np.asarray(array)
    idx = (np.abs(array - value)).argmin()
    return idx

我们在以下代码中计算所需的十三个卷积。

(请随意尝试不同的幂参数 \(p\) 值，我们用它来设置网格中的点数，以构建离散化连续对数正态分布的概率质量函数。)

# 设置网格参数
p = 15
I = 2**p
m = 0.05

# 离散化所有组件的故障率分布
# 前6个组件使用各自独特的参数，后8个共享相同的参数
component_pmfs = []
for μ, σ in params[:6]:
    _, pmf, x = discretize_lognormal(μ, σ, I, m)
    component_pmfs.append(pmf)

# 添加8份组件类型7的副本
μ7, σ7 = params[6]
_, pmf7, x = discretize_lognormal(μ7, σ7, I, m)
component_pmfs.extend([pmf7] * 8)

# 通过依次卷积计算系统故障分布
with qe.Timer() as timer:
    system_pmf = component_pmfs[0]
    for pmf in component_pmfs[1:]:
        system_pmf = fftconvolve(system_pmf, pmf)

print(f"13次卷积所需时间: {timer.elapsed:.4f} 秒")

# 卷积结果保持相同的网格间距，但延伸得更远
system_grid = np.arange(len(system_pmf)) * m
print(f"结果中的网格点数: {len(system_pmf):,}")

2.2031 seconds elapsed
13次卷积所需时间: 2.2031 秒
结果中的网格点数: 9,175,027

在绘制累积分布函数之前，我们先来看看密度本身。

fig, ax = plt.subplots(figsize=(10, 6))
upper = 2000
ax.plot(system_grid[:int(upper/m)], system_pmf[:int(upper/m)] / m, 'b-', lw=2)
ax.set_xlabel(r'故障率 (每年 $\times 10^{-9}$)')
ax.set_ylabel('密度')
plt.show()

findfont: Failed to find font weight normal, now using 600.
_images/42197ce6fe4103684f2cebf3eebbf20e4979588f434e5df93beb5c2d9e48134b.png

图 14.7 系统故障率的密度#

该密度明显右偏：一条长长的上尾远远延伸至分布主体之外。

正是这种不对称性使得单一的故障率点估计成为一个糟糕的概括，这也是分析者转而报告分位数的原因。

现在我们绘制一个与 Greenfield and Sargent [1993] 第29页图5中的累积分布函数(CDF)相对应的图

# 计算累积分布函数
cdf = np.cumsum(system_pmf)

# 绘制累积分布函数
Nx = 1400
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(x[:int(Nx / m)], cdf[:int(Nx / m)], 'b-', lw=2)

# 添加关键分位数的参考线
quantile_levels = [0.05, 0.10, 0.50, 0.90, 0.95]
for q in quantile_levels:
    ax.axhline(q, color='gray', linestyle='--', alpha=0.5)

ax.set_xlim(0, Nx)
ax.set_ylim(0, 1)
ax.set_xlabel(r'故障率 (每年 $\times 10^{-9}$)')
ax.set_ylabel('累积概率')
plt.show()
_images/adc5ec80a4842054829ed9bebd4dee50f1f3afd8de1d2454f2f3c18c6afef993.png

图 14.8 系统故障率的累积分布函数#

我们还展示一个与 Greenfield and Sargent [1993] 第28页表11相对应的表，列出了系统故障率分布的关键分位数

# Greenfield 和 Sargent (1993) 表11中报告的百分位数，
# 及其发表的数值，单位为每年 10^-9
reference = {1.0: 77, 10.0: 130, 50.0: 263, 66.5: 341,
             85.0: 513, 95.0: 811, 99.0: 1480, 99.78: 2490}

table_data = []
for pc, published in reference.items():
    ours = system_grid[find_nearest(cdf, pc/100)]
    table_data.append([f"{pc}%", f"{ours:.1f}", published,
                       f"{100*(ours - published)/published:+.1f}%"])

print("\n系统故障率分位数 (×10^-9 每年):")
print(tabulate(table_data,
      headers=['百分位数', '本文计算值', 'Greenfield-Sargent', '差异'],
      tablefmt='grid'))

系统故障率分位数 (×10^-9 每年):
+------------+--------------+----------------------+--------+
| 百分位数   |   本文计算值 |   Greenfield-Sargent | 差异   |
+============+==============+======================+========+
| 1.0%       |         76.2 |                   77 | -1.1%  |
+------------+--------------+----------------------+--------+
| 10.0%      |        128.2 |                  130 | -1.4%  |
+------------+--------------+----------------------+--------+
| 50.0%      |        260.6 |                  263 | -0.9%  |
+------------+--------------+----------------------+--------+
| 66.5%      |        338.6 |                  341 | -0.7%  |
+------------+--------------+----------------------+--------+
| 85.0%      |        509.4 |                  513 | -0.7%  |
+------------+--------------+----------------------+--------+
| 95.0%      |        807.6 |                  811 | -0.4%  |
+------------+--------------+----------------------+--------+
| 99.0%      |       1470.2 |                 1480 | -0.7%  |
+------------+--------------+----------------------+--------+
| 99.78%     |       2474.9 |                 2490 | -0.6%  |
+------------+--------------+----------------------+--------+

我们计算出的分位数与已发表的数值相差在百分之一点五以内，且都略微偏低。

这些微小的差异反映了所报告参数 \(\mu_i, \sigma_i\) 的精度、网格间距 \(m\) 以及网格截断点的影响。

14.9.2. 解读结果#

这张表中的数字，而非其中任何单一的数值，才是分析的产出。

中位故障率约为每年 \(261 \times 10^{-9}\)，而第95百分位数约为 \(808 \times 10^{-9}\)，是中位数的三倍。

这种差异并不是在陈述提升机故障的实际频率；它陈述的是分析者对提升机故障频率所知有多么有限。

再来看看分布的均值落在何处。

mean_rate = np.sum(system_grid * system_pmf)
mean_percentile = 100 * cdf[find_nearest(system_grid, mean_rate)]

print(f"平均故障率: {mean_rate:.1f} × 10⁻⁹ 每年")
print(f"均值位于第 {mean_percentile:.1f} 百分位")

平均故障率: 338.1 × 10⁻⁹ 每年
均值位于第 66.4 百分位

由于该分布是偏斜的，均值远高于中位数，大约位于第66百分位处。

这就是为什么 Greenfield and Sargent [1993] 的表11将均值与第66.5百分位数一并记录的原因。

这一实际意义正是最初激发这项研究的动机。

# 美国能源部1990年风险评估中使用的点估计值，以相同单位表示
doe_estimate = 220    # 每年 2.2 × 10^-7

pct = 100 * cdf[find_nearest(system_grid, doe_estimate)]
print(f"美国能源部的点估计值 {doe_estimate} × 10⁻⁹ 位于"
      f"第 {pct:.0f} 百分位")
print(f"因此分析者认为真实故障率超过该值的概率为 {100-pct:.0f}%")

美国能源部的点估计值 220 × 10⁻⁹ 位于第 39 百分位
因此分析者认为真实故障率超过该值的概率为 61%

若一项分析仅报告单一数值来代替整个分布，则无法传达上述任何信息。

Greenfield and Sargent [1993] 正是提出了这一点：根据他们的图，他们将能源部的点估计值定位在第36百分位，并得出结论认为真实故障率高于该值的概率大约为64%。

14.10. 练习#

练习 14.1

我们的离散化涉及两个独立的选择：网格在哪里截断，\(I = 2^p\)，以及网格间距有多细，\(m\)。

研究这两者各自控制什么。

  1. 固定 \(m = 0.05\)，对 \(p = 10, 11, \ldots, 15\) 计算系统故障率的中位数、第95百分位数和第99.78百分位数。对每个 \(p\)，还要利用 \(\sum_i \mathbb{P}\{P(A_i) > I\}\) 计算截断所丢弃的概率质量。

  2. 固定 \(p = 14\)，对 \(m = 0.4, 0.2, 0.1, 0.05, 0.025\) 重复上述计算。

  3. 哪个统计量对哪个选择敏感，为什么？讲座中使用的 \(p = 15\)，\(m = 0.05\) 选择合理吗？

练习 14.2

稀有事件近似用 \(P(A) + P(B)\) 替代 \(P(A \cup B)\)，从而舍弃了 \(P(A \cap B)\)。

评估该近似在这里的效果如何。

  1. 以十四个组件各自的 平均 故障率作为代表值，比较 \(\sum_i p_i\) 与至少一个组件故障的精确概率 \(1 - \prod_i (1 - p_i)\)。

  2. 将所有故障率分别乘以 \(10^3\)、\(10^6\) 和 \(10^7\) 后重复上述计算，并报告每种情况下的相对误差。

  3. 在什么数量级上，该近似开始变得不可忽视？

练习 14.3

一位得知故障率第95百分位数过高的监管者会想知道应该改进哪些组件。

请通过计算在完全移除每种组件类型后系统故障率的第95百分位数，来回答这个问题（共七种组件类型）。

按对上尾部的贡献程度对各组件类型进行排序，并将该排序与各组件的平均故障率进行比较。

练习 14.4

除了通过卷积计算，我们也可以通过模拟来计算系统故障率的分布。

对全部十四个组件的故障率进行抽样、求和，并将得到的分位数与卷积方法的结果进行比较，样本量分别取 \(10^4\)、\(10^5\) 和 \(10^6\)。

比较中位数、第95百分位数和第99.78百分位数。

你会更倾向于使用哪种方法，为什么？

练习 14.5

整个计算假设十四个组件的故障率在统计上是相互独立的。

研究当它们不独立时会发生什么。

假设

\[ \log P(A_i) = \mu_i + \sigma_i \left( \sqrt{\rho}\, z_0 + \sqrt{1-\rho}\, z_i \right), \]

其中 \(z_0\) 是所有组件共同承受的冲击，\(z_1, \ldots, z_{14}\) 是各自独立的特质冲击，均服从标准正态分布。

每个组件仍然保持其原有的边缘分布，但任意两个组件在对数尺度上现在具有相关系数 \(\rho\)。

对 \(\rho = 0, 0.2, 0.5, 0.8\) 模拟系统故障率，并报告中位数、第95、第99和第99.9百分位数。

解释发生了什么，以及这对安全分析意味着什么。