求助：基于KDE Bootstrap求解分布上下限及预测区间的代码问题

问题背景 #

目标是通过以下步骤构建预测区间：

• 计算实际值与预测值的残差
• 归一化残差
• 生成残差的KDE
• 对KDE重采样
• 求解分布上下限
• 构建预测区间

现有代码 #

# In[152]
sample = kde.resample(43826)
5 ** np.percentile(sample, 5)

# In[153]
def resample_kde_percentile(kde):
    sample = kde.resample(kde.n)
    return 5 ** np.percentile(sample, 5)

# In[154]
def summarize(t, digits=2):
     table = pd.DataFrame(columns=['Estimate', 'SE', 'CI90'])
     est = np.mean(t).round(digits)
     SE = np.std(t).round(digits)
     CI95 = np.percentile(t, [5, 95]).round(digits)
     table.loc[''] = est, SE, CI95
     return table

# In[155]
t10 = [resample_kde_percentile(kde)
   for i in range(1000)]

summary10 = summarize(t10)
summary10

代码问题分析 #

1. 仅计算单侧分位数：resample_kde_percentile只返回了5%分位数的转换值，未计算95%分位数，无法得到完整的分布上下限
2. 列名与计算逻辑不一致：summarize函数列名写的是CI90，但实际计算的是90%置信区间对应的[5,95]分位数，命名混淆
3. 缺失预测区间核心逻辑：当前代码仅处理了残差的分位数，没有将重采样的残差与原始预测值结合，无法生成最终的预测区间

修正后的完整代码 #

1. 完整的残差重采样与分位数计算函数 #

def resample_kde_ci(kde):
    # 从KDE重采样残差
    resampled_residuals = kde.resample(kde.n)[0].flatten()
    # 计算5%和95%分位数，再做反归一化转换（这里假设归一化时用了log5，所以用5**还原）
    lower = 5 ** np.percentile(resampled_residuals, 5)
    upper = 5 ** np.percentile(resampled_residuals, 95)
    return lower, upper

2. 修正后的结果汇总函数 #

def summarize(t, digits=2):
    # 拆分上下限数组
    lowers = [x[0] for x in t]
    uppers = [x[1] for x in t]
    # 计算统计量
    est_lower = np.mean(lowers).round(digits)
    se_lower = np.std(lowers).round(digits)
    ci_lower = np.percentile(lowers, [5, 95]).round(digits)
    
    est_upper = np.mean(uppers).round(digits)
    se_upper = np.std(uppers).round(digits)
    ci_upper = np.percentile(uppers, [5, 95]).round(digits)
    
    # 构建汇总表
    table = pd.DataFrame(
        index=['Lower Bound', 'Upper Bound'],
        columns=['Estimate', 'SE', 'CI95']
    )
    table.loc['Lower Bound'] = [est_lower, se_lower, ci_lower]
    table.loc['Upper Bound'] = [est_upper, se_upper, ci_upper]
    return table

3. 执行自助法并构建预测区间 #

# 执行1000次自助重采样
bootstrap_results = [resample_kde_ci(kde) for _ in range(1000)]

# 生成汇总表
summary = summarize(bootstrap_results)
print(summary)

# 假设你有原始的预测值数组y_pred，生成最终预测区间
# 注意：如果残差是实际值-预测值，那么预测区间应该是 y_pred ± 残差分位数
# 这里根据你的归一化逻辑调整，示例如下：
# y_pred = ... # 你的预测值数组
# lower_pred = y_pred - np.mean([x[0] for x in bootstrap_results])
# upper_pred = y_pred + np.mean([x[1] for x in bootstrap_results])

关键补充说明 #

• 确保残差的归一化/反归一化逻辑一致：如果残差是做了log5转换后再拟合KDE，那么反转换必须用5**，否则要对应调整
• 预测区间的构建需要结合原始预测值：残差的上下限是用来调整预测值的，最终预测区间应为预测值 ± 残差的分位数（或自助法得到的平均上下限）
• 自助法的重采样次数建议不少于1000次，保证结果稳定性

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