用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 項目，它用非?？茖W的數組語言以及成熟的建立下標方法擴展了 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 值?？梢杂盟鼈儊碓u估線性模型與數據的吻合程度。高的 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 項目，它用非?？茖W的數組語言以及成熟的建立下標方法擴展了 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.anti-gravitydesign.com