解决输出和概率函数缺陷的数据研究工具
本系列文章的第 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. 实现概率函数
// 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 $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. 演示输出类的通用逻辑
// 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 的函数的详细内容
// 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);
?>
' vspace='15'>
?>
数据研究脚本
该数据研究工具由单个脚本(explore.php)构成,该脚本调用 SimpleLinearRegressionHTML 类和 JpGraph 库的方法。
该脚本使用了简单的处理逻辑。该脚本的第一部分对所提交的表单数据执行基本验证。如果这些表单数据通过验证,则执行该脚本的第二部分。
该脚本的第二部分所包含的代码用于分析数据,并以 HTML 和图形格式显示汇总结果。清单 4 中显示了 explore.php 脚本的基本结构:
清单 4. explore.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 "