Technical Article

Delphi 中的 Excel 统计函数:NORM、CHISQ、BETA

在单元格中输入 =NORM.DIST(115,100,15,TRUE),Excel 会直接返回 0.8413447。该调用读起来像是个查找。其实不然。在这一个数字背后,是累积正态分布(一个没有闭合形式的积分),在 CHISQ.INV.RTBETA.DIST 背后,是严谨的库必须评估而非手动近似的特殊函数。声称与 Excel 兼容的电子表格组件必须重现这些值,直至 Excel 显示的最后一位数字,这意味着必须重现数值方法,而不只是函数名称。

HotXLS 实现了五十多个此类统计函数,使它们正确的工作几乎完全在公式栏中不可见。这是关于引擎如何计算它们的旅程:共享的特殊函数核心、保持算术稳定的分支决策,以及一个长期隐藏在尾部的逆正态分布漏洞,因为常见情况从未触及出错的代码行。

一个工作表调用,其后是五十种分布

这些函数涵盖了统计工作簿所涉及的大家族。有正态分布族(NORM.DISTNORM.S.DIST 及其逆函数);伽马和卡方分布族(GAMMA.DISTCHISQ.DISTCHISQ.DIST.RTCHISQ.INV.RT);贝塔分布族(BETA.DISTBETA.INV);抽样分布(T.DISTT.DIST.2TF.DISTF.INV);离散对(BINOM.DISTPOISSON.DIST);以及推断辅助函数(例如 CONFIDENCE.TCONFIDENCE.NORM)。从调用者的角度来看,每个都是单一的公式。您在单元格中设置输入,要求工作簿进行评估并读取结果。

var
  wb: IXLSWorkbook;
  sh: IXLSWorksheet;
begin
  wb := TXLSWorkbook.Create;
  sh := wb.Sheets.Add;
  sh.Range['A1', 'A1'].Value := 115;   // observation
  sh.Range['A2', 'A2'].Value := 100;   // mean
  sh.Range['A3', 'A3'].Value := 15;    // standard deviation

  // The XLS formula parser uses ';' as the argument separator.
  Writeln(wb.Calculate('=NORM.DIST(A1;A2;A3;TRUE())'));   // 0.8413447
  Writeln(wb.Calculate('=CHISQ.INV.RT(0.05;10)'));        // 18.3070381
  Writeln(wb.Calculate('=BETA.DIST(0.5;2;3;TRUE())'));    // 0.6875
end;

工作簿上的 Calculate 方法编译并针对活动工作表评估临时公式,并返回 Variant。有一个细节会在初次尝试时让人跌倒:Calculate 背后的公式解析器采用分号作为其参数分隔符,因此是 =SUM(A1;B1),而不是 =SUM(A1,B1)。存储的单元格公式保留 Excel 标准的逗号。相同的评估器调度下面的每个统计函数,因此一旦其中一个在 Calculate 中工作,其余的就遵循相同的路径。

其他一切赖以构建的两个函数

此集合中的大多数累积分布不是通过对其自身的定义进行求和或积分来计算的。它们是从两个特殊函数计算出来的:正则化下不完全伽马(记作 P(a, x))和正则化不完全贝塔(记作 Ix(a, b))。在内部,这些是调度器依赖的辅助函数,且逻辑链很短。卡方累积分布函数(CDF)是形状为 df/2 且尺度为 2 的伽马 CDF。伽马 CDF 直接就是 P(a, x)。t、F 和二项累积函数都是正则化不完全贝塔在合适参数下的值。泊松的 CDF 是上不完全伽马 Q。很好地实现伽马和贝塔函数,十几种分布就会免费继承它们的精确性。

