RDKit 机器学习数据统计模块 `rdkit.ML.Data.Stats` 源码级解析

• 科学计算
• 科研
• 机器学习

rdkit.ML.Data.Stats是 RDKit 机器学习子包（rdkit.ML）中用于数据矩阵统计分析的核心纯 Python 模块，提供数据标准化、协方差/相关系数矩阵、主成分分析（PCA）、均值标准差、拟合优度 R² 与置信区间估计等经典统计工具。本文以 Docs/Book/source/rdkit.ML.Data.Stats.rst 所索引的模块文档为骨架，结合 rdkit/ML/Data/Stats.py 源码、rdkit/ML/Data/UnitTestStats.py 单元测试以及仓库内的真实调用点，逐函数讲解其数学原理、实现细节与实战用法，帮助读者在 RDKit 的 ML 数据流水线中正确使用这套统计工具箱。

模块定位：RDKit ML 数据流水线中的统计基础设施

在 RDKit 中，rdkit.ML.Data包负责机器学习建模前的数据准备，其子模块划分见 Docs/Book/source/rdkit.ML.Data.rst：

• rdkit.ML.Data.MLData：MLDataSet/MLQuantDataSet数据集容器，管理样本点、变量、结果列与量化边界；
• rdkit.ML.Data.Quantize 与cQuantize：基于信息增益自动搜索量化边界；
• rdkit.ML.Data.SplitData：训练集/留出集划分；
• rdkit.ML.Data.Stats（本文主角）：对数据矩阵实施标准化、协方差/相关分析、主成分分析等统计变换；
• rdkit.ML.Data.Transforms：声明 Center / Normalize / Standardize 三类数据变换入口。

从源码结构看，Stats.py从 2001 年起由 Greg Landrum 维护（文件头版权声明），全部基于numpy实现，不依赖任何 RDKit C++ 扩展，因此可以脱离分子处理独立使用，也天然适合对分子描述符矩阵（descriptor matrix）做建模前的预处理。模块顶层共暴露 8 个函数：

数据标准化：StandardizeMatrix(mat)

StandardizeMatrix是标准的"减去均值、除以标准差"（z-score）变换，对应 rdkit/ML/Data/Stats.py：

def StandardizeMatrix(mat): nObjs = len(mat) avgs = sum(mat, 0) / float(nObjs) mat -= avgs devs = math.sqrt(sum(mat * mat, 0) / (float(nObjs - 1))) try: newMat = mat / devs except OverflowError: newMat = numpy.zeros(mat.shape, 'd') for i in range(mat.shape[1]): if devs[i] != 0.0: newMat[:, i] = mat[:, i] / devs[i] return newMat

实现要点：

1. 按列（变量）计算：矩阵按行存放样本，sum(mat, 0)得到每个变量的和，除以样本数得到列均值；
2. 标准差使用n - 1分母：即样本标准差（sample standard deviation），这是无偏估计的常用选择；
3. 原地修改与返回值并存：docstring 明确警告mat会被原地减去均值（mat -= avgs），返回值newMat是标准化后的新数组。调用方如果不希望原始数据被破坏，务必先拷贝；
4. 异常兜底：向量化除法mat / devs可能触发OverflowError（例如数据包含极大值或零标准差时），此时退化为逐列除法，且跳过标准差为 0 的列（保持 0 值不除），避免产生 NaN/Inf。

实战调用：聚类前的标准化

仓库内rdkit.ML.Cluster包的聚类标准化实现直接复用了本函数，见 rdkit/ML/Cluster/Standardize.py：

from rdkit.ML.Data import Stats def StdDev(mat): """ the standard deviation classifier This uses _ML.Data.Stats.StandardizeMatrix()_ to do the work """ return Stats.StandardizeMatrix(mat) methods = [ ("None", lambda x: x, "No Standardization"), ("Standard Deviation", StdDev, "Use the standard deviation"), ]

