PHP實現簡單線性回歸之數據研究工具
概念
簡單線性回歸建模背后的基本目標是從成對的 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;;; };
}; ?>
輸出方法
既然您已經用 PHP 實現了概率函數,那么開發基于 PHP 的數據研究工具剩下的唯一難題就是設計用于顯示分析結果的方法。
簡單的解決方案是根據需要將所有實例變量的值都顯示到屏幕上。在第一篇文章中,當顯示燃耗研究(Burnout Study)的線性方程、 T值和 T 概率時,我就是這么做的。能根據特定目的而訪問特定值是很有幫助的, SimpleLinearRegression 支持此類用法。
然而,另一種用于輸出結果的方法是將輸出的各部分系統化地進行分組。如果研究用于回歸分析的主要統計軟件包的輸出,就會發現它們往往是用同樣的方式對輸出進行分組的。它們往往有 摘要表(Summary Table)、 偏離值分析(Analysis Of Variance)表、 參數估計值(Parameter Estimate)表和 R 值(R Value)。類似地,我創建了一些輸出方法,名稱如下:
showSummaryTable() showAnalysisOfVariance() showParameterEstimates() showRValues() 我還有一個用于顯示線性預測公式的方法( getFormula() )。許多統計軟件包不輸出公式,而是希望用戶根據上述方法的輸出構造公式。部分是由于您最后用來對數據建模的公式的最終形式可能由于下列原因而與缺省公式不同:
Y軸截距沒有有意義的解釋,或者 輸入值可能是經過轉換的,而您可能需要取消對它們的轉換以獲取最終的解釋。
所有這些方法都假定輸出媒介是網頁。考慮到您有可能希望用非網頁的其它媒介輸出這些匯總值,所以我決定將這些輸出方法包裝在一個繼承了 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='http://www.cgvv.com.cn/bcjs/<?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.phpSimpleLinearRegressionHTML.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 庫的圖形生成代碼擴展該類。 從學習的角度來看,簡單線性回歸建模是值得進一步研究的,因為事實證明,它是理解更高級形式的統計建模的必由之路。在深入學習更高級的技術(如多次回歸或多變量方差分析)之前,對于簡單線性回歸的透徹理解將使您受益匪淺。即使簡單線性回歸只用一個變量來說明或預測另一個變量的偏離值,在所有的研究變量之間尋找簡單線性關系仍然常常是研究性數據分析的第一步。僅因為數據是多元的并不意味著就必須使用多元工具研究它。實際上,在開始時使用簡單線性回歸這樣的基本工具是著手探究數據模式的好方法。