“正则化”这个词是全部重点。原始的不完全伽马呈阶乘式增长,而原始贝塔积分在答案出现之前很久就可能下溢或溢出。正则化形式除以完全伽马或贝塔,因此它们完全存在于零到一的区间内,这正是概率所占的范围。这种归一化使得同一个例程可以服务于两个自由度的卡方分布以及两百个自由度的卡方分布,而中间项不会超出双精度浮点数的范围。它还解释了为什么您不通过将长尾的密度项加在一起来计算 CDF:每一项都带有其自身的四舍五入误差,误差随着级数的运行而累积,而正则化特殊函数通过评估快速收敛的级数或连分数,完全避开了求和。

对角线下方用级数,对角线上方用连分数

不完全伽马例程在计算任何内容之前会做出一个决定:它将 x 与 a + 1 进行比较。该边界不是任意的。当 x 相对于 a 较小时,P(a, x) 的幂级数展开收敛很快,而当 x 较大时,收敛很慢,最终会毫无用处。连分数具有相反的特征。因此,引擎对低于 a + 1 的 x 使用幂级数,对等于或高于 a + 1 的 x 使用 Lentz 连分数,每个分支只被要求做它擅长的工作。

连分数需要一个保护。Lentz 的方法通过携带运行的分子和分母并在每一步反转分母来工作,如果其中任何一个接近零,反转就会爆炸。解决方法是微小的下限:每当中间项的大小降至大约 1e-30 以下时,它就会被钳制到 1e-30,这使递推保持有限而不会干扰收敛值。出于相同的原因,相同的夹钳也出现在不完全贝塔的连分数中。这是一个做着承重工作的小常数,是稳定评估与除以与零无法区分的值之间的区别。

上尾 Q(a, x) 只是 1 减去 P(a, x),这就是泊松累积分支的计算方式:均值为 λ 的至多 k 次事件的概率是 Q(k + 1, λ)。再次,将其路由通过上不完全伽马而不是对 k + 1 个泊松项求和,是选择评估一个收敛表达式而不是累积许多小表达式。

没有阶乘溢出的离散概率质量

离散分布引发了不同的危险。二项概率质量涉及二项式系数,52 选 26 的系数是一个巨大的整数。直接形式化它,分子会在进行除法将其带回合理概率之前溢出双精度浮点数。引擎从不直接构建它。它通过对数伽马函数在对数空间中计算阶乘,对数进行加减,折入成功和失败概率的对数,并最终在最后进行一次指数运算。

// Binomial probability mass, evaluated entirely in log space.
// LnGammaF(n+1) is ln(n!); the three log-factorials form ln(C(n,k)),
// and the whole exponent is built before a single Exp call.
//   ln P(X=k) = ln(n!) - ln(k!) - ln((n-k)!) + k*ln(p) + (n-k)*ln(1-p)
result := Exp(LnGammaF(nt + 1) - LnGammaF(kk + 1) - LnGammaF(nt - kk + 1)
  + kk * Ln(pp) + (nt - kk) * Ln(1 - pp));

对数伽马函数本身是 Lanczos 近似,在整个正轴上精确且评估开销低。因为直到最终的 Exp 之前每个大数量都保留为其对数,因此该例程物化出的最大数字是概率本身(最大为 1)。泊松质量函数遵循相同的配方,单个对数伽马项代替分母中的阶乘。闭合形式在边缘进行特例处理(其中 p 还是为零或一),因此代码绝不会调用 Ln(0)。HotXLS 为 BINOM.DIST(5,10,0.5,FALSE) 返回 0.2460938,为累积的 POISSON.DIST(2,2,TRUE) 返回 0.6766764,在显示的所有数字上与 Excel 相匹配。

通过包围正向曲线求逆

逆分布函数提出相反的问题:给定概率,找到其 CDF 等于它的值。此集合中只有一个逆分布具有快速的直接公式。NORM.S.INV(逆标准正态)使用 Acklam 有理近似(一对多项式比率,在整个范围内大约精确到双精度浮点数的精度),分为中心区域和两个尾部。这是没有迭代的闭合形式评估。

其他逆函数没有此类公式,因此引擎对其进行数值求逆。它用从分布支撑集中选择的上限和下限包围答案,然后进行二分搜索:在中点处评估正向 CDF,移动保持目标概率被包围的任何边界,并重复直到区间变窄。对于伽马和卡方逆分布,包围从零开始,以及根据形状和尺度构建的慷慨上限估计,如果概率尚未被包围,则将上限倍增。t 逆分布包围向外变宽的对称边界;F 逆分布在非负区间上进行二分。代价是每次调用要进行几十次 CDF 评估(在电子表格速度下不可见),好处是每个逆分布与它求逆的正向函数精确一样精确。这就是为什么诸如 CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) 的往返过程极度精确地返回 0.7 的原因。

隐藏在尾部的常用对数

这里是值得讲述的漏洞,因为它是那种能够存活很长时间的漏洞。Acklam 逆正态例程有三个分支。无论何时当概率位于大约 0.025 和 0.975 之间时,使用宽的中心分支,通过没有任何对数的多项式比率运行输入。对于非常小或非常大的概率,两个尾部分支各首先对输入取对数,因为尾部的行为就像负的 p 的自然对数的平方根。

尾部分支的早期版本在应该使用自然对数的地方使用了常用对数(以 10 为底)。两者相差大约 2.30 的常数因子,因此尾部结果出现了持续且相当大的偏差。然而,该函数在每一次日常检查中看起来都很好,因为日常检查都处于中间区域。NORM.S.INV(0.5) 为零,NORM.S.INV(0.975) 是教科书式的 1.959964,这两个都是通过完全不调用对数的中心多项式运行的。错误只有在概率越入尾部时才会出现(例如 NORM.S.INV(0.001),它必须返回 -3.0902323,结果却由于自然对数与常用对数的比例关系而发生了偏移)。任何在尾部依赖于逆正态分布的函数(包括置信区间辅助函数)都继承了相同的倾斜。教训是平凡而昂贵的:具有分支结构的函数在每个分支内部都需要测试点,因为正确的常用路径会愉快地掩盖损坏的罕见路径。修复工作是将常用对数更改为自然对数这单一代币的更改,尾部值瞬间对齐了 Excel。

x 的符号决定 t 分布的尾部

学生 t 累积函数带有一个很容易弄反的微妙细节。它的值来自在 df / (df + x²) 处评估的正则化不完全贝塔,但该贝塔值是超出 x 大小的尾部概率,而不是直至 x 的累积概率。t 分布的对称形状意味着转换取决于 x 落入零的哪一侧。

// Student t CDF. ib is the regularized incomplete beta at df/(df+x*x),
// which measures the symmetric tail. The cumulative value depends on
// the sign of x; returning ib unconverted gives the wrong tail.
ib := BetaIF(df / 2, 0.5, df / (df + x * x));
if x > 0 then
  result := 1 - 0.5 * ib        // above the mean: one minus half the tail
else if x < 0 then
  result := 0.5 * ib            // below the mean: half the tail
else
  result := 0.5;                // exactly at the mean

对于零之上的 x,累积概率是一减去对称尾部的一半;对于零以下的 x,它是该尾部的一半;在零处正好是二分之一。直接返回贝塔值,您将报告分布的错误一侧,对于任何非零 x,这与曲线的整个主体偏差了。右尾和双尾变体构建在相同的分支上,这就是为什么 T.DIST.2T(1,1) 返回为 0.5 且 T.DIST(1,1,TRUE) 为 0.75 的原因,而逆函数 T.INV 针对此纠正后的 CDF 进行二分,使得往返闭合。

在单元格中所有这些都是不可见的,这正是预期的结果。您编写公式并读取与 Excel 一致的数字。如果您正在使用自己的逻辑扩展引擎,注册函数的方法在我们关于公式引擎和自定义函数的演练中进行介绍,公式跨越表格和命名范围访问的方式在关于定义名称和跨表格公式的文章中进行介绍。所有这些都提供在适用于 Delphi 和 C++Builder 的 HotXLS 电子表格组件内部,与本博客其他地方介绍的读取、写入、图表绘制和格式化 API 一起交付。