即：在分子描述符聚类流程中，用户可通过methods选择"Standard Deviation"策略，底层即调用Stats.StandardizeMatrix，保证不同量纲的描述符（如分子量、logP、TPSA）对距离度量贡献均衡。

协方差与相关系数矩阵

FormCovarianceMatrix(mat)

rdkit/ML/Data/Stats.py 中实现：

def FormCovarianceMatrix(mat): nPts = mat.shape[0] sumVect = sum(mat) sumVect /= float(nPts) for row in mat: row -= sumVect return numpy.dot(numpy.transpose(mat), mat) / (nPts - 1)

算法即教科书式定义：先对每行减去列均值（同样会原地修改mat），再计算XᵀX / (n-1)，分母使用n-1得到样本协方差矩阵的无偏估计。对角线元素即各变量的方差，非对角线元素为变量两两间的协方差。

该函数在单元测试 rdkit/ML/Data/UnitTestStats.py 中有精确验证，测试数据取自 NIST Engineering Statistics Handbook 的经典算例：

d = numpy.array([[4., 2, 0.6], [4.2, 2.1, 0.59], [3.9, 2.0, 0.58], [4.3, 2.1, 0.62], [4.1, 2.2, 0.63]]) m = Stats.FormCovarianceMatrix(d) target = numpy.array([[0.025, 0.0075, 0.00175], [0.0075, 0.007, 0.00135], [0.00175, 0.00135, 0.00043]]) diff = abs(m - target) assert max(diff.ravel()) < FLOAT_TOL # FLOAT_TOL = 1e-4

FormCorrelationMatrix(mat)

rdkit/ML/Data/Stats.py 按Pearson 相关系数公式逐对计算，而不是简单地由协方差归一化而来：

numerator = N * sum(x * y) - sumX * sumY denom = numpy.sqrt((N * sumX2 - sumX**2) * (N * sumY2 - sumY**2)) if denom != 0.0: res[i, j] = numerator / denom res[j, i] = numerator / denom else: res[i, j] = 0 res[j, i] = 0

要点：

• 结果矩阵对称，对角线恒为 1（i == j时分子分母相等，但实现中未显式赋值对角线，需注意对角元素由公式自然算出）；
• 常数变量（方差为 0）时分母为 0，相关系数被置 0，避免除零异常；
• 该实现正是R2函数 docstring 中引用的 Excel 相关系数公式的矩阵版（见下文）。

测试 rdkit/ML/Data/UnitTestStats.py 用 10×3 的样本数据验证了相关矩阵数值，容差 1e-4。

主成分分析：PrincipalComponents与TransformPoints

PrincipalComponents(mat, reverseOrder=1)

rdkit/ML/Data/Stats.py 实现 PCA 的核心步骤：

def PrincipalComponents(mat, reverseOrder=1): covMat = FormCorrelationMatrix(mat) eigenVals, eigenVects = numpy.linalg.eig(covMat) eigenVals = getattr(eigenVals, "real", eigenVals) eigenVects = getattr(eigenVects, "real", eigenVects) ptOrder = numpy.argsort(eigenVals).tolist() if reverseOrder: ptOrder.reverse() eigenVals = numpy.array([eigenVals[x] for x in ptOrder]) eigenVects = numpy.array([eigenVects[x] for x in ptOrder]) return eigenVals, eigenVects

需要特别说明的实现细节：

1. 基于相关系数矩阵而非协方差矩阵做特征分解（covMat = FormCorrelationMatrix(mat)）。这意味着 PCA 在分析前已对各变量做了标准化，消除了量纲影响——对于分子描述符这种数值跨度很大的数据（如 TPSA 0~200、logP -5~10）尤为合适；
2. 使用numpy.linalg.eig后通过getattr(..., "real", ...)提取实部，兼容复数返回类型；
3. 按特征值排序：默认reverseOrder=1按特征值从大到小排序（第一主成分对应最大方差方向）；传reverseOrder=0则升序。返回的特征向量矩阵行向量与排序后的特征值一一对应。

