2015 年 7 月

第 30 卷，第 7 期

测试运行 - 使用 C# 实现线性回归

作者 James McCaffrey

线性回归问题的目标是根据一个或多个数字预测变量的值预测数字变量的值。例如，您可能想要根据一个人的受教育水平和工作年限以及性别（男性 = 0，女性 = 1）预测他/她的年收入。

通常，要预测的变量称为因变量。通常，预测变量称为自变量。当只有一个预测变量时，这种技术有时称为简单的线性回归。当有两个或更多预测变量时，该技术一般称为多重或多变量线性回归。

了解本文所述观点的一个好方法是查看图 1 中的演示程序。该 C# 演示程序根据教育经历、工作经历和性别预测年收入。该演示以生成 10 个综合数据项开始。受教育水平级别是指介于 12 和 16 之间的值。工作经验是指介于 10 和 30 之间的值。性别是指示变量，其中男性是引用值，编码为 0，女性编码为 1。收入（以千美元为单位）位于最后一列中。在非演示方案中，您可能使用类似于 MatrixLoad 的方法读取文本文件中的数据。

图 1 使用 C# 实现线性回归

在生成综合数据之后，该演示程序使用数据创建所谓的设计矩阵。设计矩阵就是指含有添加了所有 1.0 值的前导列的数据矩阵。对于线性回归，可以使用若干个不同的算法；一些算法可以使用原始数据矩阵，而其他算法可以使用设计矩阵。该演示使用需要设计矩阵的技术。

创建设计矩阵之后，该演示程序会查找四个系数的值（12.0157、1.0180、0.5489、-2.9566）。这些系数有时称为 b 值或 beta 值。第一个值 12.0157 通常称为截距。它是一个常量，与任何预测变量都不相关。第二、第三和第四个系数值（1.0180、0.5489、-2.9566）分别与受教育水平、工作经验和性别相关。

图 1 中输出的最后部分使用这些系数的值来预测受教育水平为 14，工作经验为 12 年以及性别为 0（男性）的假想人员的收入。预测后的收入为 32.86，其计算方式如下：

income = 12.0157 + (1.0180)(14) + (0.5489)(12) + (-2.9566)(0)
                = 12.0157 + 14.2520 + 6.5868 + 0
                = 32.86

换句话说，要使用线性回归进行预测，可将预测值乘以它们对应的系数值，然后求和。这其实非常简单。请注意，前导截距值（本例中是 12.0157）可视为与始终含有值为 1 的预测变量相关联的系数。这一事实从某种程度上解释了设计矩阵中值为 1.0 的列。

线性回归问题的本质是使用原始数据或等效的设计矩阵来计算系数的值。这并不简单。该演示使用称为关闭窗体矩阵反转的方法，该方法也称为普通最小平方法。查找这些系数的值的备用方法包括以迭代方式重新衡量最小平方法、最大似然估计、脊回归、梯度下降和其他内容。

在图 1 中，在进行预测之前，该演示程序会计算称为 R 平方值的指标，其也被称为决定系数。R 平方是介于 0 和 1 之间的值，描述了预测模型与原始数据的匹配程度。这有时表达为“由模型解释的不均匀率”。 从松散的角度来解释的话，R 平方越接近 1，预测模型越好。对于实际数据而言，演示值为 0.7207 或 72% 被视为非常高（好）。

本文假设您至少具有中级 C# 编程技能，但并不假定您完全了解线性回归。虽然该演示程序过长而无法全部展示，但是本文随附了完整的源代码可进行下载。

了解线性回归

通常，使用图表能最好解释线性回归。请看图 2 中的曲线图。图形中的数据表示只根据单个变量（工作经验年限）预测年收入。每个红点对应一个数据点。例如，最左侧的数据项含有工作 = 10 和收入 = 32.06。线性回归会找到两个系数：一个截距以及一个工作变量。结果显示，这些系数的值为 27.00 和 0.43。

图 2 含有一个自变量的线性回归

系数值决定了某行的方程，如图 2 中的蓝色部分所示。该行（系数）最小化实际数据点 (yi) 和预测数据点 (fi) 之间的平方差的总和。图 2 中用虚线显示的 10 个偏差中的两个。显示额第一个偏差是 yi - fi = 28.6 - 32.6 = -4.0。请注意，偏差可以是正数或负数。如果该偏差没有被平方，负数和正数会相互抵消。

图 2 中的图形显示只含一个自变量的简单线性回归的工作方式。多变量线性回归使用多个自变量推广相同的理念，即查找使平方差的总和最小化的系数。

通过直观的表达，线性回归在一组数据点中查找最佳行。此最佳行可用于预测。例如，在图 2 中，如果假想人员的工作经验是 25 年，则蓝色行表示她的预测收入约为 38。

求解最小平方方程

如果一个线性回归问题含有 n 个预测变量，则必须找到 n+1 个系数值，其中每个系数值对应一个预测值与截距值。该演示程序使用最基本的方法来查找系数值。通常，通过图 3 中显示的令人望而生畏的方程来提供系数的值。该方程并非像它首次显示时那样复杂。

