项目
10/17/2016

2016 年 10 月

第 31 卷，第 10 期

此文章由机器翻译。

测试运行 - 使用 C# 执行 ANOVA

James McCaffrey 方差分析 (ANOVA) 是一项经典统计技术，用于推断三个或多个组的平均值是否均相等（在仅有示例数据的情况下）。例如，假设在某大学开设三门不同的计算机导论课程。每门课由同一位教师教授，但使用不同的课本和不同的教学理念。你想要知道学生的表现是否一样。

举行一个 1 到 15 人规模的评估计算机能力的考试，但是因为考试成本和时间成本都很高，所以每门课程中可以只对 6 个随机选择的学生进行测试。管理考试并对选样执行 ANOVA 以推断所有三门课的平均值是否相同。

如果你还不太熟悉 ANOVA，可能会对该技术的名称略感困惑，因为它的目标是分析数据集的平均值。之所以命名为 ANOVA 是因为该技术在后台分析方差来推断平均值。

若要了解什么是 ANOVA 以及本文要讨论的内容，一个很好的方法是查看图 1 中的演示程序。该演示为三个组设置硬编码的分数。请注意，在组 1 中仅有 4 个分数，在组 3 中仅有 5 个分数 - 这种示例大小不相等的情况是很常见的，因为测试主体可能会退出，或数据可能会丢失或损坏。

图 1 使用 C# 执行 ANOVA 的演示程序

执行 ANOVA 主要有两个步骤。在第一步中，使用示例数据计算出 F-统计量的值和一对名为“自由度”(df) 的值。在第二步中，使用 F 和 df 值确定所有人数平均值相等的可能性（即 p-值）。第一步相对简单。第二步很难。

在该演示中，F 值为 15.884。通常情况下，F 值越大，所有人数平均值相等的可能性越小。稍后我将解释为什么 df = (2, 12)。使用 F 和 df 值计算出的 p-值为 0.000425。该值非常小，因此你可以得出结论，人数平均值很很可能不相等。这时，你可以执行额外的统计测试以确定哪个人数平均值与其他不同。对于演示数据，看起来组 1（示例平均值 = 4.50）比组 2（平均值 = 9.67）和组 3（平均值 = 10.60）更糟糕。

演示程序

为了创建演示程序，我启动 Visual Studio，单击“文件 | 新建 | 项目”并选择“C# 控制台应用程序”选项。该演示程序对 .NET 的依赖程度并不明显，因此，任何 Visual Studio 版本都可以正常运行。在“解决方案资源管理器”窗口加载模板代码后，我右键单击了文件 Program.cs 并将其重新命名为 AnovaProgram.cs 且允许 Visual Studio 自动为我重新命名类 Program。

在编辑器窗口的顶部，我使用语句删除了所有不必要的内容，仅留下引用顶级 System 命名空间的内容。程序的整体结构如图 2 中所示。该演示程序过长，无法全部展示，但是本文随附了完整的演示源代码以供下载。

图 2 演示程序结构

using System;
namespace Anova
{
  class AnovaProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("Begin ANOVA using C# demo");
      // Set up sample data
      // Use data to calculate F and df
      // Use F and df to calculate p-value
      Console.WriteLine("End ANOVA demo");
    }
    static double Fstat(double[][] data, out int[] df) { . . }
    static double LogGamma(double z) { . . }
    static double BetaIncCf(double a, double b, double x) { . . }
    static double BetaInc(double a, double b, double x) { . . }
    static double PF(double a, double b, double x) { . . }
    static double QF(double a, double b, double x) { . . }
    static void ShowData(double[][] data, string[] colNames) { . . }
  }
}

静态方法 Fstat 基于存储于数组中的数组对象的数据进行计算并返回 F-统计量。此方法也会在数组输出参数中进行计算并返回两个 df 值。函数 ShowData 只是显示示例平均值的小 helper 函数。

剩余的五种方法均用于计算 p-值。方法 QF 是主要方法。它调用方法 PF，PF 随之调用方法 BetaInc，而 BetaInc 随之调用方法 BetaIncCf 和 BetaIncCf。

在一些初设 WriteLine 消息后，Main 方法设置并显示示例数据：

double[][] data = new double[3][]; // 3 groups
data[0] = new double[] { 3, 4, 6, 5 };
data[1] = new double[] { 8, 12, 9, 11, 10, 8 };
data[2] = new double[] { 13, 9, 11, 8, 12 };
string[] colNames = new string[] {  "Group1", "Group2", "Group3" };
ShowData(data, colNames);