单元测试 rdkit/ML/Data/UnitTestStats.py 验证了特征值与特征向量的数值正确性，并特别允许特征向量整体取反（eVects[i] + tVect[i]也在容差内），这是特征向量符号不确定性的标准处理方式。

TransformPoints(tFormMat, pts)

rdkit/ML/Data/Stats.py 将原始样本点投影到主成分坐标系：

def TransformPoints(tFormMat, pts): pts = numpy.array(pts) nPts = len(pts) avgP = sum(pts) / nPts pts = pts - avgP res = [None] * nPts for i in range(nPts): res[i] = numpy.dot(tFormMat, pts[i]) return res

要点：

• 投影前先减去样本均值（centering），符合"主成分坐标系原点位于数据中心"的几何意义；
• tFormMat通常即PrincipalComponents返回的特征向量矩阵（行向量），pts可以是二维数组或 numpy 数组列表，返回与输入等长的 numpy 数组列表；
• 一个便捷用法：只取特征向量矩阵的前 k 行传入，即可实现k 维降维投影。

PCA 完整流水线示例（结合单元测试）

rdkit/ML/Data/UnitTestStats.py 展示了"计算主成分 → 投影 → 校验"的完整链路：

_, eVects = Stats.PrincipalComponents(self.d) # 得到排序后的特征向量 pts = Stats.TransformPoints(eVects, self.d) # 投影到 PCA 坐标系 # 与参考投影坐标逐点比对，容差 1e-4

而 UnitTestStats.py 的testTransform2还验证了一个重要几何性质：旋转到 PCA 坐标系不改变任意两点间的归一化点积（夹角余弦），说明该变换是保角正交变换。

描述性统计：MeanAndDev与R2

MeanAndDev(vect, sampleSD=1)

rdkit/ML/Data/Stats.py 返回向量的均值与标准差：

def MeanAndDev(vect, sampleSD=1): vect = numpy.array(vect, 'd') n = vect.shape[0] if n <= 0: return 0., 0. mean = sum(vect) / n v = vect - mean if n > 1: if sampleSD: dev = numpy.sqrt(sum(v * v) / (n - 1)) # 样本标准差 else: dev = numpy.sqrt(sum(v * v) / (n)) # 总体标准差 else: dev = 0 return mean, dev

参数语义：

• sampleSD=1（默认）：分母为n-1，返回样本标准差；
• sampleSD=0：分母为n，返回总体标准差；
• 空向量或单元素向量返回标准差 0，保证数值稳定。

该函数在仓库内有一处典型的真实调用：药效团嵌入能量统计。见 rdkit/Chem/Pharm3D/EmbedLib.py，EmbedPharmacophoreToMol对多次嵌入得到的构象能量e1s/e2s/e3s/e4s调用Stats.MeanAndDev，汇总输出形如name(failed): e1(σ1) -> e2(σ2) : e3(σ3) -> e4(σ4)的能量统计报告，用于评估药效团约束嵌入的收敛质量。

R2(orig, residSum)

rdkit/ML/Data/Stats.py 计算拟合优度：

def R2(orig, residSum): vect = numpy.array(orig) n = vect.shape[0] if n <= 0: return 0., 0. oMean = sum(vect) / n v = vect - oMean oVar = sum(v * v) return 1. - residSum / oVar

即R² = 1 - 残差平方和 / 总平方和。注意源码 docstring 中的 FIX 备注：作者明确指出该简化实现"just is not right"（并不完全正确），并给出多变量情形下正确的 Pearson 平方形式：

r2 = [n*(Sxy) - (Sx)(Sy)]^2 / ([n*(Sx2) - (Sx)^2]*[n*(Sy2) - (Sy)^2])

因此引用本函数时应将其视为"基于给定残差和的 1 - SS_res/SS_tot 定义"的启发式指标，用于预测质量的大致评估，而不应视为严格的多变量决定系数。这是源码自带的文档化警告，使用者务必留意。

置信区间：GetConfidenceInterval(sd, n, level=95)