图 3 使用指标的线性回归系数解决方案

希腊字母 β 类似于脚本 B，并且表示系数值。请注意，方程中所有字母都以粗体显示，这在数学上表示它们代表多值对象（指标或数组/向量），而非简单的标量值（纯数字）。大写字母 X 表示设计矩阵。指数为 T 的大写字母 X 表示设计矩阵的转置。* 符号表示矩阵乘法。指数为 -1 表示矩阵反转。大写字母 Y 是指因变量值的列向量（含有一个列的矩阵）。因此，求解系数的值实际上表示了解矩阵操作。

图 4 中的图表说明矩阵转置、矩阵乘法和矩阵反转。矩阵的转置就是交换行和列。例如，假设您有一个 2x3 的矩阵，即该矩阵含有 2 行和 3 列。矩阵的转置之后将变成 3x2，其中原始矩阵的行变成转置后矩阵的列。

图 4 用于查找线性回归系数的三个矩阵操作

如果您以前从未遇到矩阵乘法，则可能觉得它有点奇怪。如果您将 (n x m) 大小的矩阵乘以 (m x p) 大小的矩阵，则结果是 (n x p) 大小的矩阵。例如，3x4 * 4x2 矩阵含有 3x2 的大小。有关矩阵乘法的详细讨论不在本文的范围内，但是您看过几个示例之后，可以轻松地了解该过程并且易于在代码中实施。

对于线性回归系数值而言，第三个需要求解的矩阵操作是矩阵反转。遗憾的是，矩阵反转不易于理解且难于实施。对于本文，知道只有当矩阵含有相同的行和列（矩形矩阵）时才会定义矩阵的反转就足够了。图 4 显示了一个 3x3 矩阵及其反转。

有几种算法可用于查找矩阵的反转。该演示程序使用称为 Doolittle 分解的方法。

总之，含有 n 个预测变量的线性回归问题涉及到查找 n+1 系数的值。这可以通过含有矩阵分解、矩阵乘法和矩阵反转的矩阵来完成。虽然转置和乘法非常容易，但是查找矩阵的反转非常困难。

演示程序结构

为了创建演示程序，我启动了 Visual Studio 并选择了控制台应用程序模板。我将该项目命名为 LinearRegression。该程序对 .NET Framework 的依赖程度并不明显，因此，任何 Visual Studio 版本都可以正常运行。

当模板代码加载到编辑器之后，在“解决方案资源管理器”窗口中，我会右键单击文件 Program.cs，并将其重命名为 LinearRegressionProgram.cs。我允许 Visual Studio 自动重新命名 Program 类。在“编辑器”窗口的顶部，我使用相关语句删除了所有内容，一个引用顶级 System 命名空间除外。

图 5 显示了演示程序的整体结构（为节省空间进行了一些较小的修改）。所有程序控制逻辑都包含在 Main 方法中。该演示程序使用静态方法的方法，而不是 OOP 方法。

图 5 线性回归演示程序结构

using System;
namespace LinearRegression
{
  class LinearRegressionProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("Begin linear regression demo");
      // Generate synthetic data
      // Create design matrix
      // Solve for LR coefficients
      // Calculate R-squared value
      // Do a prediction
      Console.WriteLine("End linear regression demo");
      Console.ReadLine();
    }
    static double Income(double x1, double x2,
      double x3, double[] coef) { . . }
    static double RSquared(double[][] data,
      double[] coef) { . . }
    static double[][] DummyData(int rows,
      int seed) { . . }
    static double[][] Design(double[][] data) { . . }
    static double[] Solve(double[][] design) { . . }
    static void ShowMatrix(double[][] m, int dec) { . . }
    static void ShowVector(double[] v, int dec) { . . }
    // ----------
    static double[][] MatrixTranspose(double[][] matrix)
      { . . }
    static double[][] MatrixProduct(double[][] matrixA,
      double[][] matrixB) { . . }
    static double[][] MatrixInverse(double[][] matrix)
      { . . }
    // Other matrix routines here
  }
} // ns

Income 方法使用系数值的数组从含有受教育水平、工作经验和性别的值的输入参数中返回预测收入。RSquared 方法从数据和系数中返回模型的 R 平方值。DummyData 方法生成用于该演示的合成数据。

Design 方法接受数据的矩阵，并返回含有值为 1.0 的前导列的增强设计矩阵。Solve 方法接受设计矩阵，并使用矩阵操作查找线性回归系数。

大部分艰难的工作都由执行矩阵操作的一组静态方法来完成。该演示程序以一种可能最简单的方式定义了一个矩阵，作为很多数组中的一个数组。备用方案是创建程序定义的 Matrix 类，但是我认为该方法有些故弄玄虚。有时，普通的数组要优于程序定义的对象。

MatrixTranspose 返回矩阵的转置。MatrixProduct 方法返回乘以两个矩阵的结果。MatrixInverse 方法返回矩阵的反转。该演示含有很多帮助程序方法。特别是，MatrixInverse 方法会调用帮助程序方法 MatrixDuplicate、MatrixDecompose 和 HelperSolve。

Solve 方法

