统计/数值方法/统计软件的数值比较

介绍

统计计算需要额外的精度，并且容易出现截断和抵消误差等错误。这些错误是由于二进制表示和有限精度造成的，可能会导致不准确的结果。在本研究中，我们将讨论统计软件的准确性、可用于测量准确性的不同测试和方法，以及不同软件包的比较。

软件的准确性

准确性可以定义为结果的正确性。当使用统计软件包时，假设结果是正确的，以便对这些结果进行评论。然而，必须承认计算机有一些局限性。主要问题是，计算机系统通常提供的可用精度是有限的。很明显，统计软件无法提供超出这些限制的准确结果。但是，统计软件应该认识到其限制，并明确表明已达到这些限制。我们今天通常使用两种类型的精度

• 单精度
• 双精度

二进制表示和有限精度

二进制表示和有限精度会导致不准确的结果NIKHIL。计算机没有实数——我们用有限的近似值来表示它们。

示例：假设我们要用单精度浮点数表示 0.1。结果如下

0.1 = .00011001100110011001100110 = 0.99999964 (McCullough,1998)

很明显，我们只能以二进制形式近似 0.1。如果我们尝试减去两个仅在最低精度位不同的较大数字，则这个问题会增大。例如 100000.1-100000 = .09375

使用单精度浮点数，我们只能表示 24 个有效二进制位，换句话说，6-7 个有效十进制位。在双精度中，可以表示 53 个有效二进制位，相当于 15-17 个有效十进制位。二进制表示的局限性产生了五个不同的数值范围，会导致精度损失

• 负溢出
• 负下溢
• 零
• 正下溢
• 正溢出

溢出意味着值已变得太大，无法表示。下溢意味着操作数太小，太接近零，因此结果必须设置为零。单精度和双精度表示具有不同的范围。

二进制表示的结果

这些限制在不同情况下会导致不同的错误

• 抵消误差 是由于减去两个几乎相等的数字而产生的。
• 累积误差 是在一系列计算中累积的连续舍入误差之和。在这种类型的误差中，有可能只有结果的最右边的数字受到影响，或者结果没有一个准确的数字。
• 二进制表示和有限精度的另一个结果是，两个代数上等价的公式在数值上可能不等价。例如

$\sum _{n=1}^{10000}n^{-2}$ $\sum _{n=1}^{10000}(10001-n)^{-2}$

第一个公式按升序添加数字，而第二个公式按降序添加数字。在第一个公式中，最小的数字是在计算结束时才达到的，因此这些数字全部丢失，导致舍入误差。误差是第二个公式的 650 倍。（McCullough,1998）

• 截断误差 是由于二进制表示的限制而产生的近似误差。

示例

$\sin x=x-{\frac {x^{3}}{3!}}+{\frac {x^{5}}{5!}}-{\frac {x^{7}}{7!}}+\cdots$

sin(x) 的真实值与通过将有限个项相加而获得的结果之间的差值就是截断误差。（McCullough,1998）

• 算法误差 是另一个导致不准确的原因。计算某个量可能有多种方法，而这些不同的方法在准确性上可能会有所不同。例如，根据 Sawitzki (1994) 的说法，在单精度环境中，使用以下公式来计算方差 

$S^{2}=(1/(1-n)(\sum x_{i}^{2}-n{\bar {x}}^{2}))$

测量精度

由于计算机的限制，在计算统计值时会遇到一些问题。我们需要一个指标来衡量计算值的精度。该指标基于计算值（q）和真实值（c）之间的差异。一个常用的指标是对数相对误差（Log Relative Error，简称 LRE）（表示正确有效数字的位数）（McCullough，1998）

$LRE=-\log _{10}{\left[|q-c|/|c|\right]}$

规则

• q 应接近 c（少于两位数字）。如果它们不相等，则将 LRE 设置为零
• 如果 LRE 大于 c 中的数字位数，则将 LRE 设置为 c 中的数字位数。
• 如果 LRE 小于 1，则将其设置为零。

测试统计软件

在本部分，我们将讨论旨在衡量软件精度的两个不同测试：Wilkinson 测试（Wilkinson，1985）和 NIST StRD 基准测试。

Wilkinson 统计测试

