强度仿真-主题092-不确定性量化的先进方法

kkchenjj

84人浏览 · 2026-03-11 21:16:44

kkchenjj · 2026-03-11 21:16:44 发布

主题092：不确定性量化的先进方法

摘要

不确定性量化是现代工程仿真的核心组成部分，对于确保结构可靠性和安全性具有重要意义。本教程系统阐述不确定性量化的先进理论与方法，包括谱方法、高斯过程回归、贝叶斯推断等前沿技术。通过多项式混沌展开、Karhunen-Loève展开、贝叶斯模型更新等典型案例，深入讲解高维不确定性问题的处理策略。教程提供完整的Python仿真代码，实现先进的不确定性量化分析，为复杂工程系统的可靠性评估和决策支持提供理论指导和技术工具。

关键词： 不确定性量化；多项式混沌；高斯过程；贝叶斯推断；可靠性分析

1. 引言

1.1 不确定性的来源与分类

工程实践中，结构分析和设计过程中存在多种不确定性来源：

偶然不确定性（Aleatoric Uncertainty）：

材料参数的分散性：同一批次材料的弹性模量、屈服强度存在自然波动
几何尺寸的制造误差：加工精度限制导致的尺寸偏差
载荷的随机性：风载、地震、交通载荷等具有随机特性
环境条件的波动：温度、湿度等环境因素的变化

这类不确定性是系统固有的，无法通过增加数据来消除，只能通过概率统计方法描述。

认知不确定性（Epistemic Uncertainty）：

模型形式误差：数学模型对物理现实的近似程度
参数识别误差：通过试验数据反演参数时的不确定性
边界条件的不确定性：实际工况与假设条件的差异
数据稀缺性：样本量不足导致的统计推断误差

这类不确定性源于知识缺乏，可以通过增加数据、改进模型来降低。

1.2 不确定性量化的重要性

不确定性量化在工程中的应用价值体现在：

可靠性评估：

计算结构的失效概率
确定安全系数的合理取值
评估设计余量的充分性

设计优化：

鲁棒设计：使性能对不确定性不敏感
可靠性优化：在满足可靠性约束下优化目标
风险导向设计：考虑失效后果的设计决策

维护决策：

基于可靠性的检查计划
剩余寿命预测的不确定性
维修策略的风险评估

1.3 传统方法的局限性

蒙特卡洛模拟（MCM）：

优点：概念简单、适用范围广、收敛速度与维度无关
缺点：收敛速度慢（ $O(1/N)O(1/\sqrt{N})$ ），高维问题计算成本巨大
对于小失效概率问题（ $10^{-6}$ 量级），需要数百万次模拟

泰勒展开法（FORM/SORM）：

优点：计算效率高，适合小失效概率问题
缺点：需要计算梯度，对非线性问题精度有限
难以处理多模态、高度非线性的极限状态函数

1.4 先进方法概述

本教程介绍的先进不确定性量化方法包括：

谱方法：

多项式混沌展开（PCE）
Karhunen-Loève展开（KLE）
广义混沌多项式（gPC）

代理模型方法：

高斯过程回归（Kriging）
径向基函数（RBF）
支持向量机（SVM）

贝叶斯方法：

贝叶斯参数估计
贝叶斯模型更新
马尔可夫链蒙特卡洛（MCMC）

降维技术：

主动子空间法
充分降维子空间
主成分分析（PCA）

2. 不确定性量化的数学基础

2.1 概率空间与随机变量

2.1.1 概率空间

概率空间由三元组 $(Ω,F,P)(\Omega, \mathcal{F}, P)$ 定义：

$Ω\Omega$ ：样本空间，包含所有可能结果
$F\mathcal{F}$ ： $σ\sigma$ -代数，定义可测事件集合
$P$ ：概率测度， $\mathcal{F} \rightarrow [0,1]$

2.1.2 随机变量及其分布

随机变量 $\Omega \rightarrow \mathbb{R}$ 将样本映射到实数。

常用概率分布：

分布类型	概率密度函数	参数	应用场景
正态分布	$\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$	$μ,σ\mu, \sigma$	材料参数、测量误差
对数正态	$\frac{1}{x\sigma\sqrt{2\pi}}e^{-\frac{(\ln x-\mu)^2}{2\sigma^2}}$	$μ,σ\mu, \sigma$	疲劳寿命、裂纹尺寸
威布尔分布	$\frac{k}{\lambda}(\frac{x}{\lambda})^{k-1}e^{-(x/\lambda)^k}$	$\lambda$	强度、寿命
均匀分布	$\frac{1}{b-a}$	$a, b$	几何误差、模型误差
指数分布	$\lambda e^{-\lambda x}$	$λ\lambda$	泊松过程、等待时间

2.1.3 随机过程与随机场

随机过程：参数（通常是时间）的函数，其取值是随机变量
$\omega), \quad t \in T, \omega \in \Omega$

随机场：空间参数的函数
$\omega), \quad x \in D \subset \mathbb{R}^d$

统计特性：

均值函数： $μ(x)=E[H(x,⋅)]\mu(x) = E[H(x, \cdot)]$
协方差函数： $\cdot), H(y, \cdot)]$

2.2 不确定性传播

2.2.1 问题描述

给定输入不确定性 $X$ 和计算模型 $Y = g (X)$ ，求输出 $Y$ 的统计特性。

前向不确定性传播：

已知： $X$ 的分布 $f_X(x)$
求： $Y$ 的分布 $f_Y(y)$ 或统计矩

逆向不确定性量化：

已知：观测数据 $D = \{y_1, y_2, ..., y_n\}$
求：模型参数 $θ\theta$ 的后验分布 $p(θ∣D)p(\theta|D)$

2.2.2 统计矩计算

均值（一阶矩）：
$μY=E[Y]=∫Rdg(x)fX(x)dx\mu_Y = E[Y] = \int_{\mathbb{R}^d} g(x) f_X(x) dx$

方差（二阶中心矩）：
$σY2=Var[Y]=E[(Y−μY)2]\sigma_Y^2 = Var[Y] = E[(Y - \mu_Y)^2]$

高阶矩：

偏度（三阶标准化矩）： $γ1=E[(Y−μY)3]/σY3\gamma_1 = E[(Y-\mu_Y)^3]/\sigma_Y^3$
峰度（四阶标准化矩）： $γ2=E[(Y−μY)4]/σY4−3\gamma_2 = E[(Y-\mu_Y)^4]/\sigma_Y^4 - 3$

2.3 可靠性分析基础

2.3.1 失效概率定义

$Pf=P(g(X)≤0)=∫g(x)≤0fX(x)dxP_f = P(g(X) \leq 0) = \int_{g(x) \leq 0} f_X(x) dx$

其中， $g (X)$ 是极限状态函数：

$g (X) > 0$ ：安全域
$g (X) = 0$ ：极限状态面
$g (X) < 0$ ：失效域

2.3.2 可靠性指标

Hasofer-Lind可靠性指标：
$β=min⁡g(x)=0(x−μ)TΣ−1(x−μ)\beta = \min_{g(x)=0} \sqrt{(x-\mu)^T \Sigma^{-1} (x-\mu)}$

失效概率近似：
$Pf≈Φ(−β)P_f \approx \Phi(-\beta)$

其中， $Φ\Phi$ 是标准正态分布的累积分布函数。

3. 多项式混沌展开方法

3.1 多项式混沌理论基础

3.1.1 Wiener混沌展开

N. Wiener（1938）提出用Hermite多项式展开高斯随机过程：

$Y(ξ)=∑α∈JyαHα(ξ)Y(\xi) = \sum_{\alpha \in \mathcal{J}} y_\alpha H_\alpha(\xi)$

其中：

$ξ=(ξ1,ξ2,...)\xi = (\xi_1, \xi_2, ...)$ 是标准正态随机变量
$HαH_\alpha$ 是多维Hermite多项式
$yαy_\alpha$ 是混沌系数

3.1.2 广义多项式混沌（gPC）

Xiu & Karniadakis（2002）将PCE推广到非高斯随机变量：

分布类型	正交多项式	支撑集
正态分布	Hermite	$(−∞,+∞)(-\infty, +\infty)$
均匀分布	Legendre	$[a, b]$
伽马分布	Laguerre	$+\infty)$
Beta分布	Jacobi	$[a, b]$
泊松分布	Charlier	${0, 1, 2, ...\}$

Askey方案：根据输入分布选择最优正交多项式基。

3.2 多项式混沌展开的实现

3.2.1 一维PCE展开

对于随机变量 $ξ∼N(0,1)\xi \sim N(0,1)$ ：

$Y(ξ)≈∑i=0PyiHi(ξ)Y(\xi) \approx \sum_{i=0}^{P} y_i H_i(\xi)$

Hermite多项式：

$H0(ξ)=1H_0(\xi) = 1$
$H1(ξ)=ξH_1(\xi) = \xi$
$H2(ξ)=ξ2−1H_2(\xi) = \xi^2 - 1$
$H3(ξ)=ξ3−3ξH_3(\xi) = \xi^3 - 3\xi$
$H4(ξ)=ξ4−6ξ2+3H_4(\xi) = \xi^4 - 6\xi^2 + 3$

递归关系：
$Hn+1(ξ)=ξHn(ξ)−nHn−1(ξ)H_{n+1}(\xi) = \xi H_n(\xi) - n H_{n-1}(\xi)$

3.2.2 多维PCE展开

对于 $d$ 维随机输入 $ξ=(ξ1,...,ξd)\xi = (\xi_1, ..., \xi_d)$ ：

$Y(ξ)≈∑α∈Ap,dyαΨα(ξ)Y(\xi) \approx \sum_{\alpha \in \mathcal{A}_{p,d}} y_\alpha \Psi_\alpha(\xi)$

其中， $α=(α1,...,αd)\alpha = (\alpha_1, ..., \alpha_d)$ 是多重指标， $∣α∣=∑αi≤p|\alpha| = \sum \alpha_i \leq p$ 。

总项数：
$NPCE=(d+pp)=(d+p)!d!p!N_{PCE} = \binom{d+p}{p} = \frac{(d+p)!}{d!p!}$

维数灾难：当 $d$ 较大时， $N_{PCE}$ 急剧增加。

3.2.3 混沌系数的计算

投影法（Galerkin投影）：
$yα=E[Y(ξ)Ψα(ξ)]E[Ψα2(ξ)]y_\alpha = \frac{E[Y(\xi)\Psi_\alpha(\xi)]}{E[\Psi_\alpha^2(\xi)]}$

配点法（Collocation）：

选择一组配点 ${ξ(i)}i=1N\{\xi^{(i)}\}_{i=1}^N$
在每个配点计算 $Y(i)=g(ξ(i))Y^{(i)} = g(\xi^{(i)})$
通过最小二乘拟合混沌系数：
$min⁡yα∑i=1N(Y(i)−∑αyαΨα(ξ(i)))2\min_{y_\alpha} \sum_{i=1}^N \left(Y^{(i)} - \sum_\alpha y_\alpha \Psi_\alpha(\xi^{(i)})\right)^2$

3.3 PCE的统计矩计算

3.3.1 均值和方差

均值：
$μY=E[Y]=y0\mu_Y = E[Y] = y_0$

方差：
$σY2=Var[Y]=∑α≠0yα2E[Ψα2]\sigma_Y^2 = Var[Y] = \sum_{\alpha \neq 0} y_\alpha^2 E[\Psi_\alpha^2]$

对于标准正交多项式， $E[Ψα2]=α!E[\Psi_\alpha^2] = \alpha!$ 。

3.3.2 敏感性分析（Sobol指数）

