用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 月）。
解决输出和概率函数缺陷的数据研究工具

本系列文章的第 1 部分结尾处提到了简单线性回归（Simple Linear Regression）类中缺少的三个元素。在本文中，作者 Paul Meagher 用基于 PHP 的概率函数弥补了这些缺陷，演示了如何将输出方法集成到 SimpleLinearRegression 类中并创建了图形输出。他通过构建数据研究工具解决了这些问题，该工具旨在深层次地研究中小规模的数据集所包含的信息。（在第 1 部分中，作者演示了如何用 PHP 作为实现语言来开发和实现简单线性回归算法包的核心部分。）


在这个由两部分组成的系列文章的第 1 部分（“用 PHP 实现的简单线性回归”）中，我说明了数学库对 PHP 有用的原因。我还演示了如何用 PHP 作为实现语言来开发和实现简单线性回归算法的核心部分。

本文的目标是向您展示如何使用第 1 部分中讨论的 SimpleLinearRegression 类来构建一个重要的数据研究工具。

简要回顾：概念
简单线性回归建模背后的基本目标是从成对的 X 值和 Y 值（即 X 和 Y 测量值）组成的二维平面中找到最吻合的直线。一旦用最小方差法找到这条直线，就可以执行各种统计测试，以确定这条直线与观测到的 Y 值的偏离量吻合程度。

线性方程（y = mx + b）有两个参数必须根据所提供的 X 和 Y 数据估算出来，它们是斜率（m）和 y 轴截距（b）。一旦估算出这两个参数，就可以将观测值输入线性方程，并观察方程所生成的 Y 预测值。

要使用最小方差法估算出 m 和 b 参数，就要找到 m 和 b 的估计值，使它们对于所有的 X 值得到的 Y 值的观测值和预测值最小。观测值和预测值之差称为误差（yi - (mxi + b)），并且，如果对每个误差值都求平方，然后求这些残差的和，其结果是一个被称为预测平方差的数。使用最小方差法来确定最吻合的直线涉及寻找使预测方差最小的 m 和 b 的估计值。

可以用两种基本方法来找到满足最小方差法的估计值 m 和 b。第一种方法，可以使用数值搜索过程设定不同的 m 和 b 值并对它们求值，最终决定产生最小方差的估计值。第二种方法是使用微积分找到用于估算 m 和 b 的方程。我不打算深入讨论推导出这些方程所涉及的微积分，但我确实在 SimpleLinearRegression 类中使用了这些分析方程，以找到 m 和 b 的最小平方估计值（请参阅 SimpleLinearRegression 类中的 getSlope() 和 getYIntercept 方法）。

即使拥有了可以用来找到 m 和 b 的最小平方估计值的方程，也并不意味着只要将这些参数代入线性方程，其结果就是一条与数据良好吻合的直线。这个简单线性回归过程中的下一步是确定其余的预测方差是否可以接受。

可以使用统计决策过程来否决“直线与数据吻合”这个备择假设。这个过程基于对 T 统计值的计算，使用概率函数求得随机大的观测值的概率。正如第 1 部分所提到的，SimpleLinearRegression 类生成了为数众多的汇总值，其中一个重要的汇总值是 T 统计值，它可以用来衡量线性方程与数据的吻合程度。如果吻合良好，则 T 统计值往往是一个较大的值；如果 T 值很小，就应该用一个缺省模型代替您的线性方程，该模型假定 Y 值的平均值是最佳预测值（因为一组值的平均值通常可以是下一个观测值的有用的预测值）。

要测试 T 统计值是否大到可以不用 Y 值的平均值作为最佳预测值，需要计算随机获得 T 统计值的概率。如果概率很低，那就可以不采用平均值是最佳预测值这一无效假设，并且相应地可以确信简单线性模型是与数据良好吻合的。（有关计算 T 统计值概率的更多信息，请参阅第 1 部分。）

回过头讨论统计决策过程。它告诉您何时不采用无效假设，却没有告诉您是否接受备择假设。在研究环境中，需要通过理论参数和统计参数来建立线性模型备择假设。