Wilkinson 统计测试中使用的 Wilkinson 数据集“NASTY”是由 Leland Wilkinson（1985）创建的数据集。该数据集包含不同的变量，例如“Zero”，它只包含零，“Miss”，它包含所有缺失值，等等。NASTY 在它包含的值方面是一个合理的数据集。例如，“NASTY”中的“Big”的值小于美国人口，“Tiny”与工程中许多值相当。另一方面，“统计测试”的练习并非旨在合理。这些测试旨在检查统计计算中的一些特定问题。Wilkinson 的统计测试是一个入门级测试。

NIST StRD 基准测试

这些基准测试由美国国家标准与技术研究院设计，包含不同难度级别的数据集。目的是测试统计软件在统计学不同主题和不同难度级别方面的准确性。“统计参考数据集”项目的网页上有五组数据集

• 方差分析
• 线性回归
• 马尔科夫链蒙特卡洛
• 非线性回归
• 单变量汇总统计

在所有基准测试组中，有三种不同类型的数据集：低难度级别数据集、平均难度级别数据集和高难度级别数据集。通过使用这些数据集，我们将探索统计软件是否能够为一些统计计算提供 15 位有效数字的准确结果。

NIST 提供了 11 个数据集，其中有 6 个低难度级别数据集，2 个平均难度级别数据集，1 个高难度级别数据集。为每个数据集提供了 15 位有效数字的认证值，包括均值 (μ)、标准差 (σ) 和一阶自相关系数 (ρ)。

在方差分析数据集组中，有 11 个不同难度级别的数据集，其中有 4 个低难度级别数据集，4 个平均难度级别数据集和 3 个高难度级别数据集。对于每个数据集，都提供了 15 位有效数字的认证值，包括组间自由度、组内自由度、平方和、均方差、F 统计量、$R^{2}$ 和残差标准差。由于大多数认证值用于计算 F 统计量，因此只将其 LRE $\lambda _{F}$ 与统计软件包的结果进行比较。

为了测试统计软件的线性回归结果，NIST 提供了 11 个不同难度级别的数据集，其中有 2 个低难度级别数据集，2 个平均难度级别数据集和 7 个高难度级别数据集。对于每个数据集，我们都具有 15 位有效数字的认证值，包括系数估计、系数的标准误差、残差标准差、$R^{2}$、线性回归方差分析表，其中包括残差平方和。将最不准确的系数 $\lambda _{\beta }$、标准误差 $\lambda _{\sigma }$ 和残差平方和 $\lambda _{r}$ 进行比较。在非线性回归数据集组中，NIST 设计了 27 个数据集，其中有 8 个低难度级别数据集，11 个平均难度级别数据集和 8 个高难度级别数据集。对于每个数据集，我们都具有 NIST 提供的 11 位有效数字的认证值，包括系数估计、系数的标准误差、残差平方和、残差标准差和自由度。

在非线性回归计算的情况下，我们应用曲线拟合方法。在该方法中，我们需要起始值来初始化等式中的每个变量。然后，我们生成曲线并计算收敛标准（例如，平方和）。然后，我们调整变量以使曲线更接近数据点。有几种算法可以用来调整变量

• Marquardt 和 Levenberg 方法
• 线性下降方法
• Gauss-Newton 方法

反复应用其中一种方法，直到收敛标准的差异小于收敛容差。

NIST 还提供了两组起始值：起始值 I（远离解决方案的值）、起始值 II（接近解决方案的值）。使用起始值 II 作为初始值使更容易获得准确的解决方案。因此，将优先考虑起始值 I 解决方案。

其他重要设置如下

• 收敛容差（例如，1E-6）
• 求解方法（例如，Gauss Newton 或 Levenberg Marquardt）
• 收敛准则（例如，残差平方和 (RSS) 或参数差异最大值的平方）

我们也可以在数值导数和解析导数之间进行选择。

测试示例

测试软件包：SAS、SPSS 和 S-Plus

在本部分中，我们将讨论 M.D. McCullough 应用的三种统计软件包的测试结果。在 McCullough 的工作中，SAS 6.12、SPSS 7.5 和 S-Plus 4.0 进行了测试和比较，以了解 NIST 提供的认证 LRE 值。比较将根据以下部分进行

• 单变量统计
• 方差分析
• 线性回归
• 非线性回归

单变量统计

表 1：SAS 单变量统计结果（McCullough，1998）

在 SAS 中计算的所有值似乎都或多或少地准确。对于数据集 NumAcc1，由于观察次数不足，无法计算 p 值。计算数据集 NumAcc3（平均难度）和 NumAcc 4（高难度）的标准差似乎给 SAS 带来了压力。

表 2：SPSS 单变量统计结果（McCullough，1998）