rdkit/ML/Data/Stats.py 依据t 分布估计均值的置信区间半宽：

def GetConfidenceInterval(sd, n, level=95): col = tConfs[level] dofs = n - 1 sem = sd / numpy.sqrt(n) idx = 0 while idx < len(tTable) and tTable[idx][0] < dofs: idx += 1 if idx < len(tTable): t = tTable[idx][col] else: t = tTable[-1][col] return t * sem

实现机制：

• t 临界值表内嵌在模块中：rdkit/ML/Data/Stats.py 定义了tConfs置信水平映射与tTable查找表（自由度 1~200+，单尾显著性水平 0.10~0.0005）；
• tConfs = {80: 1, 90: 2, 95: 3, 98: 4, 99: 5, 99.8: 6, 99.9: 7}，即支持 80%、90%、95%、98%、99%、99.8%、99.9% 七档置信水平；
• 计算过程：自由度dofs = n - 1，标准误sem = sd / sqrt(n)，在表中查取对应自由度的 t 临界值（超过表末自由度时回退到表末行，即大样本近似），返回t * sem；
• 返回的是区间半宽（margin of error），真实置信区间为(mean - 半宽, mean + 半宽)。

用法示例（假设某描述符的样本均值为 m、样本标准差 sd、样本量 n）：

from rdkit.ML.Data import Stats mean = 2.5 sd, _ = Stats.MeanAndDev([2.1, 2.6, 2.3, 2.7, 2.5, 2.4, 2.8, 2.2, 2.6, 2.5]) half = Stats.GetConfidenceInterval(sd, n=10, level=95) print(f"95% CI: ({mean - half:.4f}, {mean + half:.4f})")

测试验证与版本沿革：从rdkit.ML.MatOps到rdkit.ML.Data.Stats

模块的测试与历史迁移情况可以从仓库源码中确认：

单元测试：rdkit/ML/Data/UnitTestStats.py 以 1e-4 容差覆盖了协方差矩阵（NIST 参考算例）、相关系数矩阵、PCA 特征值与特征向量（允许整体符号翻转）、PCA 投影坐标以及保角性质共 5 组用例，是理解各函数数值行为的权威参考。

历史迁移：rdkit/ML/MatOps.py 的 docstring 明确写道"the two functions defined here have been moved to ML.Data.Stats"，并在文件底部以兼容性别名方式重新导出：

from rdkit.ML.Data import Stats FormCovarianceMatrix = Stats.FormCovarianceMatrix PrincipalComponents = Stats.PrincipalComponents

这说明当前Stats.py是 RDKit 矩阵统计功能的权威实现，rdkit.ML.MatOps仅是保留旧接口的转发层。相关文档见 Docs/Book/source/rdkit.ML.MatOps.rst。

与数据流水线的衔接：在完整建模流程中，Stats通常与 rdkit.ML.Data.MLData 的MLDataSet/MLQuantDataSet容器、rdkit.ML.Data.Quantize 的量化以及 rdkit.ML.Data.SplitData 的训练/留出划分配合使用，典型顺序为：构建描述符矩阵 → 标准化/相关分析/PCA 降维 → 量化 → 划分数据集 → 建模评估（R²）。

结语

rdkit.ML.Data.Stats是 RDKit 机器学习栈中一块轻量而完整的统计基础设施：StandardizeMatrix支撑聚类前的量纲归一，FormCovarianceMatrix/FormCorrelationMatrix提供变量关系诊断，PrincipalComponents/TransformPoints构成基于相关矩阵的 PCA 降维链路，MeanAndDev与R2服务回归评估，GetConfidenceInterval借助内嵌 t 表给出区间估计。使用时注意三点：StandardizeMatrix与FormCovarianceMatrix会原地修改输入矩阵、PrincipalComponents基于相关系数矩阵（量纲无关）、R2的实现存在源码自注的正确性限制。深入理解这些细节，能让你在 RDKit 描述符建模管线中少踩数据预处理与统计推断的坑。

• 科学计算
• 科研
• 机器学习