您将构建的数据研究工具实现了用于线性模型（T 测试）的统计决策过程，并提供了可以用来构造理论和统计参数的汇总数据，这些参数是建立线性模型所需要的。数据研究工具可以归类为决策支持工具，供知识工作者在中小规模的数据集中研究模式。

从学习的角度来看，简单线性回归建模值得研究，因为它是理解更高级形式的统计建模的必由之路。例如，简单线性回归中的许多核心概念为理解多次回归（Multiple Regression）、要素分析（Factor Analysis）和时间序列（Time Series）等建立了良好的基础。

简单线性回归还是一种多用途的建模技术。通过转换原始数据（通常用对数或幂转换），可以用它来为曲线数据建模。这些转换可以使数据线性化，这样就可以使用简单线性回归来为数据建模。所生成的线性模型将被表示为与被转换值相关的线性公式。

概率函数

在前一篇文章中，我通过交由 R 来求得概率值，从而避开了用 PHP 实现概率函数的问题。我对这个解决方案并非完全满意，因此我开始研究这个问题：开发基于 PHP 的概率函数需要些什么。

我开始上网查找信息和代码。一个两者兼有的来源是书籍 ［url=http://www.library.cornell.edu/nr/bookcpdf.html］Numerical Recipes in C [/url] 中的概率函数。我用 PHP 重新实现了一些概率函数代码（gammln.c 和 betai.c 函数），但我对结果还是不满意。与其它一些实现相比，其代码似乎多了些。此外，我还需要反概率函数。

幸运的是，我偶然发现了 John Pezzullo 的 Interactive Statistical Calculation。John 关于概率分布函数的网站上有我需要的所有函数，为便于学习，这些函数已用 JavaScript 实现。

我将 Student T 和 Fisher F 函数移植到了 PHP。我对 API 作了一点改动，以便符合 Java 命名风格，并将所有函数嵌入到名为 Distribution 的类中。该实现的一个很棒的功能是 doCommonMath 方法，这个库中的所有函数都重用了它。我没有花费力气去实现的其它测试（正态测试和卡方测试）也都使用 doCommonMath 方法。

这次移植的另一个方面也值得注意。通过使用 JavaScript，用户可以将动态确定的值赋给实例变量，譬如：

var PiD2 = pi() / 2


在 PHP 中不能这样做。只能把简单的常量值赋给实例变量。希望在 PHP5 中会解决这个缺陷。

请注意清单 1 中的代码并未定义实例变量 — 这是因为在 JavaScript 版本中，它们是动态赋予的值。

清单 1. 实现概率函数
<?php

// Distribution.php

// Copyright John Pezullo
// Released under same terms as PHP.
// PHP Port and OO'fying by Paul Meagher

class Distribution {

function doCommonMath($q, $i, $j, $b) {

$zz = 1;
$z = $zz;
$k = $i;


while($k <= $j) {
$zz = $zz * $q * $k / ($k - $b);
$z = $z + $zz;
$k = $k + 2;
}
return $z;
}

function getStudentT($t, $df) {

$t = abs($t);
$w = $t / sqrt($df);
$th = atan($w);

if ($df == 1) {
return 1 - $th / (pi() / 2);
}

$sth = sin($th);
$cth = cos($th);

if( ($df % 2) ==1 ) {
return
1 - ($th + $sth * $cth * $this->doCommonMath($cth * $cth, 2, $df - 3, -1))
/ (pi()/2);
} else {
return 1 - $sth * $this->doCommonMath($cth * $cth, 1, $df - 3, -1);
}

}

function getInverseStudentT($p, $df) {

$v = 0.5;
$dv = 0.5;
$t = 0;

while($dv > 1e-6) {
$t = (1 / $v) - 1;
$dv = $dv / 2;
if ( $this->getStudentT($t, $df) > $p) {
$v = $v - $dv;
} else {
$v = $v + $dv;
}
}
return $t;
}


function getFisherF($f, $n1, $n2) {
// implemented but not shown
}

function getInverseFisherF($p, $n1, $n2) {
// implemented but not shown
}

}
?>




输出方法
既然您已经用 PHP 实现了概率函数，那么开发基于 PHP 的数据研究工具剩下的唯一难题就是设计用于显示分析结果的方法。

