用 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值的觀測值和預測值最小。觀測值和預測值之差稱為誤差（ y i- (mx i+ 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 的概率函數需要些什麼。

我開始上網查找信息和代碼。一個兩者兼有的來源是書籍 Numerical Recipes in C 中的概率函數。我用 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 } } ?>
 

回頁首

圖形輸出

迄今為止，您已經實現的輸出方法都以 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('nav