拆解格物统计的统计内核：自研数值算法、语法引擎与 388 项断言

YoLo♪

294人浏览 · 2026-11-10 11:06:06

YoLo♪ · 2026-11-10 11:06:06 发布

欢迎加入开源鸿蒙 PC 社区：https://harmonypc.csdn.net/

欢迎在 PC 社区平台申请新建项目：https://atomgit.com/OpenHarmonyPCDeveloper

开源地址（AtomGit 仓库）：https://atomgit.com/OpenHarmonyPCDeveloper/ohos_PSPP.git

开源组织（AtomGit）：https://atomgit.com/OpenHarmonyPCDeveloper

上一篇《格物统计：在鸿蒙 PC 上原生重建一套统计分析工作台》讲的是产品形态、界面分层、数据层与 NAPI 边界。这一篇往下钻一层，专门拆 libpsppcore.so 里的统计内核：为什么数值层必须自己写、矩统计与分布函数怎么实现、回归为什么用扫描算子、Logistic 的 IRLS 怎么收敛、因子分析的 Jacobi 与 Varimax 怎么写、K-Means 与 Ward 怎么做确定性输出，以及 388 项断言是怎么和 numpy / scipy 对齐的。

一句话概括：界面可以重写，统计公式不能重造；但公式背后的数值基础设施，必须能在 OHOS 上跑起来。

一、为什么数值层必须自己写

上游 GNU PSPP 的统计能力建立在一组成熟的 C 依赖上：GSL 提供分布函数、特殊函数与线性代数，GLib 提供容器与基础设施，GTK 提供界面。鸿蒙 PC 的 OHOS NDK 里没有 GSL，也没有 GTK；把上游 src/math 直接交叉编译过来，会拖进一整条依赖链。

所以这个项目的策略是：

算法与公式：以上游的统计定义、数据结构、SAV 文件格式为参照；
实现：用 C++17 重新写一遍，只依赖 C++ 标准库；
数值依赖：用自研 math/ 替代 GSL 中真正用到的那部分；
平台依赖：native 内核不调用任何 OHOS API，平台相关能力全部留给 ArkTS。

这样做的直接收益是可测试：同一套 C++ 源码可以在开发机上用 clang++ 编译成宿主端测试程序，不需要模拟器、不需要 HDC，也不需要每次打包 HAP。388 项断言之所以能快速迭代，靠的就是这条「内核与 UI 解耦」的边界。

上游模块与本项目实现的对应关系：

上游能力	上游依赖	本项目实现
矩统计、描述统计	`src/math/moments.c`	`math/moments.cpp`（两遍法 + 在线法）
分位数、中位数、众数	`src/math/percentiles.c`	`math/order_stats.cpp`
分布函数、特殊函数	GSL	`math/math_compat.cpp`（Lanczos、不完全伽马 / 贝塔）
相关、回归	`src/math/correlation.c`、`lib/linreg`	`math/correlation.cpp`、`math/linreg.cpp`、`math/matrix.h`
因子分析	GSL 特征分解	`math/factor.cpp`（Jacobi + Varimax）
聚类	上游 QUICK CLUSTER	`math/cluster.cpp`（K-Means + Ward）
Logistic	上游 Logistic 过程	`math/logistic.cpp`（IRLS）
数据与文件格式	`src/data/`	`core/`、`io/data_io.cpp`、`io/sav_io.cpp`

整条调用链是：

ArkTS 分析面板 / SyntaxEditor
        ↓ NAPI
napi_init.cpp（参数解析 + 全局 Dataset）
        ↓
stats/engine（cmd_descriptives / cmd_ttest / cmd_factor …）
        ↓