简单的解决方案是根据需要将所有实例变量的值都显示到屏幕上。在第一篇文章中，当显示燃耗研究（Burnout Study）的线性方程、T 值和 T 概率时，我就是这么做的。能根据特定目的而访问特定值是很有帮助的，SimpleLinearRegression 支持此类用法。

然而，另一种用于输出结果的方法是将输出的各部分系统化地进行分组。如果研究用于回归分析的主要统计软件包的输出，就会发现它们往往是用同样的方式对输出进行分组的。它们往往有摘要表（Summary Table）、偏离值分析（Analysis Of Variance）表、参数估计值（Parameter Estimate）表和 R 值（R Value）。类似地，我创建了一些输出方法，名称如下：

showSummaryTable()
showAnalysisOfVariance()
showParameterEstimates()
showRValues()

我还有一个用于显示线性预测公式的方法（getFormula()）。许多统计软件包不输出公式，而是希望用户根据上述方法的输出构造公式。部分是由于您最后用来对数据建模的公式的最终形式可能由于下列原因而与缺省公式不同：

1.Y 轴截距没有有意义的解释
2.或者输入值可能是经过转换的，而您可能需要取消对它们的转换以获取最终的解释。

所有这些方法都假定输出媒介是网页。考虑到您有可能希望用非网页的其它媒介输出这些汇总值，所以我决定将这些输出方法包装在一个继承了 SimpleLinearRegression 类的类中。清单 2 中的代码旨在演示输出类的通用逻辑。为了使通用逻辑更突出，所以除去了实现各种 show 方法的代码。

清单 2. 演示输出类的通用逻辑
<?php

// HTML.php

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

include_once "slr/SimpleLinearRegression.php";

class SimpleLinearRegressionHTML extends SimpleLinearRegression {

function SimpleLinearRegressionHTML($X, $Y, $conf_int) {
SimpleLinearRegression::SimpleLinearRegression($X, $Y, $conf_int);
}

function showTableSummary($x_name, $y_name) { }

function showAnalysisOfVariance() { }

function showParameterEstimates() { }

function showFormula($x_name, $y_name) { }

function showRValues() {}
}

?>




这个类的构造函数只是 SimpleLinearRegression 类构造函数的包装器。这意味着如果您想显示 SimpleLinearRegression 分析的 HTML 输出，则应该实例化 SimpleLinearRegressionHTML 类，而不是直接实例化 SimpleLinearRegression 类。其优点是不会有许多未使用的方法充斥 SimpleLinearRegression 类，并且可以更自由地定义用于其它输出媒介的类（也许会对不同媒介类型实现同一 API）。

图形输出
迄今为止，您已经实现的输出方法都以 HTML 格式显示汇总值。它也适合于用 GIF、JPEG 或 PNG 格式显示这些数据的分布图（scatter plot）或线图（line plot）。

与其亲自编写生成线图和分布图的代码，我认为最好使用名为 JpGraph 的基于 PHP 的图形库。JpGraph 正由 Johan Persson 积极开发，其项目网站这样描述它：


无论是对于只有最少代码的“以快捷但不恰当方式获得的”图形，还是对于需要非常细粒度控制的复杂专业图形，JpGraph 都可以使它们的绘制变得简单。JpGraph 同样适用于科学和商业类型的图形。

JpGraph 分发版中包含大量可以根据特定需求进行定制的示例脚本。将 JpGraph 用于数据研究工具非常简单，只需找到功能与我的需求类似的示例脚本，然后对该脚本进行改写以满足我的特定需求即可。

清单 3 中的脚本是从样本数据研究工具（explore.php）中抽取的，它演示了如何调用该库以及如何将来自于 SimpleLinearRegression 分析的数据填入 Line 和 Scatter 类。这段代码中的注释是 Johan Persson 编写的（JPGraph 代码库的文档化工作做得很好）。

清单 3. 来自于样本数据研究工具 explore.php 的函数的详细内容<?php

// Snippet extracted from explore.php script

