Python环境下Keras LSTM处理3维时序数据的模型无关解释实现方法

核心适配逻辑 #

• 常规PDP工具默认适配2维输入（样本数, 特征数），无法直接兼容LSTM的3维时序结构（样本数, 时间步长, 特征数），核心改造思路是围绕LSTM的输入规则调整PDP的扰动逻辑，不需要额外引入特殊框架。
• PDP的核心计算逻辑是：固定目标特征的候选取值，让模型对其余特征保持原始分布的样本做批量预测，取所有样本预测值的平均值，得到目标特征和预测输出的全局边际关系。
• 针对不等长时间步的适配要点：计算过程中仅对每个样本实际有效的时间步做特征替换，masked的填充无效位置不做修改，避免引入计算噪声。
• 不依赖第三方解释库的内置接口：这类接口大多默认假设输入是2维结构化数据，直接传入3维时序数组会触发维度不匹配的报错；基于PDP核心逻辑自行实现的流程，完全匹配时序模型的输入要求，计算结果符合全局可解释性的分析目标。

参考实现代码 #

import numpy as np
import matplotlib.pyplot as plt
from tensorflow import keras

# 基础参数配置
# 加载已训练完成的模型与数据集，形状为(样本ID数, 最大时间步长, 输入特征数)
X_train = np.load("你的训练集存储路径.npy")  # 对应形状(18610, 151, 18)
trained_lstm = keras.models.load_model("你的LSTM模型存储路径.h5")
n_samples, n_timesteps, n_features = X_train.shape
feature_names = [f"输入变量{i+1}" for i in range(n_features)]  # 替换为你的实际特征名
n_grid_points = 50  # 每个特征取50个候选网格值计算PDP

# 识别有效时间步掩码，适配不等长时间步场景
valid_timestep_mask = ~np.all(X_train == 0, axis=-1)

# 计算每个特征的有效取值范围，取1%和99%分位避开极端值干扰
feature_range_dict = {}
for feat_idx in range(n_features):
    all_feat_vals = X_train[:, :, feat_idx]
    valid_feat_vals = all_feat_vals[valid_timestep_mask]
    feat_min = np.percentile(valid_feat_vals, 1)
    feat_max = np.percentile(valid_feat_vals, 99)
    feature_range_dict[feat_idx] = np.linspace(feat_min, feat_max, n_grid_points)

# 逐特征计算PDP取值
pdp_result_dict = {}
for feat_idx in range(n_features):
    grid_vals = feature_range_dict[feat_idx]
    mean_pred_list = []
    # 遍历当前特征的所有候选网格值
    for grid_val in grid_vals:
        X_modified = X_train.copy()
        # 仅对有效时间步的目标特征做值替换，masked位置保持原始状态
        X_modified[valid_timestep_mask, feat_idx] = grid_val
        # 批量推理得到所有样本预测为类别1的概率
        batch_preds = trained_lstm.predict(X_modified, verbose=0)
        # 取所有样本预测概率的平均值作为当前网格值的PDP结果
        mean_pred_list.append(np.mean(batch_preds[:, 1]))
    pdp_result_dict[feat_idx] = (grid_vals, np.array(mean_pred_list))

# 绘制所有特征的PDP图
fig, axes = plt.subplots(6, 3, figsize=(15, 20))
axes = axes.flatten()
for feat_idx in range(n_features):
    x_grid, y_pdp = pdp_result_dict[feat_idx]
    axes[feat_idx].plot(x_grid, y_pdp)
    axes[feat_idx].set_title(feature_names[feat_idx])
    axes[feat_idx].set_xlabel("特征取值")
    axes[feat_idx].set_ylabel("预测为类别1的平均概率")
plt.tight_layout()
plt.show()

结果解读规则 #

• 影响强度量化：PDP曲线的斜率绝对值越高，代表该特征取值变化时模型预测概率的波动越大，对预测结果的影响越强。
• 影响方向量化：若PDP曲线随特征取值升高呈上升趋势，说明该特征取值越高，模型预测为类别1的概率越高，为正向影响；反之则为负向影响。

其他可适配的模型无关方法 #

• Shapley值方法：复用上述的3维特征替换逻辑，将PDP中对单特征的固定值替换，调整为Shapley值计算过程中的特征联盟扰动逻辑，即可在不改动模型的前提下，计算所有样本的全局Shapley值，量化每个特征的平均贡献度。
• LIME方法：针对单样本的3维时序序列，对目标特征在所有有效时间步做符合原始分布的随机扰动，生成局部解释的训练样本集，即可训练得到单样本的特征贡献权重，实现局部可解释性。

注意事项 #

• 计算过程中禁止对masked的填充无效时间步做特征替换，这类位置对LSTM的预测输出无贡献，替换后会引入计算偏差。
• 若特征之间存在强相关性，可将PDP的固定值替换逻辑调整为联合分布采样，避免生成实际场景中不存在的特征组合，提升结果的合理性。
• 整个计算流程共需执行18个特征×50个网格点=900次批量推理，基于常规GPU的推理速度，全流程计算时长约30-60分钟，无性能瓶颈。

内容的提问来源于stack exchange，提问作者Danique