用 PHP 实现的简单线性回归：(一)

数据库在 PHP 中的重要性


PHP 领域中缺少了一个功能强大的工具：基于语言的数学库。在这个由两部分组成的系列文章中，Paul Meagher 希望通过提供一个如何开发分析模型库的示例来启发 PHP 开发人员去开发和实现基于 PHP 的数学库。在第 1 部分中，他演示了如何使用 PHP 作为实现语言来开发和实现简单线性回归（Simple Linear Regression）算法包的核心部分。在第 2 部分中，作者在该包中添加了一些功能：针对中小规模数据集的有用的数据分析工具。

简介
与其它开放源码语言（比如 Perl 和 Python）相比，PHP 社区缺少强有力的工作来开发数学库。

造成这种状况的一个原因可能是由于已经存在大量成熟的数学工具，这可能阻碍了社区自行开发 PHP 工具的工作。例如，我曾研究过一个功能强大的工具 S System，它拥有一组令人印象深刻的统计库，专门被设计成用来分析数据集，并且在 1998 年由于其语言设计而获得了 ACM 奖。如果 S 或者其开放源码同类 R 仅仅是一个 exec_shell 调用，那么为何还要麻烦用 PHP 实现相同的统计计算功能呢？有关 S System、它的 ACM 奖或 R 的更多信息，请参阅参考资料。

难道这不是在浪费开发人员的精力吗？如果开发 PHP 数学库的动机是出自节省开发人员的精力以及使用最好的工具来完成工作，那么 PHP 现在的课题是很有意义的。

另一方面，出于教学动机可能会鼓励对 PHP 数学库的开发。对于大约 10% 的人来说，数学是个值得探索的有趣课题。对于那些同时还熟练应用 PHP 的人来说，PHP 数学库的开发可以增强数学学习过程，换句话说，不要只阅读有关 T 测试的章节，还要实现一个能计算相应的中间值并用标准格式显示它们的类。

通过指导和训练，我希望证明开发 PHP 数学库并不是一项很难的任务，它可能代表一项有趣的技术和学习难题。在本文中，我将提供一个 PHP 数学库示例，名为 SimpleLinearRegression，它演示了一个可以用来开发 PHP 数学库的通用方法。让我们从讨论一些通用的原则开始，这些原则指导我开发这个 SimpleLinearRegression 类。

指导原则
我使用了六个通用原则来指导 SimpleLinearRegression 类的开发。

1.每个分析模型建立一个类。
2.使用逆向链接来开发类。
3.预计有大量的 getter。
4.存储中间结果。
5.为详细的 API 制定首选项。
6.尽善尽美并非目标。
7.让我们更详细地逐条研究这些指导方针。

每个分析模型建立一个类
每种主要的分析测试或过程应当有一个名称与测试或过程名相同的 PHP 类，这个类包含了输入函数、计算中间值和汇总值的函数和输出函数（将中间值和汇总值用文本或图形格式全部显示在屏幕上）。

使用逆向链接来开发类
在数学编程中，编码的目标通常是分析过程（比如 MultipleRegression、TimeSeries 或 ChiSquared）所希望生成的标准输出值。从解决问题的角度出发，这意味着您可以使用逆向链接来开发数学类的方法。

例如，汇总输出屏幕显示了一个或多个汇总统计结果。这些汇总统计结果依赖于中间统计结果的计算，这些中间统计结果又可能会涉及到更深一层的中间统计结果，以此类推。这个基于逆向链接的开发方法导出了下一个原则。

预计有大量的 getter
数学类的大部分类开发工作都涉及到计算中间值和汇总值。实际上，这意味着，如果您的类包含许多计算中间值和汇总值的 getter 方法，您不应当感到惊讶。

存储中间结果
将中间计算结果存储在结果对象内，这样您就可以将中间结果用作后续计算的输入。在 S 语言设计中实施了这一原则。在当前环境下，通过选择实例变量来表示计算得到的中间值和汇总结果，从而实施了该原则。