include ("jpgraph/jpgraph.php");
include ("jpgraph/jpgraph_scatter.php");
include ("jpgraph/jpgraph_line.php");

// Create the graph
$graph = new Graph(300,200,'auto');
$graph->SetScale("linlin");

// Setup title
$graph->title->Set("$title");
$graph->img->SetMargin(50,20,20,40);
$graph->xaxis->SetTitle("$x_name","center");
$graph->yaxis->SetTitleMargin(30);
$graph->yaxis->title->Set("$y_name");

$graph->title->SetFont(FF_FONT1,FS_BOLD);

// make sure that the X-axis is always at the
// bottom at the plot and not just at Y=0 which is
// the default position
$graph->xaxis->SetPos('min');

// Create the scatter plot with some nice colors
$sp1 = new ScatterPlot($slr->Y, $slr->X);
$sp1->mark->SetType(MARK_FILLEDCIRCLE);
$sp1->mark->SetFillColor("red");
$sp1->SetColor("blue");
$sp1->SetWeight(3);
$sp1->mark->SetWidth(4);

// Create the regression line
$lplot = new LinePlot($slr->PredictedY, $slr->X);
$lplot->SetWeight(2);
$lplot->SetColor('navy');

// Add the pltos to the line
$graph->Add($sp1);
$graph->Add($lplot);

// ... and stroke
$graph_name = "temp/test.png";
$graph->Stroke($graph_name);
?>
<img src=><?php echo $graph_name ?>' vspace='15'>

?>


数据研究脚本
该数据研究工具由单个脚本（explore.php）构成，该脚本调用 SimpleLinearRegressionHTML 类和 JpGraph 库的方法。

该脚本使用了简单的处理逻辑。该脚本的第一部分对所提交的表单数据执行基本验证。如果这些表单数据通过验证，则执行该脚本的第二部分。

该脚本的第二部分所包含的代码用于分析数据，并以 HTML 和图形格式显示汇总结果。清单 4 中显示了 explore.php 脚本的基本结构：

清单 4. explore.php 的结构
<?php

// explore.php

if (!empty($x_values)) {
$X = explode(",", $x_values);
$numX = count($X);
}

if (!empty($y_values)) {
$Y = explode(",", $y_values);
$numY = count($Y);
}

// display entry data entry form if variables not set

if ( (empty($title)) OR (empty($x_name)) OR (empty($x_values)) OR
(empty($y_name)) OR (empty($conf_int)) OR (empty($y_values)) OR
($numX != $numY) ) {

// Omitted code for displaying entry form

} else {

include_once "slr/SimpleLinearRegressionHTML.php";
$slr = new SimpleLinearRegressionHTML($X, $Y, $conf_int);

echo "<h2>$title</h2>";

$slr->showTableSummary($x_name, $y_name);
echo "<br><br>";

$slr->showAnalysisOfVariance();
echo "<br><br>";

$slr->showParameterEstimates($x_name, $y_name);
echo "<br>";

$slr->showFormula($x_name, $y_name);
echo "<br><br>";

$slr->showRValues($x_name, $y_name);
echo "<br>";

include ("jpgraph/jpgraph.php");
include ("jpgraph/jpgraph_scatter.php");
include ("jpgraph/jpgraph_line.php");

// The code for displaying the graphics is inline in the
// explore.php script. The code for these two line plots
// finishes off the script:

// Omitted code for displaying scatter plus line plot
// Omitted code for displaying residuals plot

}

?>


火灾损失研究
为了演示如何使用数据研究工具，我将使用来自假想的火灾损失研究的数据。这个研究将主要住宅区火灾损失的金额与它们到最近消防站的距离关联起来。例如，出于确定保险费的目的，保险公司会对这种关系的研究感兴趣。

该研究的数据如图 1 中的输入屏幕所示。

图 1. 显示研究数据的输入屏幕



数据被提交之后，会对它进行分析，并显示这些分析的结果。第一个显示的结果集是 Table Summary，如图 2 所示。

图 2. Table Summary 是所显示的第一个结果集



Table Summary 以表格形式显示了输入数据和其它列，这些列指出了对应于观测值 X 的预测值 Y、Y 值的预测值和观测值之间的差以及预测 Y 值置信区间的下限和上限。