对于均值和标准差计算的所有值似乎都或多或少地准确。对于数据集 NumAcc1，由于观察次数不足，无法计算 p 值。计算数据集 NumAcc3 和 -4 的标准差似乎也给 SPSS 带来了压力。对于 p 值，SPSS 只用 3 位小数表示结果，这低估了第一个 p 值的精度，高估了最后一个 p 值的精度。

表 3：S-Plus 单变量统计结果（McCullough，1998）

对于均值和标准差计算的所有值似乎都或多或少地准确。S-Plus 在计算数据集 NumAcc3 和 -4 的标准差时也存在问题。S-Plus 在计算 p 值方面表现不佳。

方差分析

表 4：SAS 方差分析结果（McCullough，1998）

结果

• SAS 只能解决低级别难度的方差分析问题。
• 对于平均或更高难度的的数据集，F 统计量可以计算出来，但性能很差，精度为 0 位。
• SPSS 可以像 SAS 一样，显示对于低级别难度的数据集的准确结果。
• SPSS 在计算方差分析方面的表现不佳。
• 对于数据集“AtmWtAg”，SPSS 没有显示 F 统计量，这似乎比显示 0 精度结果更合理。
• S-Plus 比其他软件包更能处理方差分析问题。
• 即使对于高难度的数据集，S-Plus 也可以显示比其他软件包更准确的结果。但对于高难度数据集的结果仍然不够准确。
• S-Plus 可以解决中等难度的问题，其精度足够高。

线性回归

表 5：SAS 线性回归结果（McCullough，1998）

SAS 无法为 Filip 数据集提供解决方案，该数据集是十阶多项式。除了 Filip 之外，SAS 显示的结果或多或少地准确。但对于更高难度的 datasets，性能似乎下降，尤其是在计算系数方面。

表 6：SPSS 线性回归结果（McCullough，1998）

SPSS 也存在“Filip”问题，这是一个 10 阶多项式。许多软件包无法计算其值。与 SAS 一样，SPSS 在处理高难度 datasets 时，精度较低。

表 7：S-Plus 线性回归结果（McCullough，1998）

S-Plus 是唯一一个为“Filip”数据集提供结果的软件包。Filip 结果的精度似乎不是差，而是中等。即使对于更高难度的 datasets，S-Plus 也能比其他软件包计算出更准确的结果。只有数据集“Wrampler4”和“-5”的系数精度低于平均水平。

非线性回归

表 8：SAS 非线性回归结果（McCullough，1998）

对于非线性回归，每个软件包测试了两种设置组合，因为不同的设置会影响结果。从表格中可以看出，SAS 中首选的组合比默认组合产生的结果更好。在此表格中，使用默认组合产生的结果用括号表示。由于 NIST 提供了 11 位的认证值，因此我们寻找的是 11 的 LRE 值。

首选组合

• 方法：高斯-牛顿
• 准则：PARAM
• 容差：1E-6

表 9：SPSS 非线性回归结果（McCullough，1998）

在 SPSS 中，首选组合也比默认选项表现更好。所有问题都在初始值“开始 I”下解决，而在 SAS 中，更高难度的 datasets 则使用开始 II 值解决。

首选组合

• 方法：列文伯格-马夸特
• 准则：PARAM
• 容差：1E-12

表 10：S-Plus 非线性回归结果（McCullough，1998）

从表格中可以看出，首选组合在 S-Plus 中也比默认组合更好。除了“MGH10”之外，所有问题都在初始值“开始 I”下解决。可以说，S-Plus 在计算非线性回归方面表现比其他软件更好。

首选组合

• 方法：高斯-牛顿
• 准则：RSS
• 容差：1E-6

比较结果

所有软件包都在单变量统计中为均值和标准差提供了准确的结果。测试的统计软件包之间没有太大差异。在方差分析计算中，SAS 和 SPSS 无法通过中等难度的问题，而 S-Plus 比其他软件包提供了更准确的结果。但对于高难度数据集，它也产生了较差的结果。关于线性回归问题，所有软件包似乎都可靠。如果我们检查所有软件包的结果，可以说，非线性回归结果计算的成功很大程度上取决于所选择的选项。

其他重要结果如下

• S-Plus 一次从开始 II 开始解决。
• SPSS 从未使用开始 II 作为初始值，但有一次产生了 0 精度位。
• SAS 使用开始 II 三次，每次都产生了 0 精度位。

不同版本的 SPSS 比较