总方差分解：
$σY2=∑iσi2+∑i<jσij2+...\sigma_Y^2 = \sum_{i} \sigma_i^2 + \sum_{i<j} \sigma_{ij}^2 + ...$

一阶Sobol指数：
$Si=σi2σY2S_i = \frac{\sigma_i^2}{\sigma_Y^2}$

表示第 $i$ 个输入变量对输出方差的贡献。

总Sobol指数：
$SiT=σi2+∑j≠iσij2+...σY2S_i^T = \frac{\sigma_i^2 + \sum_{j \neq i}\sigma_{ij}^2 + ...}{\sigma_Y^2}$

表示第 $i$ 个输入变量及其所有交互作用的贡献。

基于PCE的Sobol指数计算：
$Si=∑α∈Iiyα2E[Ψα2]∑α≠0yα2E[Ψα2]S_i = \frac{\sum_{\alpha \in \mathcal{I}_i} y_\alpha^2 E[\Psi_\alpha^2]}{\sum_{\alpha \neq 0} y_\alpha^2 E[\Psi_\alpha^2]}$

其中， $j≠i}\mathcal{I}_i = \{\alpha: \alpha_i > 0, \alpha_j = 0 \text{ for } j \neq i\}$ 。

4. Karhunen-Loève展开

4.1 随机场的谱表示

4.1.1 KL展开的定义

Karhunen-Loève展开将随机场表示为确定性函数与随机系数的线性组合：

$\omega) = \mu(x) + \sum_{i=1}^{\infty} \sqrt{\lambda_i} \phi_i(x) \xi_i(\omega)$

其中：

$λi\lambda_i$ 和 $ϕi(x)\phi_i(x)$ 是协方差函数 $C (x, y)$ 的特征值和特征函数
$ξi\xi_i$ 是标准正态随机变量

4.1.2 特征值问题

协方差核的特征值问题：
$∫DC(x,y)ϕi(y)dy=λiϕi(x)\int_D C(x, y) \phi_i(y) dy = \lambda_i \phi_i(x)$

Mercer定理：
$\sum_{i=1}^{\infty} \lambda_i \phi_i(x) \phi_i(y)$

4.1.3 常用协方差核

指数型（Markov核）：
$\sigma^2 \exp\left(-\frac{|x-y|}{l_c}\right)$

高斯型（平方指数）：
$\sigma^2 \exp\left(-\frac{|x-y|^2}{2l_c^2}\right)$

Matérn类：
$\sigma^2 \frac{2^{1-\nu}}{\Gamma(\nu)} \left(\frac{\sqrt{2\nu}|x-y|}{l_c}\right)^\nu K_\nu\left(\frac{\sqrt{2\nu}|x-y|}{l_c}\right)$

其中， $l_c$ 是相关长度， $ν\nu$ 控制平滑性。

4.2 KL展开的数值实现

4.2.1 离散化方法

Nystrom方法：

将积分域离散为网格 ${x_i\}_{i=1}^n$
构造协方差矩阵 $C_{ij} = C(x_i, x_j)$
求解特征值问题： $\phi = \lambda \phi$

有限元方法：

使用有限元基函数近似特征函数
求解广义特征值问题

4.2.2 截断准则

能量截断：
选择 $M$ 使得前 $M$ 个特征值的能量占比达到阈值：
$∑i=1Mλi∑i=1∞λi≥ϵ\frac{\sum_{i=1}^M \lambda_i}{\sum_{i=1}^{\infty} \lambda_i} \geq \epsilon$

通常取 $ϵ=0.95\epsilon = 0.95$ 或 $0.99$ 。

误差估计：
$E[(H(x)−HM(x))2]=∑i=M+1∞λiϕi2(x)E\left[\left(H(x) - H_M(x)\right)^2\right] = \sum_{i=M+1}^{\infty} \lambda_i \phi_i^2(x)$

4.3 KL-PCE耦合方法

对于具有空间变异性的随机输入，可以结合KL展开和PCE：

KL展开降维：将随机场 $\omega)$ 展开为有限个随机变量 $ξ1,...,ξM\xi_1, ..., \xi_M$
PCE展开：将输出 $Y$ 展开为 $ξ\xi$ 的多项式混沌

优势：

KL展开处理空间相关性
PCE处理非线性映射
组合方法可有效处理高维问题

5. 高斯过程回归

5.1 高斯过程基础

5.1.1 定义

高斯过程是随机变量的集合，其中任意有限个随机变量的联合分布都是多元高斯分布。