在非演示方案中，你的数据可能存储在文本文件中，你将编写 helper 函数以读取数据，并将数据加载到数组中的数组。

F-统计量和 df 按以下方式进行计算：

int[] df = null;
double F = Fstat(data, out df);

对于 ANOVA，数据集的 df 为一对值。第一个值为 K - 1，其中 K 是组数，而第二个值为 N - K，其中 N 是示例值的总数。因此，对于演示数据，df = (K-1, N-K) = (3-1, 15-3) = (2, 12)。

会计算 p-值并显示如下：

double pValue = QF(df[0], df[1], F);
Console.Write("p-value = ");

简言之，在执行 ANOVA 时，调用语句非常简单。但是在后台需执行大量工作。

计算 F-统计量

计算 F-统计量的值有几个子步骤。假设示例数据值是来自演示程序的数据值：

Group1: 3.00, 4.00, 6.00, 5.00
Group2: 8.00, 12.00, 9.00, 11.00, 10.00, 8.00
Group3: 13.00, 9.00, 11.00, 8.00, 12.00

第一个子步骤是计算每个组的平均值以及所有示例值的总体平均值。对于演示数据：

means[0] = (3.0 + 4.0 + 6.0 + 5.0) / 4 = 4.50
means[1] = (8.0 + 12.0 + 9.0 + 11.0 + 10.0 + 8.0) / 6 = 9.67
means[2] = (13.0 + 9.0 + 11.0 + 8.0 + 12.0) / 5 = 10.60
gMean = (3.0 + 4.0 + . . . + 12.0) / 15 = 8.60

可以从下面开始定义方法 Fstat：

static double Fstat(double[][] data, out int[] df)
{
  int K = data.Length; // Number groups
  int[] n = new int[K]; // Number items each group
  int N = 0; // total number data points
  for (int i = 0; i < K; ++i) {
    n[i] = data[i].Length;
    N += data[i].Length;
  }
...

此时，本地数组 n 表示每个组中的值数，K 表示组数，而 N 是所有组中的值的总数。接下来，会将组平均值计算为数组命名的平均值，而总体的总计平均值则计算为变量 gMean：

double[] means = new double[K];
double gMean = 0.0;
for (int i = 0; i < K; ++i) {
  for (int j = 0; j < data[i].Length; ++j) {
    means[i] += data[i][j];
    gMean += data[i][j];
  }
  means[i] /= n[i];
}
gMean /= N;

下一个子步骤是计算“组间的平方和”(SSb) 和“组间的均方”(MSb)。SSb 是每组平均值和整体平均值间的方差加权和。MSb = SSb / (K-1)，其中 K 是组数。对于演示数据：

SSb = (4 * (4.50 - 8.60)^2) +
 (6 * (9.67 - 8.60)^2) +
 (5 * (10.60 - 8.60)^2) = 94.07
MSb = 94.07 / (3-1) = 47.03

计算 SSb 和 MSb 的代码为：

double SSb = 0.0;
for (int i = 0; i < K; ++i)
  SSb += n[i] * (means[i] - gMean) * (means[i] - gMean);
double MSb = SSb / (K - 1);

下一个子步骤是计算“组内的平方和”(SSw) 和“组内的均方”(MSw)。SSw 是每个示例值及其组平均值之间的方差和。MSw = SSw / (N-K)。对于演示数据：

SSw = (3.0 - 4.50)^2 + . . + (8.0 - 9.67)^2 +
 . . + (12.0 - 10.60)^2 = 35.53
MSw = 35.53 / (15-3) = 2.96

计算 SSw 和 MSw 的代码为：

double SSw = 0.0;
for (int i = 0; i < K; ++i)
  for (int j = 0; j < data[i].Length; ++j)
    SSw += (data[i][j] - means[i]) * (data[i][j] - means[i]);
double MSw = SSw / (N - K);

最后一个子步骤是计算两个 df 值和 F-统计量。两个 df 值为 K - 1、 N - K。并且 F = MSb / MSw。对于演示数据：

df = (K-1, N-K) = (3-1, 15-3) = (2, 12)
F = 47.03 / 2.96 = 15.88.

计算 df 和 F 的演示代码为：

...
  df = new int[2];
  df[0] = K - 1;
  df[1] = N - K;
  double F = MSb / MSw;
  return F;
} // Fstat

我想你不会否认：如果知道数学等式，通过一组数据计算 F-统计量和 df 值是机械的，且相对简单。

计算 p-值