在本部分中，我们将比较 SPSS 的旧版本和新版本，以查看旧版本中的问题是否在新版本中得到解决。在本部分中，我们将 SPSS 版本 7.5 与 SPSS 版本 12.0 进行了比较。版本 7.5 的 LRE 值来自 B.D. McCullough 的文章（参见参考文献）。我们还将这些测试应用于版本 12.0，并计算了相应的 LRE 值。我们从每个难度组中选择一个数据集，并在版本 12.0 中应用了单变量统计、方差分析和线性回归。数据集的来源是 NIST 统计参考数据集档案。然后，我们使用 NIST 提供的认证值计算每个数据集的 LRE 值，以便比较两个版本的 SPSS。

单变量统计

难度：低

我们的第一个数据集是 PiDigits，它具有较低级别的难度，由 NIST 设计，用于检测计算单变量统计值的缺陷。

PiDigits 的认证值为：

• 样本均值：4.53480000000000
• 样本标准差：2.86733906028871

文件：Version12pd.jpg

表 11：SPSS 版本 12.0 对 PiDigits 数据集的结果

如表 13 所示，SPSS 12.0 的结果与 NIST 提供的认证值相符。因此，我们针对均值和标准差的 LRE 为 $\lambda _{\mu }$: 15，$\lambda _{\delta }$: 15。在 7.5 版本中，LRE 值为 $\lambda _{\mu }$: 14.7，$\lambda _{\delta }$: 15。（McCullough，1998）

难度：中等

第二个数据集是 NumAcc3，其难度来自 NIST 数据集的单变量统计平均值。NumAcc3 的认证值为：

• 样本均值：1000000.2
• 样本标准差：0.1

文件：Version12na3.jpg

表 12：SPSS 12.0 版本针对 NumAcc3 数据集的结果

从表 14 中我们可以看到，计算出的均值与 NIST 的认证值相同。因此，我们针对均值的 LRE 为 $\lambda _{\mu }$: 15。然而，标准差值与认证值不同。因此，针对标准差的 LRE 计算如下：

$\lambda _{\delta }$ : -log10 |0,10000000003464-0,1|/|0,1| = 9.5

SPSS v 7.5 的 LRE 为 $\lambda _{\mu }$: 15，$\lambda _{\delta }$: 9.5。（McCullough，1998）

难度：高

单变量统计中的最后一个数据集是 NumAcc4，其难度较高。NumAcc4 的认证值为：

• 样本均值：10000000.2
• 样本标准差：0.1

文件：Version12na4.jpg

表 13：SPSS 12.0 版本针对 NumAcc4 数据集的结果

对于这个数据集，我们也没有遇到任何关于计算出的均值的问题。因此，LRE 为 $\lambda _{\mu }$: 15。然而，标准差值与认证值不匹配。因此，我们应该按如下方式计算标准差的 LRE：

$\lambda _{\delta }$ : -log10 |0,10000000056078-0,1|/|0,1| = 8.3

SPSS v 7.5 的 LRE 为 $\lambda _{\mu }$: 15，$\lambda _{\delta }$ : 8.3 (McCullough，1998)

对于测试的这一部分，我们可以说 SPSS 的两个版本之间没有区别。对于中等和高难度数据集，提供的标准差结果仍然具有中等准确性。

方差分析

难度：低

我们用于测试 SPSS 12.0 针对较低难度级别问题的 数据集是 SiRstv。SiRstv 的认证 F 统计量为 1.18046237440255E+00

文件：Version12ss.jpg

表 14：SPSS 12.0 版本针对 SiRstv 数据集的结果

• LRE：$\lambda _{F}$: -log10 | 1,18046237440224- 1,18046237440255|/ |1,18046237440255| = 12,58
• SPSS v 7.5 的 LRE：$\lambda _{F}$ : 9,6 (McCullough，1998)

难度：中等

我们用于中等难度级别问题的 数据集是 AtmWtAg。AtmWtAg 的认证 F 统计量值为 1.59467335677930E+01。

文件：Version12awa.jpg

表 15：SPSS 12.0 版本针对 AtmWtAg 数据集的结果

• LREs : $\lambda _{F}$ : -log10 | 15,9467336134506- 15,9467335677930|/| 15,9467335677930| = 8,5
• SPSS v 7.5 的 LREs : $\lambda _{F}$: 丢失

难度：高

我们使用 SmnLsg07 数据集来测试高难度问题。SmnLsg07 的认证 F 值为 2.10000000000000E+01

文件:Version12sl07.jpg