图 3 显示了 Table Summary 之后的三个高级别数据汇总表。

图 3. 显示了 Table Summary 之后的三个高级别数据汇总表



Analysis of Variance 表显示了如何将 Y 值的偏离值归为两个主要的偏离值来源，由模型解释的方差（请看 Model 行）和模型不能解释的方差（请看 Error 行）。较大的 F 值意味着该线性模型捕获了 Y 测量值中的大多数偏离值。这个表在多次回归环境中更有用，在那里每个独立变量都在表中占有一行。

Parameter Estimates 表显示了估算的 Y 轴截距（Intercept）和斜率（Slope）。每行都包括一个 T 值以及观测到极限 T 值的概率（请看 Prob > T 列）。斜率的 Prob > T 可用于否决线性模型。

如果 T 值的概率大于 0.05（或者是类似的小概率），那么您可以否决该无效假设，因为随机观测到极限值的可能性很小。否则您就必须使用该无效假设。

在火灾损失研究中，随机获得大小为 12.57 的 T 值的概率小于 0.00000。这意味着对于与该研究中观测到的 X 值区间相对应的 Y 值而言，线性模型是有用的预测器（比 Y 值的平均值更好）。

最终报告显示了相关性系数或 R 值。可以用它们来评估线性模型与数据的吻合程度。高的 R 值表明吻合良好。

