简介：面向机器学习与特征工程场景的 K 近邻互信息计算程序，使用 MATLAB 实现，将 KNN 分类算法与互信息测量结合在一起，帮助数据科学与机器学习从业者评估特征关联性、筛选高信息量特征，从而提升 KNN 模型性能。资源共 2 个文件，包括 1 个 m 脚本与 1 个 txt 说明文件，压缩包整体约 2KB，结构精简，适合快速理解算法核心。目前已有 332 人学习下载。程序围绕数据预处理、距离计算、K 值选择、互信息计算、特征选择及 KNN 应用等环节组织，其中包含数据导入、清洗与标准化流程，还可通过交叉验证确定最优 K 值，比较不同特征子集对分类精度的影响。既便于动手复现 KraskovMI 估计思路，也方便在此基础上扩展不同距离度量或改进计算效率。对正在学习近邻算法和信息论特征选择的读者，这是一份轻量、可直接运行且易于调试的参考示例。
1. 互信息不是相关系数，K 近邻是它的无参估计器
做特征筛选时，很多人习惯先看 Pearson 相关系数，但线性关系之外的关联它完全看不见。两个变量呈正弦、分段或“X 型”结构时，相关系数可能趋近于零，互信息（Mutual Information）却能稳定地给出非零度量。互信息不假设概率分布形式，也不限制变量维度，因此被广泛用在 KNN 分类前的特征选择里。问题在于，真实数据集的联合概率密度未知，互信息的积分只能靠估计。K 近邻方法绕开了直方图分箱，直接用每个样本到第 k 个近邻的距离来近似局部密度，这就是 Kraskov-Stögbauer-Grassberger（KSG）估计器的基本思想。本资源中的KraskovMI.m正是这一思想的 MATLAB 实现，适合做特征相关性分析、聚类前变量筛选和 KNN 模型解释，尤其适合处理连续型特征和分类标签共存的场景。
2. KraskovMI.m 的算法拆解：从双重 k 近邻到互信息估计
2.1 KSG 估计为什么比传统分箱法更稳
传统的互信息计算先把变量取值范围切成网格，再统计每个格子里的样本频数。网格宽度设得太大，密度细节被抹平；设得太小，大量格子为空，估计方差飙升。KSG 方法不再依赖固定网格，而是对每个样本点寻找其在联合空间中的 k 个近邻，用这些近邻的距离反推局部概率密度。这个做法在样本量不大时依然能保持较低偏差，在连续型变量上尤其明显。
KSG 有两个常用形式。KSG1 使用最大范数（切比雪夫距离）衡量联合空间中的邻域，然后分别统计两个单变量空间里落在该邻域内的点数。KSG2 则用欧氏距离并引入高阶修正项，对边界效应更敏感。KraskovMI.m通常默认走 KSG1 路线，因为它实现简单、对 MATLAB 的knnsearch支持最友好，计算代价也低。实际比较中，KSG1 在特征维度低、样本量在几百到几千时表现已经足够好；KSG2 在变量强相关且样本量过万时稍优，但耗时明显增加。
2.2 KraskovMI.m 的核心实现结构
一个典型的KraskovMI.m函数会做四件事：读取两个输入向量或两个特征矩阵、用knnsearch找联合空间的 k 近邻、统计两个边际空间内的近邻数、用 digamma 函数做偏差校正。伪代码形式如下：
function mi = KraskovMI(x, y, k, varargin) % x, y: N x 1 列向量 % k: 近邻数量，默认 3 % varargin: 可传入距离度量 if nargin < 3 || isempty(k) k = 3; end N = length(x); xy = [x(:), y(:)]; % 1. 对每个点找联合空间中的第 k 近邻距离 [idxJK, distK] = knnsearch(xy, xy, 'K', k+1); distK = distK(:, end); % 第 k 近邻距离 % 2. 统计单变量空间内的点数 nx, ny nx = zeros(N,1); ny = zeros(N,1); for i = 1:N nx(i) = sum(abs(x - x(i)) <= distK(i)) - 1; ny(i) = sum(abs(y - y(i)) <= distK(i)) - 1; end % 3. 计算互信息 mi = psi(k) + psi(N) - mean(psi(nx + 1) + psi(ny + 1)); end逻辑说明：第一步里knnsearch(xy, xy, 'K', k+1)返回每个点按距离升序排列的邻居索引和距离，去掉自身后取第 k 个距离作为联合邻域半径。第二步对每个点分别计算x和y方向落在该半径内的点数。第三步中psi是 digamma 函数，用来代替对数的期望，减小有限样本偏差。参数k控制邻域大小，过小会使估计方差增大，过大则会把局部结构平滑掉。
2.3 距离度量和 k 值对估计值的影响
knnsearch默认使用欧氏距离，但在 KSG1 中应改用切比雪夫距离，因为联合空间的矩形邻域在边际投影后才是对称区间。若使用欧氏距离，nx和ny的统计会出现不对称，导致互信息被系统性低估。实际调试时建议在调用函数前先固定距离度量：
[idxDist, distK] = knnsearch(xy, xy, 'K', k+1, 'Distance', 'chebychev');这里'chebychev'对应最大范数，和 KSG 原始论文设定一致。修改后nx和ny的计数才具有可比性。k 值的选择可以参考floor(sqrt(N/2))作为起点，再向两侧各取 2 到 5 个值观察估计值是否稳定。如果 MI 结果随 k 剧烈波动，说明数据中存在离群点或分布过于尖锐，需要先做标准化或裁剪。
3. MATLAB 环境里的数据准备与 k 值选择
3.1 输入数据的格式与预处理
KraskovMI.m期望的输入是两列长度相同的数值向量。如果特征是类别型（如性别、地区），需要先做数值编码，但编码后的间距没有实际含义，此时互信息的结果只能作为粗糙的参考；更稳妥的做法是分别计算该特征与目标变量的条件概率表，再用离散互信息公式验证。连续特征则需要注意量纲差异对近邻搜索的干扰。例如一个特征取值在 0～1，另一个在 0～10000，后者会直接主导距离计算，使前者的局部密度形同虚设。
常见做法是先做 Z-score 标准化，再进入互信息计算：
data = readmatrix('features.csv'); X = data(:, 1:end-1); Y = data(:, end); X = (X - mean(X)) ./ std(X); Y = (Y - mean(Y)) ./ std(Y);标准化后每个维度的尺度一致，knnsearch的距离计算才不会被某一列支配。执行上述操作后，建议用any(isnan(X), 'all')检查缺失值，缺失样本直接删除，不要用均值填充，因为填充值会产生虚假的紧密邻域，导致互信息偏高。数据量在 10 万行以上时，knnsearch的最近邻搜索会占用大量内存，此时可以先用datasample随机抽样 1 万到 2 万行做初次估计。
3.2 k 值选择的经验法则与交叉验证
k 值直接决定局部邻域的半径。k 太小，每个点的邻域只覆盖到最近的一两个点，估计结果对噪声敏感；k 太大，邻域跨越多个密度区域，互信息倾向于被平滑到接近零。经验法则分两种：按样本量设定，常用的取值范围是 1 到 10，或用k = floor(sqrt(N))做粗调；按使用场景设定，特征选择阶段用较小 k 值，侧重保留强相关特征，建模阶段用大一点的 k 值，减少方差。
下面是一段对 k 做网格扫描的代码：
kList = [1, 2, 3, 5, 8, 13]; miHistory = zeros(length(kList), size(X,2)); for j = 1:size(X,2) for t = 1:length(kList) miHistory(t,j) = KraskovMI(X(:,j), Y, kList(t)); end end plot(kList, miHistory, 'o-');运行后把趋势线画出来，观察所有特征在 k 增大时是否收敛到相近的排序。如果某条曲线上下跳动超过 0.3，说明该特征分布带有过多离群点，需要检查是否混入了错误采集的样本。下表是不同数据规模下的经验 k 值参考：
| 样本量 N | 推荐 k | 备注 |
|---|---|---|
| 50～200 | 1～3 | 样本少，只能捕捉强相关 |
| 200～1000 | 3～5 | 推荐 k=3 起步 |
| 1000～10000 | 5～8 | 可同时比较 k=5 和 k=8 |
| >10000 | 8～15 | 注意计算时间，可先抽样 |
3.3 与 MATLAB 内置函数协同工作
KraskovMI.m依赖的核心函数是knnsearch，它来自 Statistics and Machine Learning Toolbox。如果环境中没有该工具箱，可以用暴力计算距离矩阵的方式替代：distMat = sqrt((x-x').^2 + (y-y').^2)，然后对每行排序取第 k 个值。距离矩阵法在 N 小于 5000 时速度尚可，超过后内存占用呈平方增长，不建议使用。MATLAB 2022b 之后knnsearch默认采用 KD 树策略，在高维特征下会自动退化为穷举搜索，因此不用刻意设置'NSMethod'参数，除非特征维度超过 32 且样本量极大。
4. 用互信息做特征选择：一个可复现的 MATLAB 流程
4.1 计算所有特征与目标变量的互信息
这里的目标变量可以是分类标签，也可以是连续值。KraskovMI.m对两类目标都适用，因为 KSG 估计本身不关心变量的类型，只看联合分布。计算前需要对分类标签做编码，常见做法是标量编码成一个整数列。获取全部特征重要性排名的代码如下：
numFeatures = size(X, 2); miScores = zeros(numFeatures, 1); kFixed = 5; for f = 1:numFeatures miScores(f) = KraskovMI(X(:, f), Y, kFixed); end [B, idx] = sort(miScores, 'descend'); fprintf('Top 5 features: %s\n', mat2str(idx(1:5)));逻辑说明：for循环把每个特征单独取出来与目标列组成两列矩阵，传入KraskovMI计算互信息。sort返回降序排名和对应的索引。输出结果后，不要只看分数最高的特征，还要对比分数最低的几个特征，若两者差距很小，说明该数据集中互信息的区分度有限，需要补充组合特征或改用条件互信息。参数kFixed设为 5 是折中值；若特征数量超过 50，可先把 k 调到 3，加快迭代速度。
4.2 对比互信息排序与 KNN 分类精度
特征选择的目标不只是找到相关特征，还要验证下游 KNN 模型的收益。常见做法是取排名前 m 的特征子集，训练 KNN 分类器，并与使用全部特征时的精度对比。下面脚本实现这一流程：
rng(42); numTop = [1, 3, 5, 10, size(X,2)]; accuracy = zeros(length(numTop), 1); for t = 1:length(numTop) topM = idx(1:numTop(t)); Xsub = X(:, topM); mdl = fitcknn(Xsub, Y, 'NumNeighbors', 7, 'Standardize', true); cvmdl = crossval(mdl, 'KFold', 5); accuracy(t) = 1 - kfoldLoss(cvmdl); end这里的fitcknn使用欧氏距离和多数投票，Standardize设为true会再次标准化，但不会改变此前互信息排序结果，因为标准化是单调变换，不会影响近邻的秩。kfoldLoss返回交叉验证损失，1 - loss即为精度。运行后观察精度曲线，如果只取前 3 个特征时精度已经接近全特征，说明其余特征是冗余的。若精度不升反降，问题通常出在 k 值选择上，应返回第 4.1 节用更小的 k 重新计算互信息。
4.3 结果分析：互信息与 KNN 精度的联动关系
下表为一个模拟数据集的输出示例，用于理解结果形态：
| 特征子集大小 | 平均互信息 | KNN 交叉验证精度 |
|---|---|---|
| 1 | 0.42 | 0.781 |
| 3 | 0.38 | 0.839 |
| 5 | 0.31 | 0.855 |
| 10 | 0.25 | 0.861 |
| 全特征(15) | 0.19 | 0.834 |
从表中可以看到，加入第 4 到第 10 个特征时互信息平均值持续下降，但 KNN 精度仍在缓慢上升，这说明这些特征单独看来相关性弱，却提供了互补的判别信息。全部特征加入后精度回落，则说明高维噪声开始干扰近邻计算。遇到这种情况，建议改用递归特征消除或条件互信息，而不是简单按单变量排序截断。
5. 偏差修正与计算加速：让 KraskovMI.m 更可靠
5.1 有限样本偏差的修正
当样本量只有几十个时，psi(nx + 1)中的nx经常为零，digamma 函数在零附近发散，导致互信息估计出现负值。负的互信息在理论上不存在，但在有限样本估计中经常发生，尤其在 k=1 时。常见处理方式有两种：一种是把nx和ny加上一个小常数后再代入psi；另一种是改用 KSG2 公式，它带有额外的修正项，可以减少负值出现的概率。实际应用中我倾向于先检查负值比例，如果超过 10%，就把 k 从 1 提高到 3。
MATLAB 的psi函数没有对参数为零做特殊处理，因此调用前必须保证nx >= 1。在KraskovMI.m内部，可以加入一行安全约束：
nx = max(nx, 1); ny = max(ny, 1);这行代码的作用是把计数下限钳制到 1，避免 digamma 函数输入为零。代价是会在极端稀疏区域产生轻微向上的偏差，但相比负值结果更可解释。
5.2 批量计算时的缓存与向量化
循环调用knnsearch是最大的性能瓶颈。当特征数较多时，所有特征的联合坐标不同，无法一次性完成搜索。这里可以做一个折中：先对所有特征做归一化，再分别与目标变量拼接计算。如果目标变量保持不变，可以在外层循环外预先计算目标变量的自近邻距离，内层只计算特征与目标的交叉距离。另一种加速方式是使用parfor并行计算特征维度：
parfor f = 1:numFeatures miScores(f) = KraskovMI(X(:, f), Y, kFixed); end启用parfor前需要确保KraskovMI内部不依赖全局变量，knnsearch也能在并行池中正常工作。MATLAB 并行池首次启动会耗时 10 到 30 秒，特征数少于 20 时不建议使用。
5.3 验证估计稳定性的双 k 检验法
最后给一个实用技巧：用两套差异较大的 k 值分别计算互信息，再检查排序变化程度。具体做法是令k1 = 3、k2 = 10，分别得到两个分数向量，计算 Spearman 秩相关系数：
mi3 = computeMI(X, Y, 3); mi10 = computeMI(X, Y, 10); rho = corr(mi3, mi10, 'Type', 'Spearman');rho大于 0.9 说明互信息估计稳定，此时无论是特征排序还是筛选阈值都可以放心用。rho在 0.7 到 0.9 之间说明部分特征受到样本噪声干扰，建议对原始数据重新做去异常值处理后再跑一次。rho低于 0.7 则要怀疑数据中存在长尾分布或缺失编码错误，此时不要继续做特征筛选，先回到数据清洗环节。这个检验只需几分钟，却能把后续模型调参阶段的大量试错成本提前消化掉。
本文还有配套的精品资源，点击获取