为详细的 API 制定首选项
当为 SimpleLinearRegression 类中的成员函数和实例变量制定命名方案时，我发现：如果我使用较长的名称（类似于 getSumSquaredError 这样的名称，而不是 getYY2）来描述成员函数和实例变量，那么就更容易了解函数的操作内容和变量所代表的意义。

我没有完全放弃简写名称；但是，当我用简写形式的名称时，我得设法提供注释以完整阐述该名称的含义。我的看法是：高度简写的命名方案在数学编程中很常见，但它们使得理解和证明某个数学例程是否按部就班更为困难，而原本不必造成此种困难。

尽善尽美并非目标
这个编码练习的目标不是一定要为 PHP 开发高度优化和严格的数学引擎。在早期阶段，应当强调学习实现意义重大的分析测试，以及解决这方面的难题。


实例变量
当对统计测试或过程进行建模时，您需要指出声明哪些实例变量。

实例变量的选择可以通过说明由分析过程生成的中间值和汇总值来确定。每个中间值和汇总值都可以有一个相应的实例变量，将变量的值作为对象属性。

我采用这样的分析来确定为清单 1 中的 SimpleLinearRegression 类声明哪些变量。可以对 MultipleRegression、ANOVA 或 TimeSeries 过程执行类似的分析。

清单 1. SimpleLinearRegression 类的实例变量
<?php

// Copyright 2003, Paul Meagher
// Distributed under GPL

class SimpleLinearRegression {

var $n;
var $X = array();
var $Y = array();
var $ConfInt;
var $Alpha;
var $XMean;
var $YMean;
var $SumXX;
var $SumXY;
var $SumYY;
var $Slope;
var $YInt;
var $PredictedY = array();
var $Error = array();
var $SquaredError = array();
var $TotalError;
var $SumError;
var $SumSquaredError;
var $ErrorVariance;
var $StdErr;
var $SlopeStdErr;
var $SlopeVal; // T value of Slope
var $YIntStdErr;
var $YIntTVal; // T value for Y Intercept
var $R;
var $RSquared;
var $DF; // Degrees of Freedom
var $SlopeProb; // Probability of Slope Estimate
var $YIntProb; // Probability of Y Intercept Estimate
var $AlphaTVal; // T Value for given alpha setting
var $ConfIntOfSlope;

var $RPath = "/usr/local/bin/R"; // Your path here

var $format = "%01.2f"; // Used for formatting output

}
?>


构造函数
SimpleLinearRegression 类的构造函数方法接受一个 X 和一个 Y 向量，每个向量都有相同数量的值。您还可以为您预计的 Y 值设置一个缺省为 95% 的置信区间（confidence interval）。

构造函数方法从验证数据形式是否适合于处理开始。一旦输入向量通过了“大小相等”和“值大于 1”测试，就执行算法的核心部分。

执行这项任务涉及到通过一系列 getter 方法计算统计过程的中间值和汇总值。将每个方法调用的返回值赋给该类的一个实例变量。用这种方法存储计算结果确保了前后链接的计算中的调用例程可以使用中间值和汇总值。还可以通过调用该类的输出方法来显示这些结果，如清单 2 所描述的那样。

清单 2. 调用类输出方法
<?php

// Copyright 2003, Paul Meagher
// Distributed under GPL

