PHP實現簡單線性迴歸之數學庫的重要性

簡介

與其它開放源碼語言（比如 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 類的開發。

每個分析模型建立一個類。
使用逆向連結來開發類。
預計有大量的 getter。
儲存中間結果。
為詳細的 API 制定喜好設定。
盡善盡美並非目標。

讓我們更詳細地逐條研究這些指導方針。
　
　　每個分析模型建立一個類

每種主要的分析測試或過程應當有一個名稱與測試或過程名相同的 PHP 類，這個類包含了輸入函數、計算中間值和匯總值的函數和輸出函數（將中間值和匯總值用文本或圖形格式全部顯示在螢幕上）。

使用逆向連結來開發類

在數學編程中，編碼的目標通常是分析過程（比如 MultipleRegression 、 TimeSeries 或 ChiSquared ）所希望產生的標準輸出值。從解決問題的角度出發，這意味著您可以使用逆向連結來開發數學類的方法。

例如，匯總輸出螢幕顯示了一個或多個匯總統計結果。這些匯總統計結果依賴於中間統計結果的計算，這些中間統計結果又可能會涉及到更深一層的中間統計結果，以此類推。這個基於逆向連結的開發方法匯出了下一個原則。

預計有大量的 getter

數學類的大部分類開發工作都涉及到計算中間值和匯總值。實際上，這意味著，如果您的類包含許多計算中間值和匯總值的 getter 方法，您不應當感到驚訝。

儲存中間結果

將中間計算結果儲存在結果對象內，這樣您就可以將中間結果用作後續計算的輸入。在 S 語言設計中實施了這一原則。在當前環境下，通過選擇執行個體變數來表示計算得到的中間值和匯總結果，從而實施了該原則。

為詳細的 API 制定喜好設定

當為 SimpleLinearRegression 類中的成員函數和執行個體變數制定命名方案時，我發現：如果我使用較長的名稱（類似於 getSumSquaredError 這樣的名稱，而不是 getYY2 ）來描述成員函數和執行個體變數，那麼就更容易瞭解函數的操作內容和變數所代表的意義。

我沒有完全放棄簡寫名稱；但是，當我用簡寫形式的名稱時，我得設法提供注釋以完整闡述該名稱的含義。我的看法是：高度簡寫的命名方案在數學編程中很常見，但它們使得理解和證明某個數學常式是否按部就班更為困難，而原本不必造成此種困難。

盡善盡美並非目標

這個編碼練習的目標不是一定要為 PHP 開發高度最佳化和嚴格的數學引擎。在早期階段，應當強調學習實現意義重大的分析測試，以及解決這方面的難題。

執行個體變數

當對統計測試或過程進行建模時，您需要指出聲明哪些執行個體變數。

執行個體變數的選擇可以通過說明由分析過程產生的中間值和匯總值來確定。每個中間值和匯總值都可以有一個相應的執行個體變數，將變數的值作為對象屬性。

我採用這樣的分析來確定為清單 1 中的 SimpleLinearRegression 類聲明哪些變數。可以對 MultipleRegression 、 ANOVA 或 TimeSeries 過程執行類似的分析。

清單 1. SimpleLinearRegression 類的執行個體變數

// 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. 調用類輸出方法

// 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）來獲得必要的值。我還想提醒大家注意該包，因為：

R 提供了許多想法，PHP 開發人員可能會在 PHP 數學庫中類比這些想法。
有了 R，可以確定從 PHP 數學庫獲得的值與那些從成熟的免費可用的開放源碼統計包中獲得的值是否一致。
　　清單 4 中的代碼示範了交給 R 來處理以擷取一個值是多麼容易。

清單 4. 交給 R 統計計算包來處理以擷取一個值

// 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 統計值相關的機率值，而第二個反Function Compute了與給定的 alpha 設定相對應的 T 統計值。 getStudentProb 方法用來評估線性模型的吻合程度； getInverseStudentProb 方法返回一個中間值，它用來計算每個預測的 Y值的信賴區間。

由於篇幅有限，我不可能逐個詳細說明這個類中的所有函數，因此如果您想搞清楚簡單線性迴歸分析中所涉及的術語和步驟，我鼓勵您參考大學本科學生使用的統計學教科書。

燃耗研究

要示範如何使用該類，我可以使用來自公用事業中燃耗（burnout）研究中的資料。Michael Leiter 和 Kimberly Ann Meechan 研究了稱為 消耗指數（Exhaustion Index)的燃耗度量單位和稱之為 集中度（Concentration）的獨立變數之間的關係。集中度是指人們的社交接觸中來自其工作環境的那部分比例。

要研究他們樣本中個人的消耗指數值與集中度值之間的關係，請將這些值裝入適當命名的數組中，並用這些數組值對該類進行執行個體化。對類進行執行個體化後，顯示該類所產生的某些匯總值以評估線性模型與資料的吻合程度。

清單 5 顯示了裝入資料和顯示匯總值的指令碼：

清單 5. 用於裝入資料並顯示匯總值的指令碼

// 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);　　
?>

　 　 　 　 　 　

Equation:
T:
Prob > T:

通過 網頁瀏覽器運行該指令碼，產生以下輸出：

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

這張表的最後一行指出擷取這樣大 T值的 隨機機率非常低。可以得出這樣的結論：與僅僅使用消耗值的均值相比，簡單線性模型的預測能力更好。

知道了某個人的工作場所聯絡的集中度，就可以用來預測他們可能正在消耗的燃耗程度。這個方程告訴我們：集中度值每增加 1 個單位，社會服務領域中一個人的消耗值就會增加 8 個單位。這進一步證明了：要減少潛在的燃耗，社會服務領域中的個人應當考慮在其工作場所之外結交朋友。

這隻是粗略地描述了這些結果可能表示的含義。為全面研究這個資料集的含義，您可能想更詳細地研究這個資料以確信這是正確的解釋。在下一篇文章中我將討論應當執行其它哪些分析。

您學到了什嗎？

其一，要開發意義重大的基於 PHP 的數學包，您不必是一名火箭科學家。堅持標準的物件導向技術，以及明確地採用逆向連結問題解決方案，就可以相對方便地使用 PHP 實現某些較為基本的統計過程。

從教學的觀點出發，我認為：如果只是因為要求您在較高和較低的抽象層次思考統計測試或常式，那麼這個練習是非常有用的。換句話說，補充您的統計測試或過程學習的一個好辦法就是將這個過程作為演算法實現。

要實現統計測試通常需要超出所給定的資訊範圍並創造性地解決和發現問題。對於發現對某個學科認識的不足而言，它也是一個好辦法。

不利的一面，您發現 PHP 對於取樣分布缺乏內在手段，而這是實現大多數統計測試所必需的。您需要交給 R 來處理以擷取這些值，但是我擔心您會沒時間或沒興趣安裝 R。某些常見機率函數的本機 PHP 實現可以解決這個問題。

另一個問題是：該類產生許多中間值和匯總值，但是匯總輸出實際上沒有利用這一點。我提供了一些難處理的輸出，但是這既不夠充分也沒進行很好的組織，以致您無法充分地解釋分析結果。實際上，我完全不知道如何可以將輸出方法整合到該類中。這需要得到解決。

最後，要弄明白資料，不僅僅是察看匯總值就可以了。您還需要明白各個資料點是如何分布的。最好的辦法之一是將您的資料繪製成圖表。再次聲明，我對這方面不太瞭解，但是如果要用這個類來分析實際資料的話就需要解決這個問題。