math/*（moments / order_stats / math_compat / linreg / factor / cluster / logistic）
        ↓
core/Dataset（列式数据 + 可见行映射）
        ↓
output/AnalysisResult（Table[] + Chart[]）
        ↓ JSON
ArkTS OutputViewer / ChartRenderer

stats/ 只依赖 core/ 与 math/，不感知 UI、文件系统与平台；math/ 只依赖标准库。这个分层让「统计过程」与「数值算法」都能独立测试。

二、矩统计：两遍法与在线法

描述统计里的均值、方差、偏度、峰度看起来简单，但在大样本与加权场景下很容易写出数值不稳定的实现。moments.cpp 提供了两套累积器：

Moments：两遍法。第一遍累加均值，第二遍围绕均值累加中心矩；
Moments1：单遍在线更新。每个新值到来时递推更新一至四阶中心矩。

两遍法的核心公式在 calc_moments 里：

void calc_moments(Moment max_moment, double w, double d1, double d2, double d3,
                  double d4, double* variance, double* skewness, double* kurtosis) {
    if (max_moment >= Moment::Variance && w > 1.) {
        double s2 = (d2 - pow2(d1) / w) / (w - 1.);
        if (variance) *variance = s2;

        if (std::fabs(s2) >= 1e-20) {
            if (max_moment >= Moment::Skewness && skewness && w > 2.) {
                double s3 = s2 * std::sqrt(s2);
                double g1 = (w * d3) / ((w - 1.0) * (w - 2.0) * s3);
                if (std::isfinite(g1)) *skewness = g1;
            }
            if (max_moment >= Moment::Kurtosis && kurtosis && w > 3.) {
                double den = (w - 2.) * (w - 3.) * pow2(s2);
                double g2 = (w * (w + 1) * d4 / (w - 1.) / den -
                             3. * pow2(d2) / den);
                if (std::isfinite(g2)) *kurtosis = g2;
            }
        }
    }
}

几个值得注意的点：

权重不是 1 的时候走加权公式：w 是个案权重之和，不是简单计数；
偏度要求 w > 2、峰度要求 w > 3：样本量不足时返回缺失值，而不是硬算出一个无意义的数；
每一步都做有限性检查：std::isfinite 拦截溢出与 NaN，避免坏值污染整张结果表。

在线法 Moments1::add 则用递推公式维护 d1_..d4_，适合只需要扫一遍数据的场景，例如分组均值、聚合与部分描述统计。两种实现共享同一个 calc_moments 出口，保证公式口径一致。

顺序统计量在 order_stats.cpp 里实现，百分位算法与上游 percentiles.c 对齐：

算法	位置定义	用途
`HAverage`	`g = (n+1)p`，线性插值	SPSS / PSPP 默认
`WAverage`	`g = n*p`，线性插值	另一种加权口径
`Round`	`g = (n+1)p`，就近取整	离散分位
`Empirical`	`g = n*p`，向上取整	经验分布

中位数走 HAverage(0.5)，众数按排序后的等值段统计最高频，箱线图的 Tukey hinges 也在这里实现。把这些「小算法」单独收在一层，统计过程就不用各自重复实现分位数。

三、分布函数与特殊函数：不引 GSL，自己算

p 值、置信区间、卡方检验、F 检验都离不开分布函数。math_compat.cpp 实现了统计过程真正需要的那一小部分：

Γ 函数：Lanczos 近似（g=7, n=9），负数走反射公式；
不完全伽马：gser 级数展开 + gcf 连分式展开，按 x < a+1 切换；
不完全贝塔：betacf 连分式，按对称性选择正反方向；
正态 CDF：0.5 * erfc(-x / (σ√2))；
t / 卡方 / F 分布 CDF：分别归约到不完全贝塔与不完全伽马；
分位数：二分法迭代 200 次求逆。

double gamma_p(double a, double x) {
    if (x < 0.0 || a <= 0.0) return SYSMIS;
    if (x == 0.0) return 0.0;
    if (x < a + 1.0) return gser(a, x);
    return 1.0 - gcf(a, x);
}

double cdf_tdist_P(double t, double nu) {
    if (nu <= 0.0) return SYSMIS;
    double x = nu / (nu + t * t);
    double ib = 0.5 * beta_inc(nu / 2.0, 0.5, x);
    return t > 0.0 ? 1.0 - ib : ib;
}

数值实现里有几处刻意的保守处理：

kEps = 3e-15、kFpMin = 1e-300、最大迭代 500 次，避免死循环；
连分式中每次更新后都把过小的 d / c 夹回 kFpMin，防止除零；
分位数用二分法而非解析近似，牺牲一点速度换稳定；
p 值格式化时，极小的 p 统一输出 <.001，不显示一串无意义的 0.000000。

这里没有追求「实现整个 GSL」，只实现统计过程实际用到的函数。依赖面越小，交叉编译与回归测试的成本越低。

四、线性代数与回归：扫描算子

线性回归没有用「构造 X 然后求 (X'X)^{-1}X'y」这种直白写法，而是用了扫描算子（sweep operator）。math/matrix.h 里实现了一个轻量稠密矩阵 Matrix 和 sweep_first（下面为节选）：

inline void sweep_first(Matrix& A, int k) {
    for (int d = 0; d < k; ++d) {
        double piv = A.at(d, d);
        if (std::fabs(piv) < 1e-300) continue;
        Matrix B(A);
        B.at(d, d) = -1.0 / piv;
        for (int i = 0; i < A.rows(); ++i) {
            if (i == d) continue;
            for (int j = 0; j < A.rows(); ++j) {
                if (j == d) continue;
                B.at(i, j) = A.at(i, j) - A.at(i, d) * A.at(d, j) / piv;
            }
        }
        // 更新第 d 列 / 行……
        A = B;
    }
}

回归的增广矩阵是：

[ X'X   X'y ]
[ y'X   y'y ]

扫描前 k 个对角元之后：

右上角的 X'y 位置变成回归系数；
右下角变成误差平方和 SSE；
左上角取负就是 (X'X)^{-1}，用来算系数标准误。

linreg_fit 再据此计算 R²、调整 R²、模型 / 残差自由度、MSE、系数 t 值与整体 F 检验。单自变量时 R 取斜率符号，多元回归取非负——这是统计软件输出口径的细节。

Logistic 回归则用了另一套线性代数：solve_inplace 做部分主元高斯消元解 A x = b，invert 做高斯-约当求逆拿协方差矩阵。两者都在 math/logistic.cpp 里，矩阵规模是「自变量数 + 1」，不大，但要求数值稳定。

因子分析用 Jacobi 特征分解。对对称相关矩阵做一系列平面旋转，直到非对角元素平方和小于 1e-24：

const double theta = (a[q][q] - a[p][p]) / (2.0 * a[p][q]);
const double t = ((theta >= 0.0) ? 1.0 : -1.0) /
                 (std::fabs(theta) + std::sqrt(theta * theta + 1.0));
const double c = 1.0 / std::sqrt(t * t + 1.0);
const double s = t * c;

Jacobi 的优点是实现简单、对对称矩阵稳定；变量数通常是个位到几十，性能完全够用。

五、统计过程：从公式到结果表

stats/engine.h 暴露了一组统一签名的过程，全部返回 AnalysisResult：

AnalysisResult cmd_descriptives(Dataset&, const std::vector<int>&, const std::vector<std::string>&);
AnalysisResult cmd_frequencies(Dataset&, const std::vector<int>&, bool);
AnalysisResult cmd_crosstabs(Dataset&, int, int);
AnalysisResult cmd_correlations(Dataset&, const std::vector<int>&, bool);
AnalysisResult cmd_ttest_one_sample(Dataset&, const std::vector<int>&, double);
AnalysisResult cmd_ttest_independent(Dataset&, int, double, double, const std::vector<int>&);
AnalysisResult cmd_ttest_paired(Dataset&, const std::vector<std::pair<int, int>>&);
AnalysisResult cmd_oneway(Dataset&, int, int);
AnalysisResult cmd_regression(Dataset&, int, const std::vector<int>&);
AnalysisResult cmd_logistic(Dataset&, int, const std::vector<int>&);
AnalysisResult cmd_explore(Dataset&, const std::vector<int>&, int);
AnalysisResult cmd_means(Dataset&, const std::vector<int>&, int, bool);
AnalysisResult cmd_reliability(Dataset&, const std::vector<int>&);
AnalysisResult cmd_factor(Dataset&, const std::vector<int>&, int, bool);
AnalysisResult cmd_cluster(Dataset&, const std::vector<int>&, int, const std::string&);

统一签名带来两个好处：GUI 与 Syntax 可以共用同一组函数；新增过程只需要在 engine.h 加声明、在 NAPI 加一个导出、在面板加一个入口。

5.1 相关分析：Pearson 与 Spearman

相关分析输出三张表：相关系数矩阵、双尾显著性、样本量矩阵，并在前两个变量上自动生成散点图。Pearson 走协方差 / 标准差；Spearman 先对每个变量取秩（并列取平均秩），再对秩做 Pearson。显著性用 t 分布：

t = r * sqrt((n - 2) / (1 - r * r));
p = 2 * (1 - cdf_tdist_P(fabs(t), n - 2));

相关矩阵采用列表删除：只要任一变量缺失，该个案整体排除，因此矩阵各格的 N 相同。结果里仍然单独给一张样本量矩阵，是为了让用户明确知道实际参与计算的样本量。

5.2 T 检验：Levene + 合并方差 + Welch

独立样本 T 检验不是只给一个 t 值。代码里同时做了：

Levene 检验（以均值为中心）：把每组值转成 |x - mean|，再做方差分析，得到 Levene F 与显著性；
假定等方差：用合并方差 sp² 计算 t，自由度 n1+n2-2；
不假定等方差：Welch t，自由度用 Welch-Satterthwaite 公式。

double sp2 = ((n1 - 1) * s1.variance + (n2 - 1) * s2.variance) / (n1 + n2 - 2);
double se_pooled = std::sqrt(sp2 * (1.0 / n1 + 1.0 / n2));
double t_pooled = (s1.mean - s2.mean) / se_pooled;
double df_pooled = n1 + n2 - 2;

double a = s1.variance / n1, b = s2.variance / n2;
double se_w = std::sqrt(a + b);
double t_w = (s1.mean - s2.mean) / se_w;
double df_w = pow2(a + b) / (pow2(a) / (n1 - 1) + pow2(b) / (n2 - 1));

结果表里两行并列输出，让用户自己决定看哪一行。单样本检验的 95% 置信区间用 idf_tdist(0.975, df) 求临界值；配对样本则先做差值，再对差值做单样本 t 检验，并额外给一张配对相关表。

在这里插入图片描述

5.3 线性回归

回归结果按标准三表组织：

模型摘要：R、R²、调整 R²、标准误；
ANOVA：回归 / 残差 / 总计的平方和、自由度、均方、F、显著性；
系数表：B、标准误、t、显著性。

这些值都来自扫描算子的同一份增广矩阵，不需要再单独求一次逆。多自变量共线时，sweep_first 遇到极小主元会跳过，结果可能变得不稳定，由上层决定如何提示。

在这里插入图片描述

5.4 二元 Logistic：IRLS 与拟合优度

Logistic 回归用迭代重加权最小二乘（IRLS）。每一轮：

eta = Xβ
mu  = 1 / (1 + exp(-eta))
w   = mu * (1 - mu)
z   = eta + (y - mu) / w
解 (X'WX) β_new = X'Wz

代码里对 mu 做了夹紧（1e-10 ~ 1-1e-10），避免 log(0) 与权重为 0：

double clamp_mu(double eta) {
    double mu = 1.0 / (1.0 + std::exp(-eta));
    if (mu < 1e-10) mu = 1e-10;
    if (mu > 1.0 - 1e-10) mu = 1.0 - 1e-10;
    return mu;
}

收敛条件是系数最大变化 < 1e-8，最多 50 轮。收敛后再用最终信息矩阵 X'WX 求逆得到协方差：

SE = sqrt(cov[j][j])
Wald = (β / SE)²
Sig = 2 * (1 - Φ(|β / SE|))
Exp(B) = exp(β)

模型摘要给 -2LL、Cox-Snell R²、Nagelkerke R²；分类表按 0.5 阈值统计；Hosmer-Lemeshow 按预测概率排序分成 10 组，df = 组数 - 2。如果某个系数的标准误过大或 |β| > 10，结果会标记「可能存在完全分离」——这是 Logistic 在真机上最常见的异常，提前识别比事后猜错更有用。

5.5 信度分析：Cronbach α 与「删除该项后的 α」

信度分析采用列表删除：任一项缺失的个案整体排除。核心公式：

α = k / (k - 1) * (1 - Σ item_variance / total_variance)

标准化 α 则用平均项间相关：

α_std = k * rbar / (1 + (k - 1) * rbar)

项总计统计量里还有三项常用指标：删除该项后的尺度均值、删除该项后的尺度方差、校正的项总计相关性，以及删除该项后的 α。后者就是把当前项拿掉再算一次 α，用来判断「哪一项拉低了整体信度」。

5.6 非参数检验：结校正与渐近显著性

实现了四个过程：

过程	统计量	说明
卡方拟合优度	`χ²`	期望比例可自定义，默认等概率
Mann-Whitney U	`U` / `Z`	两组秩和，渐近正态
Kruskal-Wallis H	`H`	多组秩和，含结校正
Wilcoxon 符号秩	`W` / `Z`	配对差值秩，含结校正

结（tie）的处理很关键：秩相同的时候取平均秩，并在方差里减去结校正项。否则同一份数据在「有并列值」与「无并列值」两种情况下会得到不同显著性，用户很难解释。

5.7 探索性分析：Shapiro-Wilk 与茎叶图

EXPLORE 除了描述统计（含 95% 均值置信区间、四分位距），还做：

Shapiro-Wilk 正态性检验（Royston AS R94），适用 4 ≤ n ≤ 5000，方差为 0 时直接返回失败；
茎叶图：未分组时输出；
箱线图：含四分位、中位数、须线与离群点。

Shapiro-Wilk 里用 Acklam 近似求标准正态分位数，再按 Royston 的系数构造 W 统计量。这个实现不依赖任何外部库，因此探索性分析也能在宿主端测试。

5.8 因子分析：Jacobi + Kaiser + Varimax

因子分析的完整链路是：

列表删除 → 标准化 → 相关矩阵 R
  → Jacobi 特征分解（特征值降序、特征向量）
  → Kaiser 准则（特征值 > 1）或用户指定因子数
  → 载荷 = 特征向量 × sqrt(特征值)
  → 共同度 = 载荷平方和
  → 可选 Kaiser 标准化 Varimax 旋转
  → 旋转后载荷平方和
  → KMO / Bartlett

Varimax 的目标是让每个变量在尽量少的因子上有高载荷。实现上先按行做 Kaiser 标准化（每行载荷除以该行平方和开方），再对每一对因子做平面旋转，迭代到最大旋转角 < 1e-7，最后把标准化还原并按平方和降序重排因子，保证输出稳定。

KMO 用相关矩阵与逆相关矩阵计算：

KMO = Σ r_ij² / (Σ r_ij² + Σ q_ij²)

其中 q_ij 是偏相关（由逆矩阵元素归一化得到）。Bartlett 球形检验则用特征值对数行列式：

χ² = -(n - 1 - (2k + 5) / 6) * ln|R|
df = k(k - 1) / 2

结果页会给出 KMO、Bartlett、共同度、方差解释（含碎石图）、成分矩阵，旋转后再生一张旋转成分矩阵。
在这里插入图片描述

5.9 聚类分析：确定性的 K-Means 与 Ward

K-Means 最容易出现「每次跑结果不一样」的问题。这里用 maximin 确定性初始化：第一个中心取第一个个案，之后每次选「离已有中心最近距离最大」的个案作为新中心。这样同一份数据、同一组参数，中心与簇编号都是确定的。

Lloyd 迭代里额外做了空簇修复：如果某个簇没有个案，就从「离自身中心最远且所在簇成员数大于 1」的个案里移一个过去，避免空簇导致中心无法更新。

收敛后按「簇大小降序、中心字典序」重排簇编号，保证输出稳定。SSE（簇内平方和）作为聚类摘要的一部分输出。

Ward 层次聚类用欧氏距离 + Lance-Williams 递推。每次合并两个最近的簇，更新到其他簇的距离：

d(u, v)² = ((nv + ni) / t) * d(v, si)²
         + ((nv + nj) / t) * d(v, sj)²
         -  (nv / t) * d(si, sj)²
t = nv + ni + nj

合并过程记录成 MergeStep 序列，cut_tree 用并查集从完整树开始撤销最后 k-1 次合并，得到目标簇数的归属。这样既支持 K-Means，也支持 Ward，并且都能输出聚类中心、各簇个案数、凝聚顺序与归属表。

在这里插入图片描述

六、语法引擎：词法、递归下降、共用同一套过程

语法编辑器不是「玩具」：它执行的是同一组 cmd_* 函数，结果也进入同一个结果树。

6.1 词法

tokenize 把输入切成标识符、数字、字符串、标点四类：

标识符允许字母、数字、下划线，以及 @、#、$；
字符串用单引号包裹，内部 '' 表示一个单引号；
数字支持小数点与科学计数法；
注释行以 * 或 COMMENT 开头，整行跳过。

6.2 表达式求值

表达式是递归下降，优先级从低到高：

OR < AND < NOT < 比较 < 加减 < 乘除 < 一元

比较运算有一个容易忽略但要紧的语义：任一侧是系统缺失值时，结果为假。这与 PSPP 的语义一致，也避免缺失值参与比较后产生不可解释的筛选结果。

double rel() {
    double v = expr();
    if (is_punct(">")) {
        // ...
        v = (is_sysmis(v) || is_sysmis(r)) ? 0.0 : ((eq ? v >= r : v > r) ? 1.0 : 0.0);
    }
    // ...
}

6.3 命令分发

run_syntax 先按行去掉注释，再按命令分隔符拆分，取第一个词作为命令名，剩余部分交给对应的解析函数：

if (name == "DESCRIPTIVES") r = do_descriptives(ds, body);
else if (name == "FREQUENCIES") r = do_frequencies(ds, body);
else if (name == "CROSSTABS") r = do_crosstabs(ds, body);
else if (name == "CORRELATIONS") r = do_correlations(ds, body);
else if (name == "T-TEST") r = do_ttest(ds, body);
else if (name == "ONEWAY") r = do_oneway(ds, body);
else if (name == "REGRESSION") r = do_regression(ds, body);
else if (name == "LOGISTIC") {
    // 命令名为两个词：LOGISTIC REGRESSION；剥掉 REGRESSION 后交给解析器
    std::string b = trim(body);
    if (to_upper(b).rfind("REGRESSION", 0) == 0) {
        size_t sp2 = b.find_first_of(" \t");
        b = (sp2 == std::string::npos) ? std::string() : trim(b.substr(sp2 + 1));
    }
    r = do_logistic(ds, b);
}
else if (name == "EXPLORE") r = do_explore(ds, body);
else if (name == "MEANS") r = do_means(ds, body);
else if (name == "RELIABILITY") r = do_reliability(ds, body);
else if (name == "FACTOR") r = do_factor(ds, body);
else if (name == "CLUSTER") r = do_cluster(ds, body);
else if (name == "QUICK") {
    // QUICK CLUSTER：剥掉 CLUSTER 前缀
    std::string b = trim(body);
    if (to_upper(b).rfind("CLUSTER", 0) == 0) b = trim(b.substr(7));
    r = do_cluster(ds, b);
}
else if (name == "NPAR") {
    // NPAR TESTS：剥掉 TESTS 前缀
    std::string b = trim(body);
    if (to_upper(b).rfind("TESTS", 0) == 0) b = trim(b.substr(5));
    r = do_npar(ds, b);
}

已覆盖的命令分三类：

类别	命令
统计分析	`DESCRIPTIVES` `FREQUENCIES` `CROSSTABS` `CORRELATIONS` `T-TEST` `ONEWAY` `MEANS` `REGRESSION` `LOGISTIC REGRESSION` `EXPLORE` `RELIABILITY` `FACTOR` `CLUSTER` `NPAR TESTS`
数据操作	`COMPUTE` `RECODE` `COUNT` `AUTORECODE` `RANK` `AGGREGATE` `SELECT IF` `FILTER IF` `SORT CASES BY` `MISSING VALUES` `VARIABLE LABELS` `VALUE LABELS` `WEIGHT` `LIST`
兼容忽略	`COMMENT` `TITLE` `SET` `SHOW` `PRESERVE` `RESTORE` `FORMATS` `DISPLAY` `DATASET` `OUTPUT` `NEW` `ADD` `APPLY` `MATRIX` 等

COMPUTE 有一个实用细节：如果目标变量不存在，会自动新建一个数值变量，默认格式 F8.2；表达式逐行求值，出错的行跳过。SELECT IF 则是「从后向前删除未选中行」，避免删除过程中下标错位。

GUI 里的「计算变量」「选择个案」「重编码」对话框并不另写一套逻辑，而是把用户输入拼成语法交给 runSyntax——这也是「GUI 与 Syntax 共用同一引擎」的具体落地。

在这里插入图片描述

七、数据变换与二进制格式

7.1 五个变换过程

stats/transforms.cpp 实现了会直接修改当前数据集的五个过程：

过程	行为
`AGGREGATE`	按 1…n 个分组变量聚合，支持 MEAN / SUM / MIN / MAX / SD / N / MEDIAN，结果替换数据集
`COUNT`	统计多个变量命中指定值的次数，写入新变量
`AUTORECODE`	分类取值映射为 1…k 整数代码，并自动生成值标签
`RANK`	平均秩（并列取平均），支持升 / 降序
`WEIGHT`	按频数权重大小四舍五入展开个案，上限 100 万行，等价频数加权

AGGREGATE 的结果替换语义与 PSPP 的 MODE=REPLACE 一致；WEIGHT 之所以用「展开个案」而不是在每条统计公式里乘权重，是为了让后续所有统计过程无需感知权重，行为与上游保持等价。

7.2 SAV：把二进制格式拆开

.sav 的读取分四步：文件头、字典记录、数据区、压缩解码。

文件头 176 字节，包含 magic（$FL2 / $FL3）、case 数、变量数、压缩模式、字节序、浮点格式等。之后是字典记录：变量记录、值标签、文档、扩展记录、结束标记。读取时还要处理「超长字符串」：宽度超过 255 的字符串会拆成多个 252 / 255 字节的段，逻辑变量与物理段要分别计数。

压缩模式 1 是 bytecode，解码逻辑直接按 opcode 还原成「无压缩字节流」：

if (op == 0) {
    // 补齐当前 case
} else if (op == 252) {
    eof = true;
} else if (op == 253) {
    // 8 字节字面量
} else if (op == 254) {
    // 8 个空格
} else if (op == 255) {
    // SYSMIS
} else {
    double val = static_cast<double>(op) - bias;
    // 数值编码
}

压缩模式 2 是 zlib。sav_io 内置了一个 puff 风格的 inflate 实现，支持 stored / fixed / dynamic Huffman 三种块，不依赖系统 zlib。这样做的好处是 native 内核的依赖面更小，HAP 自包含；代价是代码量增加，但换来的是「拿到的 .sav 无论哪种压缩模式都能读」。

写入侧当前输出无压缩 .sav，包含变量字典、值标签与数据区，可被上游工具回读。后续如果要做体积优化，只需要在写入路径补 zlib 压缩，不影响读取路径。

7.3 CSV 与 XLSX

CSV 的分隔符探测与引号解析在 data_io.cpp；XLSX 在 ArkTS 侧完成解压与 XML 解析，再转成 CSV 交给同一条 native 导入管线。这个设计的意义是：格式适配层可以独立演进，统计内核不需要知道文件从哪来。

八、NAPI 边界、内存与性能

8.1 为什么用 JSON 字符串

NAPI 层每个统计函数都返回 JSON 字符串，而不是在 C++ 里构造 napi_value 对象树。原因有三：

边界薄：C++ 只需要一个 json_writer，不需要为每种结果类型维护一套 NAPI 构造代码；
可调试：出问题时直接在 ArkTS 侧看 JSON，定位比翻 napi_value 调用栈快；
可测试：同一套结果序列化逻辑可以在宿主端验证，不必依赖 NAPI 运行时。

代价是多一次序列化，但通过分页与结果表规模控制，实际开销可控。.sav 二进制则走 ArrayBuffer，避免把二进制塞进 JSON。

8.2 内存模型

列式存储：数值列 vector<double>、字符串列 vector<string>，按列连续；
可见行映射：筛选只重建 visible_，不拷贝数据；
分页：getPage(start, count) 只序列化当前页；
结果引用稳定：AnalysisResult 用 std::deque 存表格与图表，追加时已有引用不失效；
错误兜底：每个导出返回都带 ok / error，ArkTS 解析失败时降级为空结果并提示。

10 万行描述统计在宿主端测试里作为性能用例：N=100000，均值 49.5，用来确认列式存储与矩统计在较大样本下仍然稳定。

8.3 构建

entry/src/main/cpp/CMakeLists.txt 用 file(GLOB ... CONFIGURE_DEPENDS) 收集 core / math / output / stats / io / syntax 与 napi_init.cpp，编译成 libpsppcore.so，链接 libace_napi.z.so 与 libhilog_ndk.z.so。CONFIGURE_DEPENDS 是为了避免新增 .cpp 漏进构建导致链接期 undefined symbol——这个坑在加 Logistic、因子、聚类时都踩过。

九、测试：388 项断言是怎么来的

统计软件最怕「看起来有结果，但数不对」。这个项目的测试策略是：同一套源码，两个编译目标。

entry/src/main/cpp/*.cpp
        ├── OHOS NDK → libpsppcore.so → HAP
        └── 宿主 clang++ → test_main → 388 项断言

因为 native 内核不依赖 OHOS API，测试可以直接在 macOS / Linux 上编译运行，不需要模拟器。测试入口：

bash tests/run_all_tests.sh

用例	方法	结果
统计内核对照测试	388 项断言，覆盖描述 / 频数 / 交叉表 / 相关 / T / ANOVA / 回归 / Logistic / 探索 / 分组均值 / 信度 / 非参数 / 因子 / 聚类 / 数据变换 / 排序筛选 / CSV·SAV 回环 / 语法 / 10 万行	PASS
SAV 编解码回环	由 CSV 生成 `.sav`，再逐单元格回读比对	PASS
XLSX 解析	解析真实 xlsx，与期望 CSV 逐字节对拍	PASS
PDF 导出	内置结构校验（xref / 页数 / 中文可还原），可选 `qpdf --check`	PASS

参考值的来源：

描述统计、相关、T、ANOVA、回归、Logistic、信度、非参数：用 numpy / scipy 独立计算；
因子分析：np.linalg.eigh + SVD 版 Varimax 与 Jacobi 内核逐项对拍；
聚类：K-Means 与 numpy、Ward 与 scipy 独立对拍；
SAV：写读回环，不依赖外部工具；
PDF：内置解析器检查 xref、页数与中文还原，外部 qpdf 只作为增强校验。

需要说明边界：388 项断言覆盖的是统计内核的数值与格式正确性，不等于所有 UI 路径都有自动化覆盖；界面交互仍以人工验证与截图记录为准。测试报告与参考值生成脚本在 tests/ 下，可复现。

十、经验与边界

10.1 三条经验

第一，数值层要能脱离平台编译。 这是整个项目最划算的设计决策：统计内核只依赖 C++ 标准库，宿主端就能跑 388 项断言，迭代速度比「打包 → 安装 → 点界面」快一个数量级。

第二，确定性算法优先。 K-Means 用 maximin 初始化、按簇大小与中心字典序重排；Varimax 按平方和降序重排；排序用 stable_sort。统计结果需要可复现，同样的数据与参数必须得到同样的输出。

第三，能自包含就自包含。 SAV 的 DEFLATE、PDF 的 1.4 生成器、XLSX 的 XML 解析都没有引入第三方库。代码量增加了，但依赖面、交叉编译难度与合规成本都降低了。

10.2 当前边界

原版 PSPP 定义 141 条命令，本版覆盖核心统计与数据变换命令，GLM、CTABLES、宏等长尾命令尚未实现；
支持 .sav 读写，未支持 .por 与 .spv；
图表覆盖 6 类，原版还有误差条图、ROC 曲线等；
未实现语法自动补全、输出树右键菜单等辅助交互；
输出引擎使用自研统一结果模型，没有复刻原版 .spv 中间格式与模板体系。

这些差距不影响核心分析链路，但决定了它当前更适合「日常统计分析 + 教学演示」，而不是替代原版的全部长尾场景。

十一、总结

格物统计的统计内核做了一件不太显眼但很关键的事：把上游依赖 GSL / GLib / GTK 的统计能力，重新落到一套只依赖 C++ 标准库、可在鸿蒙上交叉编译、可在宿主端直接测试的实现上。

它的技术骨架可以概括为：

math/ 提供矩统计、顺序统计、分布函数、矩阵、特征分解、回归、Logistic、因子与聚类；
stats/ 把这些数值能力组装成 15 类过程实现（对外 19 个分析入口），统一返回 AnalysisResult；
syntax/ 用词法 + 递归下降把命令翻译成同一组过程调用；
io/ 独立处理 CSV / XLSX / SAV，格式层与统计层解耦；
napi_init.cpp 用一组导出函数 + JSON 字符串把内核交给 ArkTS；
tests/ 在宿主端编译同一套源码，用 388 项断言与 numpy / scipy 对齐。

如果说上一篇讲的是「怎么在鸿蒙 PC 上做一个像样的统计软件」，这一篇讲的就是「怎么让它的数字可信」。前者是产品问题，后者是工程问题——两者都绕不过去，也都值得被认真对待。

开源地址（AtomGit 仓库）：https://atomgit.com/OpenHarmonyPCDeveloper/ohos_PSPP.git

人工智能6S服务平台

作为“人工智能6S店”的官方数字引擎，为AI开发者与企业提供一个覆盖软硬件全栈、一站式门户。

更多推荐

DBeaver-鸿蒙PC技术解析

把 DBeaver 搬到鸿蒙 PC，真正的工作量并不在「画界面」，而在三处结构性差异上。第一，复用边界要划在协议层，而不是 UI 层。数据库客户端最容易低估的是协议实现：认证插件、TLS、字符集、类型映射、错误码，每一项都是十几年打磨的结果。这里的取舍是——UI 与工作流 100% 重写，协议 100% 复用 C 库，桥接层只做统一形状，不做业务。结果是加一个新数据库只需要写一个 provider

人工智能6S服务平台

海蒂数据库管理工具-鸿蒙PC上架介绍

海蒂数据库管理工具把「原生重写」这条路走通了：界面与领域逻辑用 ArkTS / ArkUI 原生实现，协议层交叉编译成熟的 C 客户端库，NAPI 做桥，最终打成单一 HAP，装上即用——不存在独立进程或终端命令，是真正的原生桌面 GUI 应用。能复用的协议层绝不重写：libmariadb 这类 C 库在 musl 工具链下的坑（三元组、TRY_RUN、认证插件）虽多，但全部可以收敛进一个幂等脚本