function SimpleLinearRegression($X, $Y, $ConfidenceInterval="95") {

$numX = count($X);
$numY = count($Y);

if ($numX != $numY) {
die("Error: Size of X and Y vectors must be the same.");

}
if ($numX <= 1) {
die("Error: Size of input array must be at least 2.");
}

$this->n = $numX;
$this->X = $X;
$this->Y = $Y;

$this->ConfInt = $ConfidenceInterval;
$this->Alpha = (1 + ($this->ConfInt / 100) ) / 2;

$this->XMean = $this->getMean($this->X);
$this->YMean = $this->getMean($this->Y);
$this->SumXX = $this->getSumXX();
$this->SumYY = $this->getSumYY();
$this->SumXY = $this->getSumXY();
$this->Slope = $this->getSlope();
$this->YInt = $this->getYInt();
$this->PredictedY = $this->getPredictedY();
$this->Error = $this->getError();
$this->SquaredError = $this->getSquaredError();
$this->SumError = $this->getSumError();
$this->TotalError = $this->getTotalError();
$this->SumSquaredError = $this->getSumSquaredError();
$this->ErrorVariance = $this->getErrorVariance();
$this->StdErr = $this->getStdErr();
$this->SlopeStdErr = $this->getSlopeStdErr();
$this->YIntStdErr = $this->getYIntStdErr();
$this->SlopeTVal = $this->getSlopeTVal();
$this->YIntTVal = $this->getYIntTVal();
$this->R = $this->getR();
$this->RSquared = $this->getRSquared();
$this->DF = $this->getDF();
$this->SlopeProb = $this->getStudentProb($this->SlopeTVal, $this->DF);
$this->YIntProb = $this->getStudentProb($this->YIntTVal, $this->DF);
$this->AlphaTVal = $this->getInverseStudentProb($this->Alpha, $this->DF);
$this->ConfIntOfSlope = $this->getConfIntOfSlope();

return true;
}

?>


方法名及其序列是通过结合逆向链接和参考大学本科学生使用的统计学教科书推导得出的，该教科书一步一步地说明了如何计算中间值。我需要计算的中间值的名称带有“get”前缀，从而推导出方法名。

使模型与数据相吻合
SimpleLinearRegression 过程用于产生与数据相吻合的直线，其中直线具有以下标准方程：

y = b + mx

该方程的 PHP 格式看起来类似于清单 3：

清单 3. 使模型与数据相吻合的 PHP 方程
$PredictedY[$i] = $YIntercept + $Slope * $X[$i]


SimpleLinearRegression 类使用最小二乘法准则推导出 Y 轴截距（Y Intercept）和斜率（Slope）参数的估计值。这些估计的参数用来构造线性方程（请参阅清单 3），该方程对 X 和 Y 值之间的关系进行建模。

使用推导出的线性方程，您就可以得到每个 X 值对应的预测 Y 值。如果线性方程与数据非常吻合，那么 Y 的观测值与预测值趋近于一致。

如何确定是否非常吻合
SimpleLinearRegression 类生成了相当多的汇总值。一个重要的汇总值是 T 统计值，它可以用来衡量一个线性方程与数据的吻合程度。如果非常吻合，那么 T 统计值往往很大。如果 T 统计值很小，那么应当用一个模型替换该线性方程，该模型假设 Y 值的均值是最佳预测值（也就是说，一组值的均值通常是下一个观测值有用的预测值，使之成为缺省模型）。

要测试 T 统计值是否大得足以不把 Y 值的均值作为最佳预测值，您需要计算获取 T 统计值的随机概率。如果获取 T 统计值的概率很低，那么您可以否定均值是最佳预测值这个无效假设，与此相对应，也就确信简单线性模型与数据非常吻合。

那么，如何计算 T 统计值的概率呢？

计算 T 统计值概率
由于 PHP 缺少计算 T 统计值概率的数学例程，因此我决定将此任务交给统计计算包 R（请参阅参考资料中的 www.r-project.org）来获得必要的值。我还想提醒大家注意该包，因为：

1. R 提供了许多想法，PHP 开发人员可能会在 PHP 数学库中模拟这些想法。
2. 有了 R，可以确定从 PHP 数学库获得的值与那些从成熟的免费可用的开放源码统计包中获得的值是否一致。
清单 4 中的代码演示了交给 R 来处理以获取一个值是多么容易。

清单 4. 交给 R 统计计算包来处理以获取一个值
<?php

// Copyright 2003, Paul Meagher
// Distributed under GPL

