如何在Python中通过解析法计算多项式混沌分解的留一法交叉验证Q2分数

要在OpenTURNS中无需重新拟合模型，用解析法计算多项式混沌展开（PCE）的留一法交叉验证（LOOCV）Q²分数，我们可以利用最小二乘模型的LOOCV解析性质，或者直接使用OpenTURNS实验模块的内置工具。下面分两种方法详细说明：

方法一：使用OpenTURNS Experimental模块的内置LOOCV工具（推荐） #

OpenTURNS的openturns.experimental模块专门提供了FunctionalChaosLOOCV类，它封装了解析法计算LOOCV指标的逻辑，无需手动推导公式，直接基于已有的PCE结果就能快速得到Q²分数。

在你原代码的基础上，只需添加以下几行：

# 从PCE结果创建LOOCV分析对象
loocv_analyzer = otexp.FunctionalChaosLOOCV(result)
# 计算LOOCV的Q²分数
q2_score = loocv_analyzer.computeQ2()
print(f"LOOCV Q²分数（内置工具）: {q2_score:.4f}")

这个方法的优势是：

1. 无需重复计算设计矩阵、逆矩阵等中间量，效率更高
2. 自动处理带惩罚的最小二乘场景（如你使用的PenalizedLeastSquaresAlgorithmFactory）
3. 避免手动计算可能引入的数值错误

方法二：手动实现解析法（理解原理） #

如果想深入理解底层逻辑，可以手动基于最小二乘模型的LOOCV解析公式实现计算。

核心原理 #

对于第i个训练样本：

1. 原模型的残差：$e_i = y_i - \hat{y}_i$（$\hat{y}_i$是原模型对第i个样本的预测值）
2. 帽子矩阵对角线元素$h_{ii}$：表示第i个样本对自身预测结果的影响程度，对应矩阵$H = \Phi (\Phi^T \Phi + \lambda I)^{-1}\PhiT$的第i个对角线元素（$\lambda$是L2惩罚系数）
3. LOOCV的预测残差：$y_i - \hat{y}{(-i)} = \frac{e_i}{1 - h{ii}}$（$\hat{y}_{(-i)}$是删除第i个样本后模型对第i个样本的预测值）
4. Q²分数的计算公式：
   $Q^2 = 1 - \frac{\sum_{i=1}^n (y_i - \hat{y}_{(-i)})^{2}{\sum_{i=1}}n (y_i - \bar{y})^2}$
   其中$\bar{y}$是输出样本的均值。

手动计算代码实现 #

在你原代码后添加以下代码：

# --------------------------
# 步骤1：计算原模型的预测值和残差
# --------------------------
pce_metamodel = result.getMetaModel()
y_hat = pce_metamodel(inputTrain)
residuals = outputTrain - y_hat

# --------------------------
# 步骤2：计算设计矩阵与帽子矩阵对角线
# --------------------------
# 获取PCE的基函数集合
total_degree = 8
enumerate_func = multivariateBasis.getEnumerateFunction()
basis_indices = enumerate_func.getIndicesFromTotalDegree(total_degree)
pce_basis = multivariateBasis.buildBasis(basis_indices)
n_basis = len(pce_basis)
n_samples = inputTrain.getSize()

# 构建设计矩阵Phi：每行对应一个样本在基函数上的取值
Phi = ot.Matrix(n_samples, n_basis)
for i in range(n_samples):
    input_i = inputTrain[i]
    for j in range(n_basis):
        Phi[i, j] = pce_basis[j](input_i)

# 获取L2惩罚的正则化系数（针对PenalizedLeastSquares）
pls_algorithm = projectionStrategy.getAlgorithm()
lambda_reg = pls_algorithm.getLambda()

# 计算正则化后的(Phi^T Phi + lambda*I)及其逆
PhiT_Phi = Phi.computeGram()
regularized_gram = PhiT_Phi + lambda_reg * ot.IdentityMatrix(n_basis)
inv_regularized_gram = regularized_gram.computeInverse()

# 计算帽子矩阵的对角线元素h_ii
h_diag = ot.Sample(n_samples, 1)
for i in range(n_samples):
    phi_row = Phi.getRow(i)
    h_diag[i, 0] = phi_row.dot(inv_regularized_gram).dot(phi_row)

# --------------------------
# 步骤3：计算LOOCV残差与Q²分数
# --------------------------
# 计算LOOCV残差
loocv_residuals = residuals / (1.0 - h_diag)

# 计算总平方和SST = sum((y_i - y_mean)^2)
y_mean = outputTrain.computeMean()[0]
sst = ((outputTrain - y_mean) ** 2).computeSum()[0]

# 计算LOOCV残差平方和SSE_LOOCV
sse_loocv = (loocv_residuals ** 2).computeSum()[0]

# 计算Q²
q2_manual = 1.0 - sse_loocv / sst
print(f"LOOCV Q²分数（手动计算）: {q2_manual:.4f}")

关键注意事项 #

1. 惩罚项的处理：如果使用了带惩罚的最小二乘（如示例中的PenalizedLeastSquaresAlgorithmFactory），手动计算时必须包含正则化项，否则帽子矩阵的计算会出错，导致Q²结果偏差。
2. 数值稳定性：当样本量远大于基函数数量时，$h_{ii}$通常远小于1，不会出现分母为0的情况；如果基函数数量接近或超过样本量，建议增加样本数或降低PCE的总阶数。
3. 效率对比：内置工具的计算效率更高，因为它直接利用了PCE拟合过程中的中间结果，避免了重复计算设计矩阵和逆矩阵，适合大规模样本的场景。

备注：内容来源于stack exchange，提问作者Michael Baudin