无需拟合模型的差分熵差异贡献特征识别及P值计算问询

核心思路 #

直接针对差分熵差异这个统计量，通过逐个移除单观测后重新计算差异，用差异变化量衡量观测的贡献度；再通过置换检验计算贡献的显著性P值，全程无需拟合任何预测模型。

关键步骤 #

1. 计算原始两组样本的差分熵差异 ( D = |H_{\text{group1}} - H_{\text{group2}}| )
2. 对每个观测（group1/group2各样本），移除后重新计算差分熵差异 ( D' )，得到变化量 ( \Delta = |D - D'| )，( \Delta ) 越大说明该观测对原始差异的贡献越大
3. 用置换检验验证贡献的显著性：混合两组样本后随机重分组，重复步骤1-2得到置换后的 ( \Delta ) 分布，原始 ( \Delta ) 在分布中的分位数即为P值

Python 实现示例 #

import numpy as np
from scipy.stats import gaussian_kde

# 定义差分熵计算函数（基于核密度估计）
def differential_entropy(x):
    kde = gaussian_kde(x)
    xs = np.linspace(x.min() - 0.1*(x.max()-x.min()), x.max() + 0.1*(x.max()-x.min()), 1000)
    pdf = kde(xs)
    # 避免log(0)，替换极小值为1e-10
    pdf = np.where(pdf < 1e-10, 1e-10, pdf)
    return -np.trapz(pdf * np.log(pdf), xs)

# 模拟样本（替换为你的真实数据）
np.random.seed(42)
group1 = np.random.normal(loc=0, scale=1, size=100)
group2 = np.random.normal(loc=1, scale=1, size=100)

# 计算原始差分熵差异
h1 = differential_entropy(group1)
h2 = differential_entropy(group2)
original_diff = abs(h1 - h2)

# 计算每个观测的贡献度delta
delta_list = []
# 处理group1的每个观测
for i in range(len(group1)):
    new_group1 = np.delete(group1, i)
    new_h1 = differential_entropy(new_group1)
    new_diff = abs(new_h1 - h2)
    delta_list.append(("group1", i, abs(original_diff - new_diff)))

# 处理group2的每个观测
for i in range(len(group2)):
    new_group2 = np.delete(group2, i)
    new_h2 = differential_entropy(new_group2)
    new_diff = abs(h1 - new_h2)
    delta_list.append(("group2", i, abs(original_diff - new_diff)))

# 按贡献度降序排序
delta_list.sort(key=lambda x: x[2], reverse=True)
print("Top 5贡献最大的观测：")
for item in delta_list[:5]:
    print(f"组：{item[0]}, 索引：{item[1]}, 差异变化量：{item[2]:.4f}")

# 置换检验计算P值（以top1观测为例）
top_group, top_idx, top_delta = delta_list[0]
perm_deltas = []
n_perm = 1000  # 置换次数，可根据需求调整
combined = np.concatenate([group1, group2])

for _ in range(n_perm):
    # 随机重分组
    perm_idx = np.random.permutation(len(combined))
    perm_g1 = combined[perm_idx[:len(group1)]]
    perm_g2 = combined[perm_idx[len(group1):]]
    # 计算原始差异
    perm_h1 = differential_entropy(perm_g1)
    perm_h2 = differential_entropy(perm_g2)
    perm_original_diff = abs(perm_h1 - perm_h2)
    # 移除对应组的一个观测（模拟top观测的位置逻辑）
    if top_group == "group1":
        new_perm_g1 = np.delete(perm_g1, 0)
        new_perm_h1 = differential_entropy(new_perm_g1)
        perm_new_diff = abs(new_perm_h1 - perm_h2)
    else:
        new_perm_g2 = np.delete(perm_g2, 0)
        new_perm_h2 = differential_entropy(new_perm_g2)
        perm_new_diff = abs(perm_h1 - new_perm_h2)
    perm_deltas.append(abs(perm_original_diff - perm_new_diff))

# 计算P值：原始delta大于等于置换delta的比例
p_value = np.mean(np.array(perm_deltas) >= top_delta)
print(f"Top1观测的贡献P值：{p_value:.4f}")

R 实现示例 #

# 定义差分熵计算函数
differential_entropy <- function(x) {
  kde <- density(x, n = 1000)
  pdf <- kde$y
  pdf[pdf < 1e-10] <- 1e-10  # 避免log(0)
  -sum(pdf * log(pdf)) * diff(kde$x)[1]
}

# 模拟样本（替换为真实数据）
set.seed(42)
group1 <- rnorm(100, mean = 0, sd = 1)
group2 <- rnorm(100, mean = 1, sd = 1)

# 原始差分熵差异
h1 <- differential_entropy(group1)
h2 <- differential_entropy(group2)
original_diff <- abs(h1 - h2)

# 计算每个观测的delta
delta_list <- list()

# 处理group1
for (i in 1:length(group1)) {
  new_group1 <- group1[-i]
  new_h1 <- differential_entropy(new_group1)
  new_diff <- abs(new_h1 - h2)
  delta_list[[length(delta_list)+1]] <- list(group = "group1", idx = i, delta = abs(original_diff - new_diff))
}

# 处理group2
for (i in 1:length(group2)) {
  new_group2 <- group2[-i]
  new_h2 <- differential_entropy(new_group2)
  new_diff <- abs(h1 - new_h2)
  delta_list[[length(delta_list)+1]] <- list(group = "group2", idx = i, delta = abs(original_diff - new_diff))
}

# 排序输出
delta_df <- do.call(rbind, lapply(delta_list, function(x) data.frame(x)))
delta_df <- delta_df[order(-delta_df$delta), ]
print("Top 5贡献最大的观测：")
print(head(delta_df, 5))

# 置换检验计算P值（以top1为例）
top_row <- delta_df[1, ]
n_perm <- 1000
perm_deltas <- numeric(n_perm)
combined <- c(group1, group2)

for (k in 1:n_perm) {
  perm_idx <- sample(length(combined))
  perm_g1 <- combined[perm_idx[1:length(group1)]]
  perm_g2 <- combined[perm_idx[(length(group1)+1):length(combined)]]
  perm_h1 <- differential_entropy(perm_g1)
  perm_h2 <- differential_entropy(perm_g2)
  perm_original_diff <- abs(perm_h1 - perm_h2)
  
  if (top_row$group == "group1") {
    new_perm_g1 <- perm_g1[-1]
    new_perm_h1 <- differential_entropy(new_perm_g1)
    perm_new_diff <- abs(new_perm_h1 - perm_h2)
  } else {
    new_perm_g2 <- perm_g2[-1]
    new_perm_h2 <- differential_entropy(new_perm_g2)
    perm_new_diff <- abs(perm_h1 - new_perm_h2)
  }
  perm_deltas[k] <- abs(perm_original_diff - perm_new_diff)
}

p_value <- mean(perm_deltas >= top_row$delta)
cat(sprintf("Top1观测的贡献P值：%.4f\n", p_value))

优化效率提示 #

• 对于大样本，可采用并行计算加速置换检验（Python用multiprocessing，R用parallel包）
• 若无需精确P值，可减少置换次数（比如从1000降到200）
• 差分熵计算时，可固定核密度估计的带宽，避免重复计算参数

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