class SimpleLinearRegression {

var $RPath = "/usr/local/bin/R"; // Your path here

function getStudentProb($T, $df) {
$Probability = 0.0;
$cmd = "echo @#dt($T, $df)@# | $this->RPath --slave";
$result = shell_exec($cmd);
list($LineNumber, $Probability) = explode(" ", trim($result));
return $Probability;
}

function getInverseStudentProb($alpha, $df) {
$InverseProbability = 0.0;
$cmd = "echo @#qt($alpha, $df)@# | $this->RPath --slave";
$result = shell_exec($cmd);
list($LineNumber, $InverseProbability) = explode(" ", trim($result));
return $InverseProbability;
}

}

?>


请注意，这里已经设置了到 R 可执行文件的路径，并在两个函数中使用了该路径。第一个函数根据学生的 T 分布返回了与 T 统计值相关的概率值，而第二个反函数计算了与给定的 alpha 设置相对应的 T 统计值。getStudentProb 方法用来评估线性模型的吻合程度；getInverseStudentProb 方法返回一个中间值，它用来计算每个预测的 Y 值的置信区间。

由于篇幅有限，我不可能逐个详细说明这个类中的所有函数，因此如果您想搞清楚简单线性回归分析中所涉及的术语和步骤，我鼓励您参考大学本科学生使用的统计学教科书。

燃耗研究
要演示如何使用该类，我可以使用来自公共事业中燃耗（burnout）研究中的数据。Michael Leiter 和 Kimberly Ann Meechan 研究了称为消耗指数（Exhaustion Index)的燃耗度量单位和称之为集中度（Concentration）的独立变量之间的关系。集中度是指人们的社交接触中来自其工作环境的那部分比例。

要研究他们样本中个人的消耗指数值与集中度值之间的关系，请将这些值装入适当命名的数组中，并用这些数组值对该类进行实例化。对类进行实例化后，显示该类所生成的某些汇总值以评估线性模型与数据的吻合程度。

清单 5 显示了装入数据和显示汇总值的脚本：

清单 5. 用于装入数据并显示汇总值的脚本
<?php

// BurnoutStudy.php

// Copyright 2003, Paul Meagher
// Distributed under GPL

include "SimpleLinearRegression.php";

// Load data from burnout study

$Concentration = array(20,60,38,88,79,87,
68,12,35,70,80,92,
77,86,83,79,75,81,
75,77,77,77,17,85,96);

$ExhaustionIndex = array(100,525,300,980,310,900,
410,296,120,501,920,810,
506,493,892,527,600,855,
709,791,718,684,141,400,970);

$slr = new SimpleLinearRegression($Concentration, $ExhaustionIndex);

$YInt = sprintf($slr->format, $slr->YInt);
$Slope = sprintf($slr->format, $slr->Slope);
$SlopeTVal = sprintf($slr->format, $slr->SlopeTVal);
$SlopeProb = sprintf("%01.6f", $slr->SlopeProb);

?>

<table border=@#1@# cellpadding=@#5@#>
<tr>
<th align=@#right@#>Equation:</th>
<td></td>
</tr>
<tr>
<th align=@#right@#>T:</th>
<td></td>
</tr>
<tr>
<th align=@#right@#>Prob > T:</th>
<td><td>
</tr>
</table>


通过 Web 浏览器运行该脚本，产生以下输出：

Equation: Exhaustion = -29.50 + (8.87 * Concentration)
T: 6.03
Prob > T: 0.000005


这张表的最后一行指出获取这样大 T 值的随机概率非常低。可以得出这样的结论：与仅仅使用消耗值的均值相比，简单线性模型的预测能力更好。

知道了某个人的工作场所联系的集中度，就可以用来预测他们可能正在消耗的燃耗程度。这个方程告诉我们：集中度值每增加 1 个单位，社会服务领域中一个人的消耗值就会增加 8 个单位。这进一步证明了：要减少潜在的燃耗，社会服务领域中的个人应当考虑在其工作场所之外结交朋友。