线性回归演示程序的核心在于 Solve 方法。该方法的定义以下面的代码作为开头：

static double[] Solve(double[][] design)
{
  int rows = design.Length;
  int cols = data[0].Length;
  double[][] X = MatrixCreate(rows, cols - 1);
  double[][] Y = MatrixCreate(rows, 1);
...

单个输入参数是一个设计矩阵。您可能要考虑的一个备用方案是传递源数据矩阵，然后让 Solve 调用帮助程序方法 Design 来获取设计矩阵。帮助程序方法 MatrixCreate 允许为含有特定行和列的矩阵留出空间并返回该矩阵。本机 X 矩阵保留自变量（含有前导值 1.0）的值。本机 Y 矩阵只含有一列，并且保留因变量（本演示中的年收入）的值。

接下来，使用设计矩阵中的值来填充矩阵 X 和 Y 中的单元格。

int j;
for (int i = 0; i < rows; ++i)
{
  for (j = 0; j < cols - 1; ++j)
  {
    X[i][j] = design[i][j];
  }
  Y[i][0] = design[i][j]; // Last column
}

请注意，索引变量 j 是在嵌套循环的外部声明的，因此它可以用于填充 Y 矩阵。X 和 Y 指标可用之后，可以根据图 3 中显示的方程查找线性回归系数：

...
  double[][] Xt = MatrixTranspose(X);
  double[][] XtX = MatrixProduct(Xt, X);
  double[][] inv = MatrixInverse(XtX);
  double[][] invXt = MatrixProduct(inv, Xt);
  double[][] mResult = MatrixProduct(invXt, Y);
  double[] result = MatrixToVector(mResult);
  return result;
} // Solve

在该演示中，X 矩阵含有的大小为 10x4，因此它的转置矩阵 Xt 含有的大小为 4x10。将 Xt 和 X 相乘以后得到大小为 4x4 的矩阵，反转后的矩阵 inv 的大小为 4x4。通常，对于含有 n 个独立预测变量的线性回归问题，在使用矩阵反转方法时，您需要查找含有大小为 (n+1) x (n+1) 的矩阵的反转。这表示该反转方法不适合含有大量预测变量的线性回归问题。

在该代码中，将 4x4 的反转矩阵和 4x10 的转置矩阵相乘后得到的 invXt 矩阵含有的大小为 4x10。在该代码中，将 invXt 矩阵和 10x1 的 Y 矩阵相乘得到的 mResult（“矩阵结果”）含有的大小为 4x1。这些值就是您需要的系数。为了方便起见，使用帮助程序方法 MatrixToVector 将含有单列的 Y 矩阵中的值转移到普通的数组。

计算 R 平方

如前文所述，R 平方指标是衡量实际数据点与计算得出的回归线相匹配程度的指标。在数学术语中，R 平方被定义为 R2 = 1 - (SSres / SStot)。SSres 术语通常称为“残差平方和”。 它是实际 Y 值和预测 Y 值之间的平方差的总和，如图 2 中的图形中所述。SStot 术语是指“总平方和”。 它是每个实际 Y 值和所有实际 Y 值的平均值之间的平方差的总和。

此外，线性回归中的 R 平方指标也称为决定系数，并且与另一个统计指标 r 平方（小 r 平方）相关，但有不同之处。解释 R 平方有点困难，并且取决于调查之后的特定问题领域。对于自然科学和社会科学，其中的数据通常都是混乱且不完整，但是 0.6 或大于 0.6 的 R 平方值通常被视为非常好。

该回归模型解释的方差有一个备用的衡量方法，该方法称为调整后的 R 平方。此指标将预测变量的总量和数据项的总量考虑在内。在大多数情况下，使用普通的 R 平方足以了解线性回归模型的预测能力。

总结

如果您在 Internet 上搜索有关如何使用编程语言执行线性回归的示例，您不会找到很多参考。我认为缺少此类相关信息的主要原因有两个。第一，使用矩阵操作求解线性回归系数非常困难，主要就是因为矩阵反转操作。从某种程度上来说，我认为该演示程序的 MatrixInverse 方法是我撰写的最复杂的代码例程。第二，有很多现有的独立工具可以执行线性回归，尤其是 Excel 电子表格程序及其数据分析外接程序。很少需要直接在软件系统中嵌入线性回归解决方案代码。

对于线性回归的研究已持续数十年，而且有很多方法可扩展这方面技术。例如，您可以引入所谓的交互作用效果，它将两个或更多的预测变量结合起来。这些扩展有时称为一般线性模型，以便于将它们从线性回归的基本形式中区分开来。

我认为线性回归是经典统计学的“入门”方法。虽然经典统计学和机器学习之间的区别没有明确的、统一认识，但是我倾向于将经典统计学视为由数学家在 20 世纪早期开始研究的内容。在我看来，诸如神经网络分类等机器学习技术首次出现在 20 世纪 50 年代，属于比较新的内容。经典统计学线性回归与名为逻辑回归的机器学习技术非常相关，而且后者已经是几个测试运行专栏的主题。

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

衷心感谢以下 Microsoft Research 技术专家对本文的审阅：Charles Parker