概念
簡單線性迴歸建模背後的基本目標是從成對的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) { 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 值。可以用它們來評估線性模型與資料的吻合程度。高的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/ |
將來,PHP_MATH 變數的設定將透過一個用於整個PHP 數學庫的設定檔來完成。
您學到了什麼?
在本文中,您了解如何使用SimpleLinearRegression 類別開發用於中小規模的資料集的資料研究工具。在這個過程中,我還開發了一個供SimpleLinearRegression 類別使用的本機機率函數,並用HTML 輸出方法和基於JpGraph 庫的圖形生成程式碼擴展該類別。
從學習的角度來看,簡單線性迴歸建模是值得進一步研究的,因為事實證明,它是理解更高級形式的統計建模的必經之路。在深入學習更高級的技術(如多次迴歸或多變量變異數分析)之前,對於簡單線性迴歸的透徹理解將使您受益匪淺。
即使簡單線性迴歸只用一個變數來說明或預測另一個變數的偏離值,在所有的研究變數之間尋找簡單線性關係仍然常常是研究性資料分析的第一步。僅僅因為數據是多元的並不意味著就必須使用多元工具來研究它。實際上,在開始時使用簡單線性迴歸這樣的基本工具是著手探究資料模式的好方法。