这只是粗略地描述了这些结果可能表示的含义。为全面研究这个数据集的含义，您可能想更详细地研究这个数据以确信这是正确的解释。在下一篇文章中我将讨论应当执行其它哪些分析。

您学到了什么？
其一，要开发意义重大的基于 PHP 的数学包，您不必是一名火箭科学家。坚持标准的面向对象技术，以及明确地采用逆向链接问题解决方法，就可以相对方便地使用 PHP 实现某些较为基本的统计过程。

从教学的观点出发，我认为：如果只是因为要求您在较高和较低的抽象层次思考统计测试或例程，那么这个练习是非常有用的。换句话说，补充您的统计测试或过程学习的一个好办法就是将这个过程作为算法实现。

要实现统计测试通常需要超出所给定的信息范围并创造性地解决和发现问题。对于发现对某个学科认识的不足而言，它也是一个好办法。

不利的一面，您发现 PHP 对于取样分布缺乏内在手段，而这是实现大多数统计测试所必需的。您需要交给 R 来处理以获取这些值，但是我担心您会没时间或没兴趣安装 R。某些常见概率函数的本机 PHP 实现可以解决这个问题。

另一个问题是：该类生成许多中间值和汇总值，但是汇总输出实际上没有利用这一点。我提供了一些难处理的输出，但是这既不够充分也没进行很好的组织，以致您无法充分地解释分析结果。实际上，我完全不知道如何可以将输出方法集成到该类中。这需要得到解决。

最后，要弄明白数据，不仅仅是察看汇总值就可以了。您还需要明白各个数据点是如何分布的。最好的办法之一是将您的数据绘制成图表。再次声明，我对这方面不太了解，但是如果要用这个类来分析实际数据的话就需要解决这个问题。

在本系列文章的下一篇文章中，我将使用本机 PHP 代码实现一些概率函数，用几个输出方法扩展 SimpleLinearRegression 类，并生成一个报告：用表和图形格式表示中间值和汇总值，这样更容易从数据中得出结论。且待下回分解！


参考资料

1.请参考由 James T. McClave 和 Terry Sincich 编著的广受欢迎的大学教科书 Statistics，第 9 版（Prentice-Hall，在线），本文中所使用的算法步骤和“燃耗研究”示例参考了该书。
2.请查阅 PEAR 资源库，它目前包含了少量低级别的 PHP 数学类。最终，应该会很高兴地看到 PEAR 包含实现标准的较高级别的数值方法（比如 SimpleLinearRegression、MultipleRegression、TimeSeries、ANOVA、FactorAnalysis、FourierAnalysis 及其它）的包。
3.查看作者的 SimpleLinearRegression 类的所有源代码。
4.了解一下Numerical Python 项目，它用非常科学的数组语言以及成熟的建立下标方法扩展了 Python。有了该扩展，数学操作就非常接近人们期望从编译语言所获得的功能。
5.研究可用于 Perl 的许多数学参考资料，包括 CPAN 数学模块的索引和 CPAN 中算法部分的模块，以及 Perl 数据语言（Perl Data Language），它旨在为 Perl 提供压缩存储以及快速操作大型 N 维数据数组的能力。
6.有关 John Chambers 的 S 编程语言的更多信息，请查阅关于他的出版物以及他在贝尔实验室的各项研究项目的链接。还可以了解在 1998 年因语言设计而获得的 ACM 奖。
7.R 是用于统计计算和图形的语言和环境，类似于获奖的 S System，R 提供了诸如线性和非线性建模、统计测试、时间序列分析、分类、群集之类的统计和图形技术。请在 R Project 主页上了解 R。
8.如果您刚接触 PHP，那么请阅读 Amol Hatwar 的 developerWorks 系列文章：“用 PHP 开发健壮的代码:”“第 1 部分: 高屋建瓴的介绍 ”（2002 年 8 月）、“第 2 部分: 有效地使用变量”（2002 年 9 月）和“第 3 部分: 编写可重用函数”（2002 年 11 月）。

原文转自：http://www.ltesting.net