每个汇总报告对有关线性模型和数据之间关系的各种分析问题提供了答案。请查阅 Hamilton、Neter 或 Pedhauzeur 编写的教科书，以了解更高级的回归分析处理(请参阅参考资料）。

要显示的最终报告元素是数据的分布图和线图，如图 4 所示。

图 4. 最终报告元素 — 分布图和线图



大多数人都熟悉线图（如本系列中的第一幅图）的说明，因此我将不对此进行注释，只想说 JPGraph 库可以产生用于 Web 的高质量科学图表。当您输入分布或直线数据时，它也做得很好。

第二幅图将残差（观测的 Y、预测的 Y）与您预测的 Y 值关联起来。这是研究性数据分析（Exploratory Data Analysis，EDA）的倡导者所使用的图形示例，用以帮助将分析人员对数据中的模式的检测和理解能力提到最高程度。行家可以使用这幅图回答关于下列方面的问题：

可能的非正常值或影响力过度的例子
可能的曲线关系（使用转换？）
非正态残差分布
非常量误差方差或异方差性

可以轻松地扩展这个数据研究工具，以生成更多类型的图形 — 直方图、框图和四分位数图 — 这些都是标准的 EDA 工具。

数学库体系结构
对数学的业余爱好使我在最近几个月中保持着对数学库的浓厚兴趣。此类研究推动我思考如何组织我的代码库以及使其预期在未来能不断增长。

我暂时采用清单 5 中的目录结构：

清单 5. 易于增长的目录结构
phpmath/

burnout_study.php
explore.php
fire_study.php
navbar.php

dist/
Distribution.php
fisher.php
student.php
source.php

jpgraph/
etc...

slr/
SimpleLinearRegression.php
SimpleLinearRegressionHTML.php

temp/




例如，未来有关多次回归的工作，将涉及扩展这个库以包括 matrix 目录，该目录用来容纳执行矩阵操作（这是对于更高级形式的回归分析的需求）的 PHP 代码。我还将创建一个 mr 目录，以容纳实现多次回归分析输入方法、逻辑和输出方法的 PHP 代码。

请注意这个目录结构包含一个 temp 目录。必须设置该目录的许可权，使 explore.php 脚本能够将输出图写到该目录。在尝试安装 phpmath_002.tar.gz 源代码时请牢记这一点。此外，请在 JpGraph 项目网站上阅读安装 JpGraph 的指示信息（请参阅参考资料）。

最后提一点，如果采取以下作法，可以将所有软件类移到 Web 根目录之外的文档根目录：


使某个全局 PHP_MATH 变量有权访问非 Web 根目录位置，并且
确保在所有需要或包括的文件路径前面加上这个已定义的常量作为前缀。

将来，对 PHP_MATH 变量的设置将通过一个用于整个 PHP 数学库的配置文件来完成。

您学到了什么？
在本文中，您了解了如何使用 SimpleLinearRegression 类开发用于中小规模的数据集的数据研究工具。在此过程中，我还开发了一个供 SimpleLinearRegression 类使用的本机概率函数，并用 HTML 输出方法和基于 JpGraph 库的图形生成代码扩展该类。

从学习的角度来看，简单线性回归建模是值得进一步研究的，因为事实证明，它是理解更高级形式的统计建模的必由之路。在深入学习更高级的技术（如多次回归或多变量方差分析）之前，对于简单线性回归的透彻理解将使您受益匪浅。

即使简单线性回归只用一个变量来说明或预测另一个变量的偏离值，在所有的研究变量之间寻找简单线性关系仍然常常是研究性数据分析的第一步。仅因为数据是多元的并不意味着就必须使用多元工具研究它。实际上，在开始时使用简单线性回归这样的基本工具是着手探究数据模式的好方法。

本系列研究了简单线性回归分析的两个应用。在本文中，我研究了“到消防站的距离”和“火灾损失”之间的强线性关系。在第一篇文章中，我研究了“社会集中度”和称为“消耗指数”的测量值之间的线性关系，尽管这种关系相对弱一些，但仍然十分明显。（作为练习，用本文中讨论的数据研究工具重新研究第一个研究案例中较为凌乱的数据可能会很有趣。您可能会注意到 y 轴截距是负数的情况，这意味着“社会集中度”为 0，预测消耗指数为 -29.50。这有意义吗？在对一种现象建模时，您应该问问自己：方程是否应该包含可选的 y 轴截距，如果可以，那么该 y 轴截距在线性方程中会起什么作用。）

对于简单线性回归的进一步研究可能包括对这些主题的研究：


*如果想从您的方程以及可以使用的其它计算公式中略去截距，则何时可以这样做
*何时以及如何使用幂、对数和其它转换来对数据进行线性化，以便用简单线性回归来对该数据建模
*可以用来评估您的建模假设的充分性，并可以更清晰地洞察数据中的模式的其它可视化方法

这些是有待学习简单线性回归的学生研究的一部分更高级的主题。参考资料包含了几个指向高级主题文章的链接，您可以从中参考更多关于回归分析的信息。

标准 PHP 安装提供了开发基于数学的重要应用程序所必须的许多资源。我希望这个系列的文章能启发其他开发人员出于乐趣、技术或学习挑战的目的而用 PHP 来实现数学例程。


相关附件：本文所用源代码下载

参考资料
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 月）。
9.访问 John Pezzullo 的优秀站点，该站点专门提供执行统计计算的网页。基于 PHP 的概率函数是以在 John 的概率函数页面所找到的代码为基础的。
10.到 Digital Library of Mathematical Functions 了解关于 M. Abramowitz 和 I.A. Stegun 编写的书籍 The Handbook of Mathematical Functions（也称为 AMS55）的更多信息。
11.查看 JpGraph 站点，以获取关于 PHP 的主要 OO 图形库的大量信息。
12.阅读美国国家标准与技术研究所（National Institute of Standards，NIST）出版的 The Engineering Handbook of Statistics，该手册上有几章是关于 Exploratory Data Analysis 的，非常不错。
13.如果您对于更详尽地学习关于回归的主题感兴趣的话，请尝试阅读以下有用的参考资料：

L. C. Hamilton（1992年）。Regression with Graphics。加州 Pacific Grove：Brooks/Cole Publishing Company。
J Neter、M.H. Kutner 和 W Wasserman W（1990 年）。Applied Linear Regression Models（第 3 版）。芝加哥 Irwin。
E. J. Pedhazur（1982 年）。Multiple regression in behavioral research。纽约州，纽约市：Holt，Rinehart and Winston。

14.阅读 Cameron Laird 的文章“Open source in the biosciences”。PHP 需要更好的数学工具来参与这个不断成长的市场（developerWorks，2002 年 11 月）。
15.查看 RWeb，它是基于 Web 的 R 接口。

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