$\sim GP(\mu(x), k(x, x'))$

其中：

$μ(x)\mu(x)$ ：均值函数
$k (x, x^{'})$ ：协方差函数（核函数）

5.1.2 核函数

平方指数核（SE/RBF）：
$\sigma_f^2 \exp\left(-\frac{\|x-x'\|^2}{2l^2}\right)$

Matérn核：
$\sigma_f^2 \frac{2^{1-\nu}}{\Gamma(\nu)} \left(\frac{\sqrt{2\nu}\|x-x'\|}{l}\right)^\nu K_\nu\left(\frac{\sqrt{2\nu}\|x-x'\|}{l}\right)$

有理二次核：
$\sigma_f^2 \left(1 + \frac{\|x-x'\|^2}{2\alpha l^2}\right)^{-\alpha}$

周期核：
$\sigma_f^2 \exp\left(-\frac{2\sin^2(\pi|x-x'|/p)}{l^2}\right)$

5.2 高斯过程回归（Kriging）

5.2.1 回归模型

给定训练数据 ${(x_i, y_i)\}_{i=1}^n$ ，其中 $yi=f(xi)+ϵiy_i = f(x_i) + \epsilon_i$ ， $ϵi∼N(0,σn2)\epsilon_i \sim N(0, \sigma_n^2)$ 。

预测分布：
对于新输入 $x_*$ ，预测输出 $f(x_*)$ 的后验分布为：

$f(x∗)∣X,y,x∗∼N(μ(x∗),σ2(x∗))f(x_*) | X, y, x_* \sim N(\mu(x_*), \sigma^2(x_*))$

其中：
$μ(x∗)=k∗T(K+σn2I)−1y\mu(x_*) = k_*^T (K + \sigma_n^2 I)^{-1} y$
$σ2(x∗)=k(x∗,x∗)−k∗T(K+σn2I)−1k∗\sigma^2(x_*) = k(x_*, x_*) - k_*^T (K + \sigma_n^2 I)^{-1} k_*$

$K_{ij} = k(x_i, x_j)$ ， $k_{*i} = k(x_*, x_i)$ 。

5.2.2 超参数优化

通过最大化对数似然函数估计核函数超参数：

$log⁡p(y∣X,θ)=−12yTKθ−1y−12log⁡∣Kθ∣−n2log⁡(2π)\log p(y|X, \theta) = -\frac{1}{2}y^T K_\theta^{-1} y - \frac{1}{2}\log|K_\theta| - \frac{n}{2}\log(2\pi)$

优化方法：

共轭梯度法
L-BFGS-B
遗传算法

5.3 GP在不确定性量化中的应用

5.3.1 代理模型

优势：

提供预测均值和方差
天然的不确定性估计
适合小样本学习

自适应采样：

基于预测方差选择新样本点
平衡探索与利用

5.3.2 贝叶斯优化

采集函数：

期望改进（EI）：
$EI(x) = E[\max(0, f(x) - f^+)]$
概率改进（PI）：
$PI(x) = P(f(x) > f^+)$
上置信界（UCB）：
$\mu(x) + \kappa \sigma(x)$

其中， $f^+$ 是当前最优值。

6. 贝叶斯推断与MCMC

6.1 贝叶斯理论基础

6.1.1 贝叶斯定理

$p(θ∣D)=p(D∣θ)p(θ)p(D)p(\theta|D) = \frac{p(D|\theta) p(\theta)}{p(D)}$

其中：

$p(θ)p(\theta)$ ：先验分布
$p(D∣θ)p(D|\theta)$ ：似然函数
$p(θ∣D)p(\theta|D)$ ：后验分布
$\int p(D|\theta) p(\theta) d\theta$ ：证据（边缘似然）

6.1.2 先验分布的选择

无信息先验：

均匀先验： $p(θ)∝1p(\theta) \propto 1$
Jeffreys先验： $p(θ)∝∣I(θ)∣p(\theta) \propto \sqrt{|I(\theta)|}$

共轭先验：

似然	共轭先验	后验
正态（已知方差）	正态	正态
正态（已知均值）	逆伽马	逆伽马
二项	Beta	Beta
泊松	伽马	伽马

信息先验：

基于专家知识
基于历史数据

6.2 马尔可夫链蒙特卡洛（MCMC）

6.2.1 Metropolis-Hastings算法

算法步骤：

初始化 $θ(0)\theta^{(0)}$
对于 $i = 1, 2, ..., N$ ：
- 从提议分布 $q(θ∗∣θ(i−1))q(\theta^*|\theta^{(i-1)})$ 采样候选点
- 计算接受概率：
  $α=min⁡(1,p(θ∗∣D)q(θ(i−1)∣θ∗)p(θ(i−1)∣D)q(θ∗∣θ(i−1)))\alpha = \min\left(1, \frac{p(\theta^*|D) q(\theta^{(i-1)}|\theta^*)}{p(\theta^{(i-1)}|D) q(\theta^*|\theta^{(i-1)})}\right)$
- 以概率 $α\alpha$ 接受 $θ∗\theta^*$ ，否则保持 $θ(i−1)\theta^{(i-1)}$

提议分布：

随机游走： $q(θ∗∣θ)=N(θ,Σ)q(\theta^*|\theta) = N(\theta, \Sigma)$
独立采样： $q(θ∗∣θ)=q(θ∗)q(\theta^*|\theta) = q(\theta^*)$

6.2.2 Gibbs采样

对于多维参数 $θ=(θ1,...,θd)\theta = (\theta_1, ..., \theta_d)$ ：

初始化 $θ(0)\theta^{(0)}$
对于 $i = 1, 2, ..., N$ ：
- 从 $p(θ1∣θ2(i−1),...,θd(i−1),D)p(\theta_1|\theta_2^{(i-1)}, ..., \theta_d^{(i-1)}, D)$ 采样 $θ1(i)\theta_1^{(i)}$
- 从 $p(θ2∣θ1(i),θ3(i−1),...,D)p(\theta_2|\theta_1^{(i)}, \theta_3^{(i-1)}, ..., D)$ 采样 $θ2(i)\theta_2^{(i)}$
- …
- 从 $p(θd∣θ1(i),...,θd−1(i),D)p(\theta_d|\theta_1^{(i)}, ..., \theta_{d-1}^{(i)}, D)$ 采样 $θd(i)\theta_d^{(i)}$

优点：不需要选择提议分布，接受率总是1。

6.2.3 高级MCMC方法

Hamiltonian Monte Carlo（HMC）：

利用梯度信息指导采样
引入动量变量
使用辛积分器

No-U-Turn Sampler（NUTS）：

自动确定轨迹长度
避免手动调参
在Stan、PyMC3中使用

仿射不变MCMC（emcee）：

使用系综采样
对线性变换不变
适合多模态分布

6.3 贝叶斯模型更新

6.3.1 模型参数更新

给定观测数据 $D$ ，更新模型参数 $θ\theta$ ：

$p(θ∣D)∝p(D∣θ)p(θ)p(\theta|D) \propto p(D|\theta) p(\theta)$

应用：

材料参数识别
模型修正
损伤诊断

6.3.2 模型选择

贝叶斯因子：
$B12=p(D∣M1)p(D∣M2)=∫p(D∣θ1,M1)p(θ1∣M1)dθ1∫p(D∣θ2,M2)p(θ2∣M2)dθ2B_{12} = \frac{p(D|M_1)}{p(D|M_2)} = \frac{\int p(D|\theta_1, M_1) p(\theta_1|M_1) d\theta_1}{\int p(D|\theta_2, M_2) p(\theta_2|M_2) d\theta_2}$

信息准则：

AIC（Akaike信息准则）
BIC（贝叶斯信息准则）
DIC（偏差信息准则）

7. 高维问题的处理策略

7.1 维数灾难问题

7.1.1 问题描述

当输入维度 $d$ 增加时：

PCE的项数 $NPCE=(d+pp)N_{PCE} = \binom{d+p}{p}$ 指数增长
蒙特卡洛模拟的收敛速度不变，但需要更多样本覆盖高维空间
网格法的点数 $N = n^d$ 指数增长

维数灾难的直观理解：

在单位超立方体 $0,1]^d$ 中，要达到相同的采样密度，样本量随维度指数增长
高维空间中，数据稀疏，距离度量失效

7.1.2 降维技术

主成分分析（PCA）：

寻找数据方差最大的方向
将数据投影到低维子空间

核PCA：

使用核函数将数据映射到高维特征空间
在特征空间进行PCA

流形学习：

假设数据位于低维流形上
ISOMAP、LLE、t-SNE等方法

7.2 主动子空间法

7.2.1 基本概念

Paul Constantine（2015）提出的主动子空间法：

如果函数 $f (x)$ 的梯度主要沿某些方向变化，则可以识别出"主动子空间"。

矩阵定义：
$\int \nabla f(x) \nabla f(x)^T \rho(x) dx$

其中， $ρ(x)\rho(x)$ 是输入的概率密度。

7.2.2 特征分解

对 $C$ 进行特征分解：
$\Lambda W^T$

其中， $Λ=diag(λ1,...,λd)\Lambda = diag(\lambda_1, ..., \lambda_d)$ ， $λ1≥λ2≥...≥λd≥0\lambda_1 \geq \lambda_2 \geq ... \geq \lambda_d \geq 0$ 。

主动子空间：

前 $k$ 个特征向量 $W_1 = [w_1, ..., w_k]$ 张成主动子空间
剩余特征向量 $W_2 = [w_{k+1}, ..., w_d]$ 张成非主动子空间

函数近似：
$\approx g(W_1^T x)$

将 $d$ 维问题降为 $k$ 维问题。

7.2.3 数值估计

蒙特卡洛估计：
$\approx \frac{1}{N} \sum_{i=1}^N \nabla f(x^{(i)}) \nabla f(x^{(i)})^T$

使用PCE估计：
如果 $f (x)$ 有PCE展开，可以直接计算梯度的解析表达式。

7.3 充分降维子空间

7.3.1 充分降维的概念

寻找低维子空间 $\subset \mathbb{R}^d$ ，使得：
$\perp X | P_S X$

即给定投影 $P_S X$ ，输出 $Y$ 与输入 $X$ 条件独立。

7.3.2 sliced逆回归（SIR）

算法步骤：

将数据按 $Y$ 值分片（slice）
计算每个片内 $X$ 的均值
对这些均值进行主成分分析
主成分方向即为充分降维方向

8. Python仿真实现与案例分析

8.1 仿真程序架构

本节提供完整的不确定性量化先进方法Python程序，包含以下模块：

多项式混沌展开模块：PCE构建、统计矩计算、Sobol指数
Karhunen-Loève展开模块：随机场离散化、特征值问题
高斯过程回归模块：Kriging建模、超参数优化
贝叶斯推断模块：MCMC采样、后验分析
降维模块：主动子空间、PCA

8.2 案例1：多项式混沌展开

问题描述：
悬臂梁端部受集中力 $P$ 和分布载荷 $q$ 作用，分析端部位移的不确定性。

随机输入：

$\sim N(200, 10)$ GPa（弹性模量）
$\sim N(8.33\times10^{-6}, 4\times10^{-7})$ m⁴（截面惯性矩）
$\sim N(1000, 100)$ N（集中力）
$\sim N(500, 50)$ N/m（分布载荷）

输出：端部位移 $\frac{PL^3}{3EI} + \frac{qL^4}{8EI}$

仿真目标：

构建4维PCE模型
计算均值、方差、偏度、峰度
计算一阶和总Sobol指数
与蒙特卡洛结果对比

8.3 案例2：高斯过程回归

问题描述：
使用高斯过程构建悬臂梁位移的代理模型。

训练数据：

使用拉丁超立方采样生成50个训练点
在4维输入空间 $(E, I, P, q)$ 中采样

测试与验证：

在100个测试点上评估GP预测
计算预测误差（RMSE、MAE）
分析预测方差与真实误差的关系
自适应采样：选择方差最大的点加入训练集

8.4 案例3：贝叶斯参数识别

问题描述：
通过悬臂梁的位移测量数据，识别弹性模量 $E$ 。

观测数据：

测量位移 $wobs=wtrue+ϵw_{obs} = w_{true} + \epsilon$ ， $ϵ∼N(0,σnoise)\epsilon \sim N(0, \sigma_{noise})$
生成10组观测数据

先验分布：

$\sim N(200, 50)$ GPa（较宽的先验）

MCMC采样：

使用Metropolis-Hastings算法采样后验
计算后验均值和95%置信区间
与真实值对比
分析收敛性（迹线图、自相关图）

8.5 案例4：主动子空间分析

问题描述：
分析悬臂梁位移对各输入参数的敏感性，识别主动子空间。

分析步骤：

计算梯度协方差矩阵 $C$
特征分解，绘制特征值谱
确定主动子空间维度 $k$
在主动子空间上构建响应面
对比全维模型和降维模型的精度

9. 工程应用与案例分析

9.1 航空航天结构可靠性

应用背景：
飞机机翼在气动载荷、材料分散性和制造误差影响下的可靠性评估。

不确定性来源：

材料性能：复合材料弹性常数的分散性
几何参数：翼型厚度、铺层角度的制造误差
载荷条件：阵风载荷的随机性

分析方法：

使用PCE建立气动弹性响应的代理模型
计算失效概率（强度失效、颤振边界）
识别关键不确定性来源（Sobol指数）
基于可靠性的设计优化

9.2 核设备概率安全评估

应用背景：
核反应堆压力容器在辐照脆化和热冲击下的断裂风险评估。

不确定性来源：

材料性能：断裂韧性 $K_{IC}$ 的辐照退化
缺陷尺寸：无损检测的不确定性
载荷条件：失水事故的温度和压力瞬态

分析方法：

贝叶斯更新：基于检测数据更新缺陷分布
MCMC采样：计算失效概率的后验分布
敏感性分析：识别对风险贡献最大的因素
检查间隔优化：基于风险的在役检查计划

9.3 风电叶片疲劳寿命预测

应用背景：
风力机叶片在变幅载荷、材料分散性和环境因素影响下的疲劳寿命预测。

不确定性来源：

材料性能：S-N曲线的分散性
载荷谱：风速分布、湍流强度
损伤累积：Miner准则的模型误差

分析方法：

使用GP建立疲劳寿命的代理模型
考虑认知不确定性：模型形式误差
贝叶斯更新：基于监测数据修正寿命预测
维护决策：考虑不确定性的最优更换策略

10. 总结与展望

10.1 关键技术总结

本教程系统介绍了不确定性量化的先进方法：

谱方法：

多项式混沌展开：高效计算统计矩和敏感性指数
Karhunen-Loève展开：处理空间相关随机场
组合方法：KL-PCE处理高维空间相关问题

代理模型方法：

高斯过程回归：提供预测不确定性估计
自适应采样：高效构建代理模型
贝叶斯优化：全局优化与不确定性量化结合

贝叶斯方法：

MCMC采样：从复杂后验分布中采样
贝叶斯模型更新：融合观测数据与先验知识
模型选择：贝叶斯因子与信息准则

降维技术：

主动子空间法：识别重要输入方向
充分降维：保持输入-输出关系的最小子空间
缓解维数灾难，提高计算效率

10.2 发展趋势

深度学习方法：

深度神经网络作为代理模型
物理信息神经网络（PINN）融合物理约束
生成模型（VAE、GAN）用于不确定性量化

多保真度方法：

结合高保真和低保真模型的信息
多保真度高斯过程
自适应多保真度采样

实时不确定性量化：

数字孪生中的实时更新
流式数据处理方法
边缘计算与不确定性量化

10.3 学习建议

夯实概率统计基础：理解概率论、数理统计的基本概念
掌握数值方法：熟悉蒙特卡洛、数值积分等方法
动手实践：通过Python实现各种算法，加深理解
关注前沿进展：跟踪机器学习与不确定性量化的交叉研究
工程应用：结合实际工程问题，体会方法的价值

参考文献

Xiu, D. and Karniadakis, G.E. (2002). The Wiener-Askey polynomial chaos for stochastic differential equations. SIAM Journal on Scientific Computing, 24(2), pp.619-644.
Ghanem, R.G. and Spanos, P.D. (1991). Stochastic Finite Elements: A Spectral Approach. Springer-Verlag.
Sudret, B. (2008). Global sensitivity analysis using polynomial chaos expansions. Reliability Engineering & System Safety, 93(7), pp.964-979.
Rasmussen, C.E. and Williams, C.K.I. (2006). Gaussian Processes for Machine Learning. MIT Press.
Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A. and Rubin, D.B. (2013). Bayesian Data Analysis. 3rd ed. CRC Press.
Constantine, P.G. (2015). Active Subspaces: Emerging Ideas for Dimension Reduction in Parameter Studies. SIAM.
Sullivan, T.J. (2015). Introduction to Uncertainty Quantification. Springer.
Smith, R.C. (2013). Uncertainty Quantification: Theory, Implementation, and Applications. SIAM.
Najm, H.N. (2009). Uncertainty quantification and polynomial chaos techniques in computational fluid dynamics. Annual Review of Fluid Mechanics, 41, pp.35-52.
Betz, W., Papaioannou, I. and Straub, D. (2014). Numerical methods for the discretization of random fields by means of the Karhunen-Loève expansion. Computer Methods in Applied Mechanics and Engineering, 271, pp.109-129.

附录：Python仿真程序

"""
主题092：不确定性量化的先进方法
Python仿真程序

包含以下功能模块：
1. 多项式混沌展开 (PCE)
2. Karhunen-Loève展开 (KLE)
3. 高斯过程回归 (GPR)
4. 贝叶斯推断与MCMC
5. 主动子空间分析

作者: 强度仿真专家
日期: 2026-02-27
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats, integrate, optimize, linalg
from scipy.special import factorial, comb, gamma, kv
from scipy.interpolate import interp1d
import warnings
warnings.filterwarnings('ignore')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['font.size'] = 10

# ==================== 1. 多项式混沌展开类 ====================

class PolynomialChaosExpansion:
    """
    多项式混沌展开 (PCE) 类
    
    用于不确定性传播和敏感性分析
    """
    
    def __init__(self, dim, order, dist_types=None):
        """
        初始化PCE
        
        参数:
            dim: 输入维度
            order: 多项式阶数
            dist_types: 分布类型列表，如 ['normal', 'uniform', ...]
        """
        self.dim = dim
        self.order = order
        self.dist_types = dist_types if dist_types else ['normal'] * dim
        
        # 生成多指数集合
        self.multi_indices = self._generate_multi_indices()
        self.n_terms = len(self.multi_indices)
        
        # 混沌系数
        self.coefficients = None
        
    def _generate_multi_indices(self):
        """生成总阶数不超过order的多重指标集合"""
        indices = []
        
        def recursive_generate(current, remaining_order, position):
            if position == self.dim:
                if remaining_order >= 0:
                    indices.append(tuple(current))
                return
            
            for k in range(remaining_order + 1):
                recursive_generate(current + [k], remaining_order - k, position + 1)
        
        recursive_generate([], self.order, 0)
        return indices
    
    def _hermite_polynomial(self, n, xi):
        """计算n阶Hermite多项式 (物理学家版本，对应标准正态)"""
        if n == 0:
            return np.ones_like(xi)
        elif n == 1:
            return xi
        else:
            # 使用递归关系: H_n(x) = x*H_{n-1}(x) - (n-1)*H_{n-2}(x)
            H_prev2 = np.ones_like(xi)
            H_prev1 = xi
            for i in range(2, n + 1):
                H_current = xi * H_prev1 - (i - 1) * H_prev2
                H_prev2 = H_prev1
                H_prev1 = H_current
            return H_current
    
    def _legendre_polynomial(self, n, xi):
        """计算n阶Legendre多项式 (定义在[-1,1])"""
        if n == 0:
            return np.ones_like(xi)
        elif n == 1:
            return xi
        else:
            # 递归关系
            P_prev2 = np.ones_like(xi)
            P_prev1 = xi
            for i in range(2, n + 1):
                P_current = ((2*i - 1) * xi * P_prev1 - (i - 1) * P_prev2) / i
                P_prev2 = P_prev1
                P_prev1 = P_current
            return P_current
    
    def _multivariate_polynomial(self, xi, alpha):
        """
        计算多维正交多项式
        
        参数:
            xi: 标准随机变量样本 [n_samples, dim]
            alpha: 多重指标
        
        返回:
            Psi: 多项式值 [n_samples]
        """
        n_samples = xi.shape[0] if xi.ndim > 1 else 1
        Psi = np.ones(n_samples)
        
        for i, (a, dist_type) in enumerate(zip(alpha, self.dist_types)):
            if a > 0:
                if dist_type == 'normal':
                    Psi *= self._hermite_polynomial(a, xi[:, i] if xi.ndim > 1 else xi[i])
                elif dist_type == 'uniform':
                    Psi *= self._legendre_polynomial(a, xi[:, i] if xi.ndim > 1 else xi[i])
        
        return Psi
    
    def _polynomial_norm(self, alpha):
        """计算多项式的L2范数平方 E[Psi_alpha^2]"""
        norm_sq = 1.0
        for a, dist_type in zip(alpha, self.dist_types):
            if dist_type == 'normal':
                # Hermite多项式: E[H_n^2] = n!
                norm_sq *= factorial(a)
            elif dist_type == 'uniform':
                # Legendre多项式: E[P_n^2] = 1/(2n+1)
                norm_sq *= 1.0 / (2*a + 1)
        return norm_sq
    
    def fit_collocation(self, model, n_samples, sampling='lhs'):
        """
        使用配点法拟合PCE系数
        
        参数:
            model: 计算模型函数 model(xi)
            n_samples: 样本数量
            sampling: 采样方法 'lhs' 或 'random'
        """
        # 生成样本点
        if sampling == 'lhs':
            xi_samples = self._lhs_sampling(n_samples)
        else:
            xi_samples = self._random_sampling(n_samples)
        
        # 计算模型响应
        y_samples = np.array([model(xi) for xi in xi_samples])
        
        # 构建设计矩阵
        Phi = np.zeros((n_samples, self.n_terms))
        for j, alpha in enumerate(self.multi_indices):
            Phi[:, j] = self._multivariate_polynomial(xi_samples, alpha)
        
        # 最小二乘求解
        self.coefficients, residuals, rank, s = linalg.lstsq(Phi, y_samples)
        
        # 计算R2分数
        y_pred = Phi @ self.coefficients
        ss_res = np.sum((y_samples - y_pred)**2)
        ss_tot = np.sum((y_samples - np.mean(y_samples))**2)
        self.r2_score = 1 - ss_res / ss_tot if ss_tot > 0 else 0
        
        return self.coefficients, self.r2_score
    
    def _lhs_sampling(self, n_samples):
        """拉丁超立方采样"""
        samples = np.zeros((n_samples, self.dim))
        for i in range(self.dim):
            # 分层采样
            perm = np.random.permutation(n_samples)
            samples[:, i] = (perm + np.random.rand(n_samples)) / n_samples
            
            # 转换到标准分布
            if self.dist_types[i] == 'normal':
                samples[:, i] = stats.norm.ppf(samples[:, i])
            elif self.dist_types[i] == 'uniform':
                samples[:, i] = 2 * samples[:, i] - 1  # 映射到[-1, 1]
        
        return samples
    
    def _random_sampling(self, n_samples):
        """随机采样"""
        samples = np.zeros((n_samples, self.dim))
        for i in range(self.dim):
            if self.dist_types[i] == 'normal':
                samples[:, i] = np.random.randn(n_samples)
            elif self.dist_types[i] == 'uniform':
                samples[:, i] = np.random.uniform(-1, 1, n_samples)
        return samples
    
    def predict(self, xi):
        """使用PCE模型预测"""
        if self.coefficients is None:
            raise ValueError("PCE模型尚未拟合")
        
        y_pred = 0
        for coef, alpha in zip(self.coefficients, self.multi_indices):
            y_pred += coef * self._multivariate_polynomial(xi.reshape(1, -1), alpha)[0]
        return y_pred
    
    def get_mean(self):
        """获取均值 (零阶系数)"""
        if self.coefficients is None:
            raise ValueError("PCE模型尚未拟合")
        return self.coefficients[0]
    
    def get_variance(self):
        """获取方差"""
        if self.coefficients is None:
            raise ValueError("PCE模型尚未拟合")
        
        variance = 0
        for coef, alpha in zip(self.coefficients[1:], self.multi_indices[1:]):
            variance += coef**2 * self._polynomial_norm(alpha)
        return variance
    
    def get_std(self):
        """获取标准差"""
        return np.sqrt(self.get_variance())
    
    def get_skewness(self):
        """获取偏度"""
        if self.coefficients is None:
            raise ValueError("PCE模型尚未拟合")
        
        mu = self.get_mean()
        sigma = self.get_std()
        
        # 三阶中心矩
        mu3 = 0
        for i, (coef_i, alpha_i) in enumerate(zip(self.coefficients, self.multi_indices)):
            for j, (coef_j, alpha_j) in enumerate(zip(self.coefficients, self.multi_indices)):
                for k, (coef_k, alpha_k) in enumerate(zip(self.coefficients, self.multi_indices)):
                    # E[Psi_i * Psi_j * Psi_k]
                    E_ijk = self._triple_product_expectation(alpha_i, alpha_j, alpha_k)
                    mu3 += coef_i * coef_j * coef_k * E_ijk
        
        mu3 -= 3*mu*sigma**2 + mu**3  # 转换为中心矩
        return mu3 / sigma**3
    
    def get_kurtosis(self):
        """获取超额峰度"""
        if self.coefficients is None:
            raise ValueError("PCE模型尚未拟合")
        
        mu = self.get_mean()
        sigma = self.get_std()
        
        # 四阶中心矩 (简化计算)
        # 使用蒙特卡洛从PCE模型采样
        n_mc = 10000
        samples = self.sample(n_mc)
        mu4 = np.mean((samples - mu)**4)
        
        return mu4 / sigma**4 - 3
    
    def _triple_product_expectation(self, alpha_i, alpha_j, alpha_k):
        """计算E[Psi_i * Psi_j * Psi_k] (简化版本)"""
        # 对于正交多项式，使用数值积分
        # 这里使用简化的高斯积分
        if self.dim == 1:
            # 一维情况
            xi_quad, w_quad = np.polynomial.hermite.hermgauss(10)
            Psi_i = self._hermite_polynomial(alpha_i[0], xi_quad)
            Psi_j = self._hermite_polynomial(alpha_j[0], xi_quad)
            Psi_k = self._hermite_polynomial(alpha_k[0], xi_quad)
            return np.sum(w_quad * Psi_i * Psi_j * Psi_k) / np.sqrt(np.pi)
        else:
            # 高维情况使用近似
            return 0.0
    
    def sample(self, n_samples):
        """从PCE模型采样"""
        # 生成标准随机变量
        xi_samples = self._random_sampling(n_samples)
        
        # 计算PCE响应
        y_samples = np.zeros(n_samples)
        for i in range(n_samples):
            y_samples[i] = self.predict(xi_samples[i])
        
        return y_samples
    
    def sobol_indices(self):
        """
        计算Sobol敏感性指数
        
        返回:
            S1: 一阶Sobol指数
            ST: 总Sobol指数
        """
        if self.coefficients is None:
            raise ValueError("PCE模型尚未拟合")
        
        variance = self.get_variance()
        if variance < 1e-15:
            return np.zeros(self.dim), np.zeros(self.dim)
        
        S1 = np.zeros(self.dim)
        ST = np.zeros(self.dim)
        
        for i in range(self.dim):
            # 一阶效应: 只有alpha_i > 0，其他为0
            var_i = 0
            var_total = 0
            
            for coef, alpha in zip(self.coefficients[1:], self.multi_indices[1:]):
                norm_sq = self._polynomial_norm(alpha)
                
                # 一阶: 只有第i个变量
                if alpha[i] > 0 and sum(alpha) == alpha[i]:
                    var_i += coef**2 * norm_sq
                
                # 总效应: 包含第i个变量的所有项
                if alpha[i] > 0:
                    var_total += coef**2 * norm_sq
            
            S1[i] = var_i / variance
            ST[i] = var_total / variance
        
        return S1, ST


# ==================== 2. Karhunen-Loève展开类 ====================

class KarhunenLoeveExpansion:
    """
    Karhunen-Loève展开类
    
    用于随机场的谱表示和降维
    """
    
    def __init__(self, domain, mean_func, cov_kernel):
        """
        初始化KL展开
        
        参数:
            domain: 空间域 [x_min, x_max]
            mean_func: 均值函数
            cov_kernel: 协方差核函数 cov_kernel(x, y)
        """
        self.domain = domain
        self.mean_func = mean_func
        self.cov_kernel = cov_kernel
        
        self.eigenvalues = None
        self.eigenfunctions = None
        
    def discretize_nystrom(self, n_points, n_modes):
        """
        使用Nystrom方法离散化
        
        参数:
            n_points: 离散点数
            n_modes: 保留的模态数
        """
        # 离散化空间
        x = np.linspace(self.domain[0], self.domain[1], n_points)
        dx = x[1] - x[0]
        
        # 构造协方差矩阵
        C = np.zeros((n_points, n_points))
        for i in range(n_points):
            for j in range(n_points):
                C[i, j] = self.cov_kernel(x[i], x[j])
        
        # 求解特征值问题
        # 注意: 需要乘以dx来近似积分
        eigenvalues, eigenvectors = linalg.eigh(C * dx)
        
        # 按特征值降序排列
        idx = np.argsort(eigenvalues)[::-1]
        eigenvalues = eigenvalues[idx]
        eigenvectors = eigenvectors[:, idx]
        
        # 只保留正特征值
        positive_idx = eigenvalues > 1e-10
        eigenvalues = eigenvalues[positive_idx]
        eigenvectors = eigenvectors[:, positive_idx]
        
        # 归一化特征函数
        for i in range(len(eigenvalues)):
            norm = np.sqrt(np.trapezoid(eigenvectors[:, i]**2, x))
            eigenvectors[:, i] /= norm
        
        # 保留前n_modes个模态
        n_modes = min(n_modes, len(eigenvalues))
        self.eigenvalues = eigenvalues[:n_modes]
        self.eigenfunctions = eigenvectors[:, :n_modes]
        self.x_grid = x
        self.n_modes = n_modes
        
        return self.eigenvalues, self.eigenfunctions
    
    def energy_ratio(self):
        """计算保留的能量比例"""
        if self.eigenvalues is None:
            raise ValueError("KL展开尚未离散化")
        return np.sum(self.eigenvalues) / np.sum(self.eigenvalues)
    
    def reconstruct(self, xi):
        """
        重构随机场样本
        
        参数:
            xi: 标准正态随机变量 [n_modes]
        
        返回:
            H: 随机场样本
        """
        if self.eigenvalues is None:
            raise ValueError("KL展开尚未离散化")
        
        H = self.mean_func(self.x_grid)
        for i in range(self.n_modes):
            H += np.sqrt(self.eigenvalues[i]) * self.eigenfunctions[:, i] * xi[i]
        
        return H
    
    def sample(self, n_samples):
        """
        生成随机场样本
        
        参数:
            n_samples: 样本数量
        
        返回:
            samples: 随机场样本 [n_samples, n_points]
        """
        samples = np.zeros((n_samples, len(self.x_grid)))
        for i in range(n_samples):
            xi = np.random.randn(self.n_modes)
            samples[i] = self.reconstruct(xi)
        return samples


# ==================== 3. 高斯过程回归类 ====================

class GaussianProcessRegression:
    """
    高斯过程回归 (Kriging) 类
    
    用于构建代理模型和不确定性估计
    """
    
    def __init__(self, kernel_type='rbf', noise_variance=1e-5):
        """
        初始化GP回归
        
        参数:
            kernel_type: 核函数类型 'rbf', 'matern', 'rational_quadratic'
            noise_variance: 观测噪声方差
        """
        self.kernel_type = kernel_type
        self.noise_variance = noise_variance
        self.X_train = None
        self.y_train = None
        self.K = None
        self.L = None
        self.alpha = None
        
        # 超参数
        self.length_scale = 1.0
        self.signal_variance = 1.0
        
    def _kernel(self, x1, x2):
        """计算核函数"""
        if self.kernel_type == 'rbf':
            # 平方指数核
            sq_dist = np.sum((x1 - x2)**2)
            return self.signal_variance * np.exp(-0.5 * sq_dist / self.length_scale**2)
        elif self.kernel_type == 'matern':
            # Matérn 3/2核
            dist = np.sqrt(np.sum((x1 - x2)**2))
            if dist == 0:
                return self.signal_variance
            sqrt3 = np.sqrt(3)
            return self.signal_variance * (1 + sqrt3 * dist / self.length_scale) * \
                   np.exp(-sqrt3 * dist / self.length_scale)
        else:
            # 默认RBF
            sq_dist = np.sum((x1 - x2)**2)
            return self.signal_variance * np.exp(-0.5 * sq_dist / self.length_scale**2)
    
    def _kernel_matrix(self, X1, X2=None):
        """计算核矩阵"""
        if X2 is None:
            X2 = X1
        
        n1 = X1.shape[0]
        n2 = X2.shape[0]
        K = np.zeros((n1, n2))
        
        for i in range(n1):
            for j in range(n2):
                K[i, j] = self._kernel(X1[i], X2[j])
        
        return K
    
    def fit(self, X, y, optimize_hyperparams=True):
        """
        训练GP模型
        
        参数:
            X: 训练输入 [n_samples, n_features]
            y: 训练输出 [n_samples]
            optimize_hyperparams: 是否优化超参数
        """
        self.X_train = np.array(X)
        self.y_train = np.array(y)
        
        if optimize_hyperparams:
            self._optimize_hyperparameters()
        
        # 计算核矩阵
        self.K = self._kernel_matrix(self.X_train)
        self.K += self.noise_variance * np.eye(len(self.y_train))
        
        # Cholesky分解
        try:
            self.L = linalg.cholesky(self.K, lower=True)
            self.alpha = linalg.solve_triangular(
                self.L.T, 
                linalg.solve_triangular(self.L, self.y_train, lower=True),
                lower=False
            )
        except np.linalg.LinAlgError:
            # 如果Cholesky失败，使用直接求解
            self.alpha = linalg.solve(self.K, self.y_train)
            self.L = None
    
    def _optimize_hyperparameters(self):
        """优化超参数"""
        def neg_log_likelihood(params):
            self.length_scale = np.exp(params[0])
            self.signal_variance = np.exp(params[1])
            
            K = self._kernel_matrix(self.X_train)
            K += self.noise_variance * np.eye(len(self.y_train))
            
            try:
                L = linalg.cholesky(K, lower=True)
                alpha = linalg.solve_triangular(
                    L.T,
                    linalg.solve_triangular(L, self.y_train, lower=True),
                    lower=False
                )
                log_likelihood = -0.5 * np.dot(self.y_train, alpha) - \
                                np.sum(np.log(np.diag(L))) - \
                                0.5 * len(self.y_train) * np.log(2 * np.pi)
                return -log_likelihood
            except:
                return 1e10
        
        # 使用简单的网格搜索
        best_nll = np.inf
        best_params = [0, 0]  # log(length_scale), log(signal_variance)
        
        for log_ls in np.linspace(-2, 2, 10):
            for log_sv in np.linspace(-2, 2, 10):
                nll = neg_log_likelihood([log_ls, log_sv])
                if nll < best_nll:
                    best_nll = nll
                    best_params = [log_ls, log_sv]
        
        self.length_scale = np.exp(best_params[0])
        self.signal_variance = np.exp(best_params[1])
    
    def predict(self, X_test, return_std=True):
        """
        预测
        
        参数:
            X_test: 测试输入 [n_test, n_features]
            return_std: 是否返回标准差
        
        返回:
            mu: 预测均值
            sigma: 预测标准差 (如果return_std=True)
        """
        X_test = np.array(X_test)
        if X_test.ndim == 1:
            X_test = X_test.reshape(1, -1)
        
        # 计算k_* = k(X_test, X_train)
        k_star = self._kernel_matrix(X_test, self.X_train)
        
        # 预测均值
        mu = k_star @ self.alpha
        
        if not return_std:
            return mu
        
        # 预测方差
        v = linalg.solve_triangular(self.L, k_star.T, lower=True) if self.L is not None else \
            linalg.solve(self.K, k_star.T)
        
        k_star_star = np.array([self._kernel(x, x) for x in X_test])
        var = k_star_star - np.sum(v**2, axis=0)
        var = np.maximum(var, 0)  # 确保非负
        
        return mu, np.sqrt(var)
    
    def sample_y(self, X_test, n_samples=1):
        """从后验分布采样"""
        mu, std = self.predict(X_test, return_std=True)
        return mu[:, None] + std[:, None] * np.random.randn(len(mu), n_samples)


# ==================== 4. 贝叶斯推断与MCMC类 ====================

class BayesianInference:
    """
    贝叶斯推断类
    
    使用MCMC方法进行参数估计
    """
    
    def __init__(self, log_prior, log_likelihood):
        """
        初始化贝叶斯推断
        
        参数:
            log_prior: 对数先验函数 log_prior(theta)
            log_likelihood: 对数似然函数 log_likelihood(theta, data)
        """
        self.log_prior = log_prior
        self.log_likelihood = log_likelihood
        self.samples = None
        
    def log_posterior(self, theta, data):
        """计算对数后验"""
        log_prior_val = self.log_prior(theta)
        if np.isinf(log_prior_val):
            return -np.inf
        return log_prior_val + self.log_likelihood(theta, data)
    
    def metropolis_hastings(self, data, theta_init, n_samples=10000, proposal_std=0.1, burn_in=2000):
        """
        Metropolis-Hastings采样
        
        参数:
            data: 观测数据
            theta_init: 初始参数值
            n_samples: 采样数量
            proposal_std: 提议分布的标准差
            burn_in: 预烧期样本数
        
        返回:
            samples: MCMC样本 (去除预烧期)
        """
        theta = np.array(theta_init)
        dim = len(theta)
        
        all_samples = np.zeros((n_samples + burn_in, dim))
        all_samples[0] = theta
        
        n_accepted = 0
        
        for i in range(1, n_samples + burn_in):
            # 提议新样本
            theta_proposal = theta + np.random.randn(dim) * proposal_std
            
            # 计算接受概率
            log_post_current = self.log_posterior(theta, data)
            log_post_proposal = self.log_posterior(theta_proposal, data)
            
            log_alpha = log_post_proposal - log_post_current
            alpha = np.exp(min(log_alpha, 0))
            
            # 接受或拒绝
            if np.random.rand() < alpha:
                theta = theta_proposal
                n_accepted += 1
            
            all_samples[i] = theta
        
        # 去除预烧期
        self.samples = all_samples[burn_in:]
        self.acceptance_rate = n_accepted / (n_samples + burn_in - 1)
        
        return self.samples
    
    def get_posterior_stats(self):
        """获取后验统计量"""
        if self.samples is None:
            raise ValueError("尚未进行MCMC采样")
        
        mean = np.mean(self.samples, axis=0)
        std = np.std(self.samples, axis=0)
        ci_lower = np.percentile(self.samples, 2.5, axis=0)
        ci_upper = np.percentile(self.samples, 97.5, axis=0)
        
        return {
            'mean': mean,
            'std': std,
            'ci_95': (ci_lower, ci_upper)
        }
    
    def geweke_diagnostic(self, first=0.1, last=0.5):
        """
        Geweke收敛诊断
        
        比较链的前段和后段的均值差异
        """
        if self.samples is None:
            raise ValueError("尚未进行MCMC采样")
        
        n = len(self.samples)
        first_samples = self.samples[:int(first * n)]
        last_samples = self.samples[int((1 - last) * n):]
        
        z_scores = []
        for i in range(self.samples.shape[1]):
            mean_first = np.mean(first_samples[:, i])
            mean_last = np.mean(last_samples[:, i])
            var_first = np.var(first_samples[:, i])
            var_last = np.var(last_samples[:, i])
            
            z = (mean_first - mean_last) / np.sqrt(var_first / len(first_samples) + var_last / len(last_samples))
            z_scores.append(z)
        
        return np.array(z_scores)


# ==================== 5. 主动子空间分析类 ====================

class ActiveSubspace:
    """
    主动子空间分析类
    
    用于高维问题的降维
    """
    
    def __init__(self, model, grad_model=None):
        """
        初始化主动子空间分析
        
        参数:
            model: 计算模型函数
            grad_model: 梯度计算函数 (可选)
        """
        self.model = model
        self.grad_model = grad_model
        self.eigenvalues = None
        self.eigenvectors = None
        
    def compute_gradient_covariance(self, X_samples, pdf_func=None):
        """
        计算梯度协方差矩阵
        
        参数:
            X_samples: 输入样本 [n_samples, dim]
            pdf_func: 概率密度函数 (可选)
        """
        n_samples, dim = X_samples.shape
        C = np.zeros((dim, dim))
        
        for x in X_samples:
            # 计算梯度
            if self.grad_model is not None:
                grad = self.grad_model(x)
            else:
                grad = self._numerical_gradient(x)
            
            # 权重 (如果提供pdf)
            weight = pdf_func(x) if pdf_func else 1.0
            
            C += np.outer(grad, grad) * weight
        
        C /= n_samples
        self.C = C
        
        # 特征分解
        self.eigenvalues, self.eigenvectors = linalg.eigh(C)
        
        # 按降序排列
        idx = np.argsort(self.eigenvalues)[::-1]
        self.eigenvalues = self.eigenvalues[idx]
        self.eigenvectors = self.eigenvectors[:, idx]
        
        return self.eigenvalues, self.eigenvectors
    
    def _numerical_gradient(self, x, h=1e-5):
        """数值计算梯度"""
        dim = len(x)
        grad = np.zeros(dim)
        
        for i in range(dim):
            x_plus = x.copy()
            x_minus = x.copy()
            x_plus[i] += h
            x_minus[i] -= h
            grad[i] = (self.model(x_plus) - self.model(x_minus)) / (2 * h)
        
        return grad
    
    def project_to_active(self, x, n_active):
        """
        投影到主动子空间
        
        参数:
            x: 输入向量
            n_active: 主动子空间维度
        
        返回:
            y: 主动变量
        """
        W1 = self.eigenvectors[:, :n_active]
        return W1.T @ x
    
    def reconstruct_from_active(self, y, n_active):
        """
        从主动子空间重构
        
        参数:
            y: 主动变量
            n_active: 主动子空间维度
        
        返回:
            x: 重构的输入向量
        """
        W1 = self.eigenvectors[:, :n_active]
        return W1 @ y
    
    def energy_ratio(self, n_active):
        """计算前n_active个特征值的能量占比"""
        return np.sum(self.eigenvalues[:n_active]) / np.sum(self.eigenvalues)


# ==================== 6. 仿真案例 ====================

def cantilever_beam_deflection(E, I, P, q, L=2.0):
    """
    悬臂梁端部位移计算
    
    参数:
        E: 弹性模量 (Pa)
        I: 截面惯性矩 (m^4)
        P: 端部集中力 (N)
        q: 分布载荷 (N/m)
        L: 梁长度 (m)
    
    返回:
        w: 端部位移 (m)
    """
    w_point = P * L**3 / (3 * E * I)
    w_dist = q * L**4 / (8 * E * I)
    return w_point + w_dist


def case1_pce_analysis():
    """案例1: 多项式混沌展开不确定性分析"""
    print("\n" + "="*70)
    print("案例1: 多项式混沌展开不确定性分析")
    print("="*70)
    
    # 问题参数
    L = 2.0  # 梁长度 (m)
    
    # 随机输入参数 (均值和标准差)
    E_mean, E_std = 200e9, 10e9  # 弹性模量 (Pa)
    I_mean, I_std = 8.33e-6, 4e-7  # 截面惯性矩 (m^4)
    P_mean, P_std = 1000, 100  # 集中力 (N)
    q_mean, q_std = 500, 50  # 分布载荷 (N/m)
    
    # 标准化函数
    def to_standard(x, mean, std):
        """转换到标准正态变量"""
        return (x - mean) / std
    
    def from_standard(xi, mean, std):
        """从标准正态变量转换"""
        return mean + xi * std
    
    # 定义计算模型 (输入为标准正态变量)
    def model(xi):
        E = from_standard(xi[0], E_mean, E_std)
        I = from_standard(xi[1], I_mean, I_std)
        P = from_standard(xi[2], P_mean, P_std)
        q = from_standard(xi[3], q_mean, q_std)
        return cantilever_beam_deflection(E, I, P, q, L)
    
    print("\n  问题描述:")
    print(f"    - 悬臂梁长度: {L} m")
    print(f"    - 弹性模量: 均值={E_mean/1e9:.1f} GPa, 标准差={E_std/1e9:.1f} GPa")
    print(f"    - 截面惯性矩: 均值={I_mean*1e6:.3f}e-6 m⁴, 标准差={I_std*1e6:.3f}e-6 m⁴")
    print(f"    - 集中力: 均值={P_mean} N, 标准差={P_std} N")
    print(f"    - 分布载荷: 均值={q_mean} N/m, 标准差={q_std} N/m")
    
    # 创建PCE模型
    dim = 4
    order = 3
    pce = PolynomialChaosExpansion(dim, order, dist_types=['normal']*dim)
    
    print(f"\n  PCE模型:")
    print(f"    - 输入维度: {dim}")
    print(f"    - 多项式阶数: {order}")
    print(f"    - 混沌项数: {pce.n_terms}")
    
    # 拟合PCE模型
    n_samples = 200
    coeffs, r2 = pce.fit_collocation(model, n_samples, sampling='lhs')
    
    print(f"\n  拟合结果:")
    print(f"    - 配点数: {n_samples}")
    print(f"    - R² 分数: {r2:.6f}")
    
    # 计算统计矩
    mean_pce = pce.get_mean()
    std_pce = pce.get_std()
    skew_pce = pce.get_skewness()
    kurt_pce = pce.get_kurtosis()
    
    print(f"\n  PCE统计矩:")
    print(f"    - 均值: {mean_pce*1000:.4f} mm")
    print(f"    - 标准差: {std_pce*1000:.4f} mm")
    print(f"    - 变异系数: {std_pce/mean_pce*100:.2f}%")
    print(f"    - 偏度: {skew_pce:.4f}")
    print(f"    - 超额峰度: {kurt_pce:.4f}")
    
    # 蒙特卡洛对比
    n_mc = 100000
    mc_samples = np.zeros(n_mc)
    for i in range(n_mc):
        xi = np.random.randn(dim)
        mc_samples[i] = model(xi)
    
    print(f"\n  蒙特卡洛对比 (N={n_mc}):")
    print(f"    - 均值: {np.mean(mc_samples)*1000:.4f} mm")
    print(f"    - 标准差: {np.std(mc_samples)*1000:.4f} mm")
    print(f"    - 偏度: {stats.skew(mc_samples):.4f}")
    print(f"    - 超额峰度: {stats.kurtosis(mc_samples):.4f}")
    
    # 计算Sobol指数
    S1, ST = pce.sobol_indices()
    
    print(f"\n  Sobol敏感性指数:")
    var_names = ['E (弹性模量)', 'I (惯性矩)', 'P (集中力)', 'q (分布载荷)']
    print(f"    {'变量':<20} {'一阶指数':<12} {'总指数':<12}")
    print(f"    {'-'*44}")
    for i, name in enumerate(var_names):
        print(f"    {name:<20} {S1[i]:<12.4f} {ST[i]:<12.4f}")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    
    # 1. PCE收敛性分析
    ax1 = axes[0, 0]
    sample_sizes = [50, 100, 150, 200, 300, 500]
    means_conv = []
    stds_conv = []
    
    for n in sample_sizes:
        pce_temp = PolynomialChaosExpansion(dim, order, dist_types=['normal']*dim)
        pce_temp.fit_collocation(model, n, sampling='lhs')
        means_conv.append(pce_temp.get_mean() * 1000)
        stds_conv.append(pce_temp.get_std() * 1000)
    
    ax1.plot(sample_sizes, means_conv, 'bo-', label='Mean (mm)')
    ax1.axhline(y=np.mean(mc_samples)*1000, color='b', linestyle='--', label='MC Mean')
    ax1.set_xlabel('Number of Samples', fontsize=11)
    ax1.set_ylabel('Mean Displacement (mm)', fontsize=11, color='b')
    ax1.tick_params(axis='y', labelcolor='b')
    ax1.grid(True, alpha=0.3)
    ax1.legend(loc='upper left')
    
    ax1_twin = ax1.twinx()
    ax1_twin.plot(sample_sizes, stds_conv, 'rs-', label='Std (mm)')
    ax1_twin.axhline(y=np.std(mc_samples)*1000, color='r', linestyle='--', label='MC Std')
    ax1_twin.set_ylabel('Std Displacement (mm)', fontsize=11, color='r')
    ax1_twin.tick_params(axis='y', labelcolor='r')
    ax1_twin.legend(loc='upper right')
    ax1.set_title('PCE Convergence Analysis', fontsize=12, fontweight='bold')
    
    # 2. Sobol指数
    ax2 = axes[0, 1]
    x_pos = np.arange(len(var_names))
    width = 0.35
    ax2.bar(x_pos - width/2, S1, width, label='First-order', color='steelblue')
    ax2.bar(x_pos + width/2, ST, width, label='Total', color='coral')
    ax2.set_xlabel('Input Variables', fontsize=11)
    ax2.set_ylabel('Sobol Index', fontsize=11)
    ax2.set_title('Sensitivity Analysis (Sobol Indices)', fontsize=12, fontweight='bold')
    ax2.set_xticks(x_pos)
    ax2.set_xticklabels(['E', 'I', 'P', 'q'], fontsize=10)
    ax2.legend()
    ax2.grid(True, alpha=0.3, axis='y')
    
    # 3. 概率密度对比
    ax3 = axes[1, 0]
    pce_samples = pce.sample(10000)
    ax3.hist(mc_samples*1000, bins=50, density=True, alpha=0.5, label='Monte Carlo', color='blue')
    ax3.hist(pce_samples*1000, bins=50, density=True, alpha=0.5, label='PCE Sample', color='red')
    ax3.set_xlabel('Displacement (mm)', fontsize=11)
    ax3.set_ylabel('Probability Density', fontsize=11)
    ax3.set_title('PDF Comparison', fontsize=12, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 4. 混沌系数
    ax4 = axes[1, 1]
    coef_abs = np.abs(coeffs)
    sorted_idx = np.argsort(coef_abs)[::-1][:20]  # 前20个
    ax4.barh(range(len(sorted_idx)), coef_abs[sorted_idx], color='green', alpha=0.7)
    ax4.set_xlabel('|Coefficient|', fontsize=11)
    ax4.set_ylabel('Polynomial Term Index', fontsize=11)
    ax4.set_title('PCE Coefficients (Top 20)', fontsize=12, fontweight='bold')
    ax4.set_yticks(range(len(sorted_idx)))
    ax4.set_yticklabels([f'{pce.multi_indices[i]}' for i in sorted_idx], fontsize=8)
    ax4.grid(True, alpha=0.3, axis='x')
    
    plt.tight_layout()
    plt.savefig('d:\\文档\\强度仿真\\主题092\\PCE不确定性分析.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("\n  可视化结果已保存: PCE不确定性分析.png")


def case2_gaussian_process():
    """案例2: 高斯过程回归代理模型"""
    print("\n" + "="*70)
    print("案例2: 高斯过程回归代理模型")
    print("="*70)
    
    # 问题参数
    L = 2.0
    E_mean, E_std = 200e9, 10e9
    I_mean, I_std = 8.33e-6, 4e-7
    P_mean, P_std = 1000, 100
    q_mean, q_std = 500, 50
    
    def to_standard(x, mean, std):
        return (x - mean) / std
    
    def from_standard(xi, mean, std):
        return mean + xi * std
    
    def model(xi):
        E = from_standard(xi[0], E_mean, E_std)
        I = from_standard(xi[1], I_mean, I_std)
        P = from_standard(xi[2], P_mean, P_std)
        q = from_standard(xi[3], q_mean, q_std)
        return cantilever_beam_deflection(E, I, P, q, L)
    
    print("\n  训练GP代理模型...")
    
    # 生成训练数据 (LHS)
    n_train = 30
    dim = 4
    X_train = np.zeros((n_train, dim))
    
    for i in range(dim):
        perm = np.random.permutation(n_train)
        X_train[:, i] = (perm + np.random.rand(n_train)) / n_train
        X_train[:, i] = stats.norm.ppf(X_train[:, i])
    
    y_train = np.array([model(xi) for xi in X_train])
    
    # 训练GP模型
    gp = GaussianProcessRegression(kernel_type='rbf', noise_variance=1e-10)
    gp.fit(X_train, y_train, optimize_hyperparams=True)
    
    print(f"\n  GP模型参数:")
    print(f"    - 训练样本数: {n_train}")
    print(f"    - 长度尺度: {gp.length_scale:.4f}")
    print(f"    - 信号方差: {gp.signal_variance:.4f}")
    
    # 生成测试数据
    n_test = 100
    X_test = np.random.randn(n_test, dim)
    y_true = np.array([model(xi) for xi in X_test])
    y_pred, y_std = gp.predict(X_test, return_std=True)
    
    # 计算误差指标
    mse = np.mean((y_true - y_pred)**2)
    rmse = np.sqrt(mse)
    mae = np.mean(np.abs(y_true - y_pred))
    r2 = 1 - np.sum((y_true - y_pred)**2) / np.sum((y_true - np.mean(y_true))**2)
    
    print(f"\n  测试集性能:")
    print(f"    - RMSE: {rmse*1000:.4f} mm")
    print(f"    - MAE: {mae*1000:.4f} mm")
    print(f"    - R²: {r2:.6f}")
    
    # 检查预测区间覆盖率
    in_interval = np.abs(y_true - y_pred) < 1.96 * y_std
    coverage = np.mean(in_interval)
    print(f"    - 95%置信区间覆盖率: {coverage*100:.1f}%")
    
    # 可视化 (仅显示前两个维度的切片)
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    
    # 1. 预测vs真实值
    ax1 = axes[0, 0]
    ax1.scatter(y_true*1000, y_pred*1000, alpha=0.6, c='blue', s=30)
    ax1.plot([y_true.min()*1000, y_true.max()*1000], 
             [y_true.min()*1000, y_true.max()*1000], 'r--', lw=2, label='Perfect Prediction')
    ax1.set_xlabel('True Displacement (mm)', fontsize=11)
    ax1.set_ylabel('Predicted Displacement (mm)', fontsize=11)
    ax1.set_title(f'GP Prediction vs True (R²={r2:.4f})', fontsize=12, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 2. 预测误差分布
    ax2 = axes[0, 1]
    errors = (y_true - y_pred) * 1000
    ax2.hist(errors, bins=20, density=True, alpha=0.7, color='green', edgecolor='black')
    ax2.axvline(x=0, color='r', linestyle='--', linewidth=2)
    ax2.set_xlabel('Prediction Error (mm)', fontsize=11)
    ax2.set_ylabel('Density', fontsize=11)
    ax2.set_title('Prediction Error Distribution', fontsize=12, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    
    # 3. 预测标准差vs绝对误差
    ax3 = axes[1, 0]
    ax3.scatter(y_std*1000, np.abs(errors), alpha=0.6, c='purple', s=30)
    ax3.plot([0, y_std.max()*1000], [0, 1.96*y_std.max()*1000], 'r--', lw=2, label='95% Bound')
    ax3.set_xlabel('Predicted Std (mm)', fontsize=11)
    ax3.set_ylabel('Absolute Error (mm)', fontsize=11)
    ax3.set_title('Uncertainty Calibration', fontsize=12, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 4. 一维切片可视化 (固定其他变量为0)
    ax4 = axes[1, 1]
    xi_range = np.linspace(-3, 3, 100)
    X_slice = np.zeros((100, dim))
    X_slice[:, 0] = xi_range  # 变化第一个变量
    
    y_slice_pred, y_slice_std = gp.predict(X_slice, return_std=True)
    
    # 计算真实值
    y_slice_true = np.array([model(xi) for xi in X_slice])
    
    ax4.plot(xi_range, y_slice_true*1000, 'b-', linewidth=2, label='True Model')
    ax4.plot(xi_range, y_slice_pred*1000, 'r--', linewidth=2, label='GP Prediction')
    ax4.fill_between(xi_range, 
                     (y_slice_pred - 2*y_slice_std)*1000,
                     (y_slice_pred + 2*y_slice_std)*1000,
                     alpha=0.3, color='red', label='95% CI')
    ax4.scatter(X_train[:, 0], y_train*1000, c='green', s=50, zorder=5, label='Training Data')
    ax4.set_xlabel('Standard Normal Variable ξ₁', fontsize=11)
    ax4.set_ylabel('Displacement (mm)', fontsize=11)
    ax4.set_title('1D Slice (Other Variables Fixed)', fontsize=12, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('d:\\文档\\强度仿真\\主题092\\GP代理模型.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("\n  可视化结果已保存: GP代理模型.png")


def case3_bayesian_inference():
    """案例3: 贝叶斯参数识别"""
    print("\n" + "="*70)
    print("案例3: 贝叶斯参数识别 (识别弹性模量)")
    print("="*70)
    
    # 真实参数
    E_true = 210e9  # Pa
    I_true = 8.33e-6  # m^4
    L = 2.0  # m
    P = 1000  # N
    q = 500  # N/m
    
    # 计算真实位移
    w_true = cantilever_beam_deflection(E_true, I_true, P, q, L)
    
    # 生成观测数据 (添加噪声)
    np.random.seed(42)
    noise_std = 0.001  # 1mm噪声
    n_obs = 10
    w_obs = w_true + np.random.randn(n_obs) * noise_std
    
    print(f"\n  真实参数:")
    print(f"    - 弹性模量 E: {E_true/1e9:.1f} GPa")
    print(f"    - 真实位移: {w_true*1000:.4f} mm")
    
    print(f"\n  观测数据:")
    print(f"    - 观测次数: {n_obs}")
    print(f"    - 观测位移均值: {np.mean(w_obs)*1000:.4f} mm")
    print(f"    - 观测位移标准差: {np.std(w_obs)*1000:.4f} mm")
    print(f"    - 噪声标准差: {noise_std*1000:.1f} mm")
    
    # 定义先验 (对数正态先验，因为E必须为正)
    E_prior_mean = 200e9
    E_prior_std = 50e9
    
    def log_prior(E):
        """对数先验 (正态分布)"""
        if E <= 0:
            return -np.inf
        return -0.5 * ((E - E_prior_mean) / E_prior_std)**2
    
    def log_likelihood(E, data):
        """对数似然"""
        w_model = cantilever_beam_deflection(E, I_true, P, q, L)
        n = len(data)
        log_lik = -0.5 * n * np.log(2 * np.pi * noise_std**2) - \
                  0.5 * np.sum((data - w_model)**2) / noise_std**2
        return log_lik
    
    # 创建贝叶斯推断对象
    bi = BayesianInference(log_prior, log_likelihood)
    
    # MCMC采样
    print("\n  MCMC采样...")
    E_init = [E_prior_mean]
    samples = bi.metropolis_hastings(w_obs, E_init, n_samples=20000, 
                                      proposal_std=5e9, burn_in=5000)
    
    # 后验统计
    stats = bi.get_posterior_stats()
    
    print(f"\n  MCMC结果:")
    print(f"    - 接受率: {bi.acceptance_rate*100:.1f}%")
    print(f"    - 后验均值: {stats['mean'][0]/1e9:.2f} GPa")
    print(f"    - 后验标准差: {stats['std'][0]/1e9:.2f} GPa")
    print(f"    - 95%置信区间: [{stats['ci_95'][0][0]/1e9:.2f}, {stats['ci_95'][1][0]/1e9:.2f}] GPa")
    
    # Geweke诊断
    z_scores = bi.geweke_diagnostic()
    print(f"    - Geweke Z-score: {z_scores[0]:.3f} (|Z| < 2 表示收敛)")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    
    # 1. 迹线图
    ax1 = axes[0, 0]
    ax1.plot(samples/1e9, linewidth=0.5, alpha=0.7)
    ax1.axhline(y=E_true/1e9, color='r', linestyle='--', linewidth=2, label=f'True E={E_true/1e9:.1f} GPa')
    ax1.axhline(y=stats['mean'][0]/1e9, color='g', linestyle='--', linewidth=2, label=f'Posterior Mean={stats["mean"][0]/1e9:.1f} GPa')
    ax1.set_xlabel('MCMC Iteration', fontsize=11)
    ax1.set_ylabel('E (GPa)', fontsize=11)
    ax1.set_title('MCMC Trace Plot', fontsize=12, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 2. 后验分布
    ax2 = axes[0, 1]
    ax2.hist(samples/1e9, bins=50, density=True, alpha=0.7, color='steelblue', edgecolor='black')
    ax2.axvline(x=E_true/1e9, color='r', linestyle='--', linewidth=2, label='True Value')
    ax2.axvline(x=stats['mean'][0]/1e9, color='g', linestyle='--', linewidth=2, label='Posterior Mean')
    ax2.axvline(x=stats['ci_95'][0][0]/1e9, color='orange', linestyle=':', linewidth=2, label='95% CI')
    ax2.axvline(x=stats['ci_95'][1][0]/1e9, color='orange', linestyle=':', linewidth=2)
    ax2.set_xlabel('E (GPa)', fontsize=11)
    ax2.set_ylabel('Posterior Density', fontsize=11)
    ax2.set_title('Posterior Distribution of E', fontsize=12, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 3. 自相关函数
    ax3 = axes[1, 0]
    max_lag = 100
    autocorr = np.correlate(samples[:, 0] - np.mean(samples[:, 0]), 
                            samples[:, 0] - np.mean(samples[:, 0]), mode='full')
    autocorr = autocorr[len(autocorr)//2:]
    autocorr = autocorr / autocorr[0]
    lags = np.arange(min(max_lag, len(autocorr)))
    ax3.plot(lags, autocorr[:len(lags)], 'b-', linewidth=1.5)
    ax3.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
    ax3.axhline(y=0.05, color='r', linestyle='--', linewidth=1, label='5% threshold')
    ax3.axhline(y=-0.05, color='r', linestyle='--', linewidth=1)
    ax3.set_xlabel('Lag', fontsize=11)
    ax3.set_ylabel('Autocorrelation', fontsize=11)
    ax3.set_title('Autocorrelation Function', fontsize=12, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 4. 先验vs后验对比
    ax4 = axes[1, 1]
    E_range = np.linspace(100e9, 300e9, 500)
    from scipy import stats as sp_stats
    prior_pdf = sp_stats.norm.pdf(E_range, E_prior_mean, E_prior_std)
    
    # 核密度估计后验
    from scipy.stats import gaussian_kde
    kde = gaussian_kde(samples[:, 0])
    posterior_pdf = kde(E_range)
    
    ax4.plot(E_range/1e9, prior_pdf/np.max(prior_pdf), 'b-', linewidth=2, label='Prior (normalized)')
    ax4.plot(E_range/1e9, posterior_pdf/np.max(posterior_pdf), 'r-', linewidth=2, label='Posterior (normalized)')
    ax4.axvline(x=E_true/1e9, color='g', linestyle='--', linewidth=2, label=f'True E={E_true/1e9:.1f} GPa')
    ax4.set_xlabel('E (GPa)', fontsize=11)
    ax4.set_ylabel('Normalized Density', fontsize=11)
    ax4.set_title('Prior vs Posterior', fontsize=12, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('d:\\文档\\强度仿真\\主题092\\贝叶斯参数识别.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("\n  可视化结果已保存: 贝叶斯参数识别.png")


def case4_active_subspace():
    """案例4: 主动子空间分析"""
    print("\n" + "="*70)
    print("案例4: 主动子空间分析")
    print("="*70)
    
    # 问题参数
    L = 2.0
    E_mean, E_std = 200e9, 10e9
    I_mean, I_std = 8.33e-6, 4e-7
    P_mean, P_std = 1000, 100
    q_mean, q_std = 500, 50
    
    def to_standard(x, mean, std):
        return (x - mean) / std
    
    def from_standard(xi, mean, std):
        return mean + xi * std
    
    def model(xi):
        E = from_standard(xi[0], E_mean, E_std)
        I = from_standard(xi[1], I_mean, I_std)
        P = from_standard(xi[2], P_mean, P_std)
        q = from_standard(xi[3], q_mean, q_std)
        return cantilever_beam_deflection(E, I, P, q, L)
    
    print("\n  主动子空间分析...")
    
    # 创建主动子空间分析对象
    as_analysis = ActiveSubspace(model)
    
    # 生成梯度样本
    n_samples = 100
    dim = 4
    X_samples = np.random.randn(n_samples, dim)
    
    # 计算梯度协方差矩阵
    eigenvalues, eigenvectors = as_analysis.compute_gradient_covariance(X_samples)
    
    print(f"\n  特征值分析:")
    var_names = ['E', 'I', 'P', 'q']
    print(f"    {'变量':<5} {'特征值':<15} {'能量占比':<12} {'累计占比':<12}")
    print(f"    {'-'*50}")
    cumsum = 0
    for i in range(dim):
        cumsum += eigenvalues[i]
        print(f"    {i+1:<5} {eigenvalues[i]:<15.6e} {eigenvalues[i]/np.sum(eigenvalues)*100:<12.2f}% {cumsum/np.sum(eigenvalues)*100:<12.2f}%")
    
    # 特征向量分析
    print(f"\n  特征向量 (主动方向):")
    print(f"    {'变量':<10}", end='')
    for i in range(dim):
        print(f"{'w'+str(i+1):<10}", end='')
    print()
    print(f"    {'-'*50}")
    for i, name in enumerate(var_names):
        print(f"    {name:<10}", end='')
        for j in range(dim):
            print(f"{eigenvectors[i, j]:<10.4f}", end='')
        print()
    
    # 确定主动子空间维度
    n_active = 2
    energy_ratio = as_analysis.energy_ratio(n_active)
    print(f"\n  主动子空间:")
    print(f"    - 维度: {n_active}")
    print(f"    - 能量占比: {energy_ratio*100:.2f}%")
    
    # 在主动子空间上构建响应面
    print("\n  构建降阶响应面...")
    
    # 生成主动变量样本
    n_train = 50
    y_train = np.random.randn(n_train, n_active)
    
    # 映射回全空间 (非主动方向设为0)
    X_train_full = np.zeros((n_train, dim))
    for i in range(n_train):
        X_train_full[i] = as_analysis.reconstruct_from_active(y_train[i], n_active)
    
    # 计算响应
    f_train = np.array([model(X_train_full[i]) for i in range(n_train)])
    
    # 使用GP构建响应面
    gp_as = GaussianProcessRegression(kernel_type='rbf')
    gp_as.fit(y_train, f_train)
    
    # 测试降阶模型
    n_test = 100
    y_test = np.random.randn(n_test, n_active)
    X_test_full = np.zeros((n_test, dim))
    for i in range(n_test):
        X_test_full[i] = as_analysis.reconstruct_from_active(y_test[i], n_active)
    
    f_true = np.array([model(X_test_full[i]) for i in range(n_test)])
    f_pred, _ = gp_as.predict(y_test, return_std=True)
    
    r2_as = 1 - np.sum((f_true - f_pred)**2) / np.sum((f_true - np.mean(f_true))**2)
    
    print(f"\n  降阶模型性能:")
    print(f"    - R²: {r2_as:.6f}")
    print(f"    - RMSE: {np.sqrt(np.mean((f_true - f_pred)**2))*1000:.4f} mm")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    
    # 1. 特征值谱
    ax1 = axes[0, 0]
    x_pos = np.arange(1, dim+1)
    ax1.semilogy(x_pos, eigenvalues, 'bo-', linewidth=2, markersize=8)
    ax1.axhline(y=eigenvalues[n_active-1], color='r', linestyle='--', 
                linewidth=2, label=f'Cutoff (n_active={n_active})')
    ax1.set_xlabel('Eigenvalue Index', fontsize=11)
    ax1.set_ylabel('Eigenvalue (log scale)', fontsize=11)
    ax1.set_title('Eigenvalue Spectrum', fontsize=12, fontweight='bold')
    ax1.set_xticks(x_pos)
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 2. 特征向量热图
    ax2 = axes[0, 1]
    im = ax2.imshow(eigenvectors, cmap='RdBu_r', aspect='auto', vmin=-1, vmax=1)
    ax2.set_xticks(range(dim))
    ax2.set_yticks(range(dim))
    ax2.set_xticklabels([f'w{i+1}' for i in range(dim)])
    ax2.set_yticklabels(var_names)
    ax2.set_xlabel('Eigenvector', fontsize=11)
    ax2.set_ylabel('Input Variable', fontsize=11)
    ax2.set_title('Eigenvector Components', fontsize=12, fontweight='bold')
    plt.colorbar(im, ax=ax2)
    
    # 3. 降阶响应面 (2D可视化)
    ax3 = axes[1, 0]
    if n_active >= 2:
        y1_range = np.linspace(-3, 3, 50)
        y2_range = np.linspace(-3, 3, 50)
        Y1, Y2 = np.meshgrid(y1_range, y2_range)
        
        F_pred = np.zeros_like(Y1)
        for i in range(len(y1_range)):
            for j in range(len(y2_range)):
                y = np.array([Y1[j, i], Y2[j, i]])
                pred_result = gp_as.predict(y.reshape(1, -1), return_std=True)
                F_pred[j, i] = pred_result[0][0] if isinstance(pred_result[0], np.ndarray) else pred_result[0]
        
        contour = ax3.contourf(Y1, Y2, F_pred*1000, levels=20, cmap='viridis')
        ax3.scatter(y_train[:, 0], y_train[:, 1], c='red', s=30, 
                   edgecolors='white', linewidth=0.5, label='Training Data')
        ax3.set_xlabel('Active Variable 1', fontsize=11)
        ax3.set_ylabel('Active Variable 2', fontsize=11)
        ax3.set_title('Response Surface in Active Subspace', fontsize=12, fontweight='bold')
        ax3.legend()
        plt.colorbar(contour, ax=ax3, label='Displacement (mm)')
    
    # 4. 降阶模型验证
    ax4 = axes[1, 1]
    ax4.scatter(f_true*1000, f_pred*1000, alpha=0.6, c='blue', s=30)
    ax4.plot([f_true.min()*1000, f_true.max()*1000], 
             [f_true.min()*1000, f_true.max()*1000], 'r--', lw=2, label='Perfect Prediction')
    ax4.set_xlabel('True Model (mm)', fontsize=11)
    ax4.set_ylabel('Reduced Model (mm)', fontsize=11)
    ax4.set_title(f'Reduced Model Validation (R²={r2_as:.4f})', fontsize=12, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('d:\\文档\\强度仿真\\主题092\\主动子空间分析.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("\n  可视化结果已保存: 主动子空间分析.png")


def case5_kl_expansion():
    """案例5: Karhunen-Loève展开随机场"""
    print("\n" + "="*70)
    print("案例5: Karhunen-Loève展开随机场")
    print("="*70)
    
    # 定义随机场参数
    L = 1.0  # 域长度
    sigma = 1.0  # 标准差
    lc = 0.1  # 相关长度
    
    # 定义均值函数和协方差核
    mean_func = lambda x: np.zeros_like(x)
    
    def cov_kernel(x, y):
        """指数型协方差核"""
        return sigma**2 * np.exp(-np.abs(x - y) / lc)
    
    print("\n  随机场参数:")
    print(f"    - 域长度: {L} m")
    print(f"    - 标准差: {sigma}")
    print(f"    - 相关长度: {lc} m")
    
    # 创建KL展开对象
    kle = KarhunenLoeveExpansion([0, L], mean_func, cov_kernel)
    
    # 离散化
    n_points = 100
    n_modes = 20
    eigenvalues, eigenfunctions = kle.discretize_nystrom(n_points, n_modes)
    
    print(f"\n  KL展开结果:")
    print(f"    - 离散点数: {n_points}")
    print(f"    - 保留模态数: {n_modes}")
    print(f"    - 前5个特征值: {eigenvalues[:5]}")
    
    # 计算能量占比
    energy_5 = np.sum(eigenvalues[:5]) / np.sum(eigenvalues)
    energy_10 = np.sum(eigenvalues[:10]) / np.sum(eigenvalues)
    energy_20 = np.sum(eigenvalues[:20]) / np.sum(eigenvalues)
    
    print(f"\n  能量占比:")
    print(f"    - 前5个模态: {energy_5*100:.2f}%")
    print(f"    - 前10个模态: {energy_10*100:.2f}%")
    print(f"    - 前20个模态: {energy_20*100:.2f}%")
    
    # 生成随机场样本
    n_samples = 5
    samples = kle.sample(n_samples)
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    
    # 1. 特征值谱
    ax1 = axes[0, 0]
    ax1.semilogy(range(1, len(eigenvalues)+1), eigenvalues, 'bo-', linewidth=2, markersize=6)
    ax1.axvline(x=5, color='r', linestyle='--', alpha=0.7, label='5 modes')
    ax1.axvline(x=10, color='g', linestyle='--', alpha=0.7, label='10 modes')
    ax1.set_xlabel('Mode Number', fontsize=11)
    ax1.set_ylabel('Eigenvalue (log scale)', fontsize=11)
    ax1.set_title('KL Eigenvalue Spectrum', fontsize=12, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 2. 特征函数
    ax2 = axes[0, 1]
    for i in range(min(5, n_modes)):
        ax2.plot(kle.x_grid, eigenfunctions[:, i], label=f'Mode {i+1}', linewidth=1.5)
    ax2.set_xlabel('Position x', fontsize=11)
    ax2.set_ylabel('Eigenfunction', fontsize=11)
    ax2.set_title('First 5 Eigenfunctions', fontsize=12, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 3. 随机场样本
    ax3 = axes[1, 0]
    for i in range(n_samples):
        ax3.plot(kle.x_grid, samples[i], alpha=0.7, linewidth=1.5)
    ax3.set_xlabel('Position x', fontsize=11)
    ax3.set_ylabel('Field Value', fontsize=11)
    ax3.set_title('Random Field Samples', fontsize=12, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    
    # 4. 协方差核重构
    ax4 = axes[1, 1]
    x_grid = kle.x_grid
    C_reconstructed = np.zeros((len(x_grid), len(x_grid)))
    for i in range(n_modes):
        C_reconstructed += eigenvalues[i] * np.outer(eigenfunctions[:, i], eigenfunctions[:, i])
    
    # 理论协方差
    C_theory = np.zeros((len(x_grid), len(x_grid)))
    for i, xi in enumerate(x_grid):
        for j, xj in enumerate(x_grid):
            C_theory[i, j] = cov_kernel(xi, xj)
    
    im = ax4.imshow(C_reconstructed, extent=[0, L, 0, L], origin='lower', cmap='viridis')
    ax4.set_xlabel('Position x', fontsize=11)
    ax4.set_ylabel('Position y', fontsize=11)
    ax4.set_title('Reconstructed Covariance Kernel', fontsize=12, fontweight='bold')
    plt.colorbar(im, ax=ax4)
    
    plt.tight_layout()
    plt.savefig('d:\\文档\\强度仿真\\主题092\\KL随机场展开.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("\n  可视化结果已保存: KL随机场展开.png")


# ==================== 7. 主程序 ====================

if __name__ == "__main__":
    print("\n" + "="*70)
    print("主题092：不确定性量化的先进方法")
    print("Python仿真程序")
    print("="*70)
    
    # 设置随机种子
    np.random.seed(42)
    
    # 运行所有案例
    case1_pce_analysis()
    case2_gaussian_process()
    case3_bayesian_inference()
    case4_active_subspace()
    case5_kl_expansion()
    
    print("\n" + "="*70)
    print("  所有仿真案例已完成！")
    print("="*70)
    print("\n生成的可视化文件:")
    print("  1. PCE不确定性分析.png")
    print("  2. GP代理模型.png")
    print("  3. 贝叶斯参数识别.png")
    print("  4. 主动子空间分析.png")
    print("  5. KL随机场展开.png")
    print("="*70)