基于生成 F 和 df 的示例数据将 F-统计量和 df 值转换为指出所有人数平均值均相等有多大可能性的 p-值理论上很简单，但在实际操作中很困难。我将尽可能简要介绍，而剩余的大量详细信息则需要做大量的额外说明。请看图 3 中的曲线图。

图 3 通过 F-统计量和 df 值计算 p-值

每对可能的 df 值都将确定一个曲线图，称为 F-分布。基于 df 的值，F-分布的形状可能具有很大的差别。图 3 中的曲线图显示 df = (4, 12) 的 F-分布。我从演示数据中使用了 df = (4, 12) 而非 df = (2, 12)，因为 df = (2, 12) F-分布的形状是非典型的。

任何 F-分布下的总区域恰好是 1.0。如果你知道 F-统计量的值，则 p-值是从 F 到正无穷大的 F-分布下的区域。有点困惑的是，从 0 到 F-统计量的 F-分布下的区域通常被称为 PF，而从 F-统计量到正无穷大的 F-分布下的区域通常被称为 QF（表示 p-值）。因为分布下的总区域是 1，所以 PF + QF = 1。事实证明计算 PF 比计算 QF 容易一些，因此找到 p-值 (QF)，通常计算出 PF，然后用 1 减去 PF 可得到 QF。

计算 PF 难度很大，但是幸运的是，神奇的估计等式人们已熟知几十年。这些数学等式，以及数百种其他数学等式，可以在一本有名的参考资料中找到，即由 M. Abramowitz 和 I.A.Stegun 所著的“数学函数手册”。该书在科学编程人员间经常被简称为“A&S”。每个 A&S 等式均有一个 ID 号。

在演示中，方法 PF 的确只是方法 BetaInc 的包装：

static double PF(double a, double b, double x)
{
  double z = (a * x) / (a * x + b);
  return BetaInc(a / 2, b / 2, z);
}

方法 BetaInc 的名称代表“不完整的 Beta”。方法 BetaInc 使用 A&S 等式 6.6.2 和 26.5.8。这些等式调用 LogGamma 函数和 BetaIncCf 函数。LogGamma 函数的解释和执行极其困难。简单地说，数学 Gamma 函数将阶乘的概念扩展到实值数。与阶乘一样，Gamma 函数的值可以变得非常大，因此常见的处理方法是计算 Gamma 的日志以使其值更小。

计算 LogGamma 非常复杂，你可以使用多种算法进行计算。演示程序使用名为 Lanczos 的算法，其中近似值为 (g=5, n=7)。A&S 引用具有可以计算 LogGamma 的不同算法，但是 Lanczos 的近似值（在 A&S 出版时并不为人所知）给出更准确的结果。

方法 BetaIncCf 的名称代表“由连分数计算的不完整 Beta”。演示程序针对方法 BetaIncCf 使用 A&S 等式 26.5.8。

总结

ANOVA 测试进行三种数学假设：组数据项在数学上是独立的；人数数据组正常分布（如 Gaussian 分布中所示）；且人数数据组具有相等的方差。

有几种方法可以对这些假设进行测试，但是解释其结果非常具有挑战性。问题在于，尽管当数据有一些不正常或具有不等方差时 ANOVA 仍可以运行，真实数据也几乎不可能完全正常且具有完全相等的方差。最重要的是证明 ANOVA 假设非常困难，因此，你在解释结果时应非常保守。

ANOVA 与 t-检验关系密切。T-检验确定两个组的人数平均值是否相等（在仅有示例数据的情况下）。因此，如果你有三个组（如演示中所示），你可以执行三个 t-检验，将组 1 和组 2、组 1 和组 3，以及组 2 和组 3 进行比较，来替代使用 ANOVA。但是，并不推荐使用此方法，因为它将引入名为类型 1 的错误（误报）。

本文中介绍的 ANOVA 类型为单向（或单因素）ANOVA。如果有两个因素，则是双向 ANOVA，另一项不同的技术。

ANOVA 基于通过数据组计算的 F-统计量的值。还有一些使用 F-统计量的其他统计测试。例如，你可以使用 F-统计量来推断两组数据的方差是否相同。

Dr.James McCaffrey 供职于华盛顿地区雷蒙德市沃什湾的 Microsoft Research。他参与过多个 Microsoft 产品的工作，包括 Internet Explorer 和 Bing。Scripto可通过 jammc@microsoft.com 与 McCaffrey 取得联系。

衷心感谢以下 Microsoft 技术专家对本文的审阅： Chris Lee 和 Kirk Olynk

测试运行 - 使用 C# 执行 ANOVA

演示程序

计算 F-统计量

计算 p-值

总结

其他资源