表 16: SPSS 12.0 版本对数据集 SmnLsg07 的结果

• LREs : $\lambda _{F}$ : -log10 | 21,0381922055595 - 21|/| 21| = 2,7
• SPSS v 7.5 的 LREs : $\lambda _{F}$: 0

在 12.0 版本中计算的方差分析结果比 7.5 版本中计算的结果更好。但是准确度仍然太低。

线性回归

难度：低

我们的低难度数据集是 Norris，用于线性回归。Norris 的认证值为：

• 样本残差平方和 : 26.6173985294224

• 

  表 17: Norris 的系数估计（www.itl.nist.gov）

• LREs : $\lambda _{r}$ : 9,9 $\lambda _{\beta }$ : 12,3 $\lambda _{\sigma }$ : 10,2
• SPSS v 7.5 的 LREs : $\lambda _{r}$: 9,9 , $\lambda _{\beta }$ : 12,3 , $\lambda _{\sigma }$ : 10,2 (McCullough, 1998)

文件:Version12nor.jpg

表 18: SPSS 12.0 版本对数据集 Norris 的结果

难度：中等

我们使用 NoInt1 数据集来测试平均难度数据集的性能。回归模型如下：

y = B1*x + e

NoInt1 的认证值 

• 样本残差平方和 : 127,272727272727
• 系数估计 : 2.07438016528926, 标准差 : 0.16528925619834E-0(www.itl.nist.gov)
• LREs: $\lambda _{r}$:12,8 $\lambda _{\beta }$: 15 $\lambda _{\sigma }$: 12,9
• SPSS v. 7.5 的 LREs : $\lambda _{r}$: 12,8 , $\lambda _{\beta }$: 14,7 , $\lambda _{\sigma }$: 12,5 (McCullough, 1998)

文件:Version12nint.jpg

表 19: SPSS 12.0 版本对数据集 NoInt1 的结果

难度：高

我们的高难度数据集是 NIST 设计的 Longley。

• 模型: y =B0+B1*x1 + B2*x2 + B3*x3 + B4*x4 + B5*x5 + B6*x6 +e
• LREs 
  • $\lambda _{r}$: -log10 |836424,055505842-836424,055505915|/ |836424,055505915| = 13,1
  • $\lambda _{\beta }$ : 15
  • $\lambda _{\sigma }$ : -log10 | 0,16528925619836E-01 – 0,16528925619834E-01|/ |0,16528925619834E-01| = 12,9
• SPSS v. 7.5 的 LREs : $\lambda _{r}$: 12,8 , $\lambda _{\beta }$ : 14,7 , $\lambda _{\sigma }$ : 12,5 (McCullough, 1998)

从计算的 LREs 可以得出结论，两个版本的线性回归结果之间没有太大差异。

结论

通过应用这些测试，我们试图找出软件是否可靠并提供准确的结果。然而，根据结果，我们可以说不同的软件包对同一个问题的结果不同，这可能会导致我们对统计研究问题的错误解读。

具体来说，我们可以看到 SAS、SPSS 和 S-Plus 在解决线性回归问题方面优于 ANOVA 问题。所有这三者在 F 统计量计算方面都提供较差的结果。

从比较 SPSS 的两个不同版本的结果来看，我们可以得出结论，考虑到版本号之间的差异，SPSS v.12 和 v.7.5 提供的结果的准确性之间的差异并不大。另一方面，SPSS v.12 在处理 ANOVA 问题方面比旧版本好得多。然而，它在更高难度的问题上仍然存在问题。

参考文献

• McCullough, B.D. 1998, 'Assessing The Reliability of Ststistical Software: Part I',The American Statistician, Vol.52, No.4, pp.358-366.
• McCullough, B.D. 1999, 'Assessing The Reliability of Ststistical Software: Part II', The American Statistician, Vol.53, No.2, pp.149-159
• Sawitzki, G. 1994, 'Testing Numerical Reliability of Data Analysis Systems', Computational Statistics & Data Analysis, Vol.18, No.2, pp.269-286
• Wilkinson, L. 1993, 'Practical Guidelines for Testing Statistical Software' in 25th Conference on Statistical Computing at Schloss Reisenburg, ed. P. Dirschedl& R. Ostermnann, Physica Verlag
• National Institute of Standards and Technology. (1 September 2000). The Statistical Reference Datasets: Archives, [Online], Available from: <http://www.itl.nist.gov/div898/strd/general/dataarchive.html> [10 November 2005].