Machine Learning With Go 第4章:回歸

4 回歸

之前有轉載過一篇文章:容量推薦引擎:基於吞吐量和利用率的預測縮放,裏面用到了基本的線性回歸來預測容器的資源利用情況。後面打算學一下相關的知識,譯自:Machine Learning With Go

我們將探究的第一組機器學習技術通常被稱為回歸(regression),我們可以將回歸理解為一個變量(例如銷售額)的變化是如何影響到其他變量(如用戶數)的。對於機器學習技術來說,這是一個很好的開端,它們是構成其他更加複雜技術的基礎。

機器學習中的回歸技術通常會注重評估連續值(如股票價格、溫度或疾病進展等)。下一章討論的歸類(Classification)則會注重離散值,或離散的類集合(如欺詐/非欺詐、坐下/起立/跑動,或熱狗/非熱狗等)。正如上面提到的,回歸技術會貫徹到機器學習中,並作為歸類算法的一部分,但本章中,我們將會注重其基本的應用–預測連續值。

理解回歸模型術語

正如前面提到的,回歸本身是一個分析一個變量和另一個變量之間關係的過程,但在機器學習中還用到了一些術語來描述這些變量以及各種類型的回歸和與回歸有關的過程:

  • 響應(response)或因變量(dependent variable):這些術語可以互用,表示基於其他一個或多個變量來預測的變量,通常使用y表示
  • 解釋變量(Explanatory variables)、自變量(independent variables)、特徵(features)、屬性(attributes)或回歸係數(regressors):這些術語可以互用,表示用於預測響應的變量,通常使用xx1 , x2表示
  • 線性回歸:該類型的回歸會假設因變量會線性依賴自變量(即遵循直線方程)
  • 非線性回歸:該類型的回歸會假設因變量會非線性依賴自變量(如多項式或指數)
  • 多元回歸(Multiple regression:):具有超過一個自變量的回歸
  • 擬合(Fitting)或訓練(training):參數化模型的過程(如回歸模型),可以用該模型來預測一個特定的因變量
  • 預測:使用參數模型預測因變量的過程(如回歸模型)

部分術語會在回歸上下文和本書的其他上下文中使用。

線性回歸

線性回歸是最簡單的機器學習模型之一,但不能出於某些原因而忽略該模型。正如前面提到的,它是其他模型的基礎,且有一些非常重要的優勢。

正如在本書中討論的,完整性在機器學習應用中非常重要,模型越簡單,解釋性越強,則越容易維護其完整性。此外,如果模型簡單且具有解釋性,那麼就可以幫助理解變量之間的推斷關係,並簡化開發過程中的檢查工作。

Mike Lee Williams (來自Fast Forward Labs的)說過:

未來是算法,可解釋模型在人類和智能機器之間建立了一種更安全、更高效、最終更具協作性的關係。

線性回歸模型是可解釋的,因此可以為數據科學提供一種安全且高效的選項。當需要搜索一種可以預測連續變量的模型時,如果數據和相關條件具備,則應該考慮並使用線性回歸(或多元線性回歸)。

線性回歸概述

在線性回歸中,我們會嘗試使用如下線性方程,使用一個自變量x,對因變量y進行建模:

\[y = mx + b
\]

這裡,m為直線的斜率,b為截距。例如,我們想要通過每天訪問網站的userssales進行建模,為了使用線性回歸,我們會希望通過確定mb來讓預測公司的銷售額:

\[sales = m * (number~of~users) + b
\]

這樣,我們的訓練模型就是該參數化函數。通過輸入Number of Users 來預測 sales,如下:

image

線性回歸的訓練或擬合需要確定mb的值,這樣得出的公式就有預測響應的能力。有多種方式來確定mb,但最常見的是普通最小二乘法(ordinary least squares (OLS))。

為了使用OLS來確定mb,首先為mb選擇一個值來創建第一條示例線(example line)。然後測量每個已知點(如訓練集)和示例線之間的垂直距離,這些距離稱為誤差(errors)或殘差(residuals)。下圖展示了評估和驗證:

image

下面,我們計算這些誤差平方和:

\[\frac{error^2_1+error^2_2+…+error^2_N}{N}
\]

通過調整mb來最小化誤差的平方和。換句話說,我們訓練的線性回歸直線是平方和最小的直線。

有很多種方式可以找出誤差平方和最小的直線,如通過OLS可以找出並分析這條直線。但最常用的減少誤差平方和的優化技術稱為梯度下降法(gradient descent)。相比於分析法,這種方法更容易實現,且便於計算(如內存),也更加靈活。

可以說,線性回歸和其他回歸的實現都利用梯度下降來擬合或訓練線性回歸線。實際上,梯度下降法在機器學習中無處不在,由此可以產生更加複雜的模型技術,如深度學習。

梯度下降法

梯度下降法有很多變種,且在機器學習世界中無處不在。最重要的是,它們用於確定線性或邏輯回歸等算法的最佳係數,同時也在更複雜的技術中發揮着重要作用(至少部分基於線性/邏輯回歸(如神經網絡))。

梯度下降法的一般思想是確定某些參數的變化方向和幅度,這些參數將使預測曲線朝着正確的方向移動,以優化某些度量(如誤差)。想像站在某個地方,如果要向較低的位置移動,則需要朝向下的方向移動。這基本上就是梯度下降算法在優化參數時所做的事情。

讓我們看一下所謂的隨機梯度下降(SGD),這是一種增量的梯度下降,從而對這個過程有更多直覺上的了解。我們在第5章”分類”的邏輯回歸實現中使用了SGD。在該示例中,我們實現了對邏輯回歸參數的訓練或擬合,如下所示:

// logisticRegression fits a logistic regression model
// for the given data.
func logisticRegression(features *mat64.Dense, labels []float64, numSteps int, learningRate f)
// Initialize random weights.
_, numWeights := features.Dims()
weights := make([]float64, numWeights)
s := rand.NewSource(time.Now().UnixNano())
r := rand.New(s)
for idx, _ := range weights {
  weights[idx] = r.Float64()
}
// Iteratively optimize the weights.
for i := 0; i < numSteps; i++ {
  // Initialize a variable to accumulate error for this iteration.
  var sumError float64
  // Make predictions for each label and accumulate error.
  for idx, label := range labels {
      // Get the features corresponding to this label.
      featureRow := mat64.Row(nil, idx, features)
      // Calculate the error for this iteration's weights.
      pred := logistic(featureRow[0]*weights[0] + featureRow[1]*weights[1])
      predError := label - pred
      sumError += math.Pow(predError, 2)
      // Update the feature weights.
      for j := 0; j < len(featureRow); j++ {
          weights[j] += learningRate * predError * pred * (1 - pred) * featureRow[j]
      }
  }
}
return weights
}

// Iteratively optimize the weights注釋下面的循環實現了通過SGD來優化邏輯回歸參數。下面選擇這部分循環來看下到底發生了什麼。

首先,我們使用當前權重和預測值與理想值(即實際觀察值)之間的差值來計算模型的輸出:

// Calculate the error for this iteration's weights.
pred := logistic(featureRow[0]*weights[0] + featureRow[1]*weights[1])
predError := label - pred
sumError += math.Pow(predError, 2)

根據SGD,我們將根據如下公式來計算參數(在本例中為權weights)的更新:

\[update=leaning~rate\times~gradient~of~the~parameters
\]

gradient是cost函數的數學梯度。

更多參見: //mathworld.wolfram.com/Gradient.html

然後將該更新應用到參數,如下所示:

\[parameter=parameters-update
\]

在我們的邏輯回歸模型中,計算結果如下:

// Update the feature weights.
for j := 0; j < len(featureRow); j++ {
  weights[j] += learningRate * predError * pred * (1 - pred) * featureRow[j]
}

機器學習中廣泛使用了這種類型的SGD,但在某些場景下,這種梯度下降法可能導致過擬合或陷入局部最小值/最大值(而不是尋找全局最優值)。

為了解決這些問題,可以使用一個梯度下降的變種,稱為批量梯度下降(batch gradient descent)。在批量梯度下降中,可以基於所有訓練數據集中的梯度來計算每個參數更新,而不針對數據集的特定觀測值或行。這種方式有助於防止過度擬合,但它也可能很慢,並且存在內存問題,因為需要計算每個參數相對於整個數據集的梯度。微批量梯度下降(Mini-batch gradient descent)是另一個變種,它在試圖保持批量梯度下降的某些好處的同時,更易於計算。在微批量梯度下降法中,梯度是在訓練數據集的子集上計算的,而不是在整個訓練數據集上計算的。

在邏輯回歸的場景中,你可能看到過使用梯度上升或下降,梯度上升與梯度下降是一回事,只是cost函數的方向不同而已。更多參見: //stats.stackexchange.com/questions/261573/using-gradient-ascent-instead-of-gradient-descent-for-logistic-regression

gonum團隊已經實現了梯度下降法:gonum.org/v1/gonum/optimize。文檔地址://pkg.go.dev/gonum.org/v1/gonum/optimize#GradientDescent

線性回歸的假設和缺點

與所有機器學習模型一樣,線性回歸併不能適用於所有場景,它的前提是假設你的數據之間的關係是確定的:

  • 線性關係:線性回歸會假設因變量線性依賴自變量(線性方程)。如果這種關係不是線性的,則線性回歸可能會表項不佳
  • 正態性:假設變量遵循正太分佈(看起來像鐘形)。本章後面會討論這種特性以及非正態分佈變量下的取捨。
  • 非多重共線性:多重共線性是一個特別的術語,它意味着自變量並不是真正獨立的,它們會以某種形式相互依賴
  • 沒有自相關性:自相關性是另一個特別的術語,意味着變量依賴自身或自身的某個版本(如存在某些可預測的時序中)。
  • 同方差性:這可能是這一組術語中最特別的一個,但它相對比較簡單,且並不需要經常關注。線性回歸假設回歸線周圍的數據的方差與自變量值的方差大致相同。

從技術上講,為了使用線性回歸,需要滿足上述所有假設。但最重要的是我們需要知道數據是如何分佈的,以及它們是如何表現的。後續在使用線性回歸的示例中會深入討論這些假設。

作為一個數據科學家或分析師,在使用線性回歸時需要注意到線性回歸的不足:

  • 使用特定範圍的自變量來訓練線性回歸模型,在預測該範圍外的數據時應該格外小心,因為你的線性回歸直線可能並不適用(如,某些極端數值下,因變量可能並不是線性的)。
  • 可能為兩個並無關聯的變量建立了一個線性回歸模型。需要確保變量之間有邏輯上的關聯性。
  • 可能會因為擬合某些特定類型數據中的異常或極端值而偏離回歸線,如OLS。有一些方式可以讓擬合回歸不受異常值的影響,或針對異常值展示出不同的行為,如正交最小二乘法(orthogonal least squares)或嶺回歸(ridge regression)。

線性回歸的例子

為了描述線性回歸,讓我們創建第一個機器學習模型。下面數據為廣告數據,保存格式為.csv

$ head Advertising.csv
TV,Radio,Newspaper,Sales
230.1,37.8,69.2,22.1
44.5,39.3,45.1,10.4
17.2,45.9,69.3,9.3
151.5,41.3,58.5,18.5
180.8,10.8,58.4,12.9
8.7,48.9,75,7.2
57.5,32.8,23.5,11.8
120.2,19.6,11.6,13.2
8.6,2.1,1,4.8

該數據集包括一系列表示廣告媒體屬性(TVRadioNewspaper)以及對應的銷售額(Sales),本例中我們的目標是對銷售額(因變量)和廣告支出(因變量)進行模型。

分析數據

為了構建模型(或流程),並確保能夠對模型的結果進行檢查,首先需要對數據進行分析(所有機器學習模型的第一個步驟)。我們需要了解變量是如何分佈的,以及變量的範圍和可變性。

為了實現該目標,我們將計算第2章矩陣、概率和統計中討論的匯總數據。這裡,我們將使用github.com/go-gota/gota/tree/master/dataframe中的內置方法,一次性計算出數據集中的所有列的匯總信息:

// Open the CSV file.
advertFile, err := os.Open("Advertising.csv")
if err != nil {
    log.Fatal(err)
}
defer advertFile.Close()
// Create a dataframe from the CSV file.
advertDF := dataframe.ReadCSV(advertFile)
// Use the Describe method to calculate summary statistics
// for all of the columns in one shot.
advertSummary := advertDF.Describe()
// Output the summary statistics to stdout.
fmt.Println(advertSummary)

編譯並運行後得到如下結果:

$ go build
$ ./myprogram
[7x5] DataFrame
   column   TV         Radio     Newspaper  Sales
0: mean     147.042500 23.264000 30.554000  14.022500
1: stddev   85.854236  14.846809 21.778621  5.217457
2: min      0.700000   0.000000  0.300000   1.600000
3: 25%      73.400000  9.900000  12.600000  10.300000
4: 50%      149.700000 22.500000 25.600000  12.900000
5: 75%      218.500000 36.500000 45.100000  17.400000
6: max      296.400000 49.600000 114.000000 27.000000
   <string> <float>    <float>   <float>    <float>

上面以表格形式打印出所有的匯總數據,包括平均值、標準偏差、最小值、最大值、25%/75%百分位和中位數(或50%百分位)。

這些值為我們提供了良好的數值參考,後續會在訓練線性回歸模型時將看到這些數字。但缺乏直觀上的理解,為此,我們需要為每列數值創建一個直方圖:

// Open the advertising dataset file.
f, err := os.Open("Advertising.csv")
if err != nil {
    log.Fatal(err)
}
defer f.Close()
// Create a dataframe from the CSV file.
advertDF := dataframe.ReadCSV(f)
// Create a histogram for each of the columns in the dataset.
for _, colName := range advertDF.Names() {
    // Create a plotter.Values value and fill it with the
    // values from the respective column of the dataframe.
    plotVals := make(plotter.Values, advertDF.Nrow())
    for i, floatVal := range advertDF.Col(colName).Float() {
        plotVals[i] = floatVal
    }
    // Make a plot and set its title.
    p, err := plot.New()
    if err != nil {
        log.Fatal(err)
    }
    p.Title.Text = fmt.Sprintf("Histogram of a %s", colName)
    // Create a histogram of our values drawn
    // from the standard normal.
    h, err := plotter.NewHist(plotVals, 16)
    if err != nil {
        log.Fatal(err)
    }
    // Normalize the histogram.
    h.Normalize(1)
    // Add the histogram to the plot.
    p.Add(h)
    // Save the plot to a PNG file.
    if err := p.Save(4*vg.Inch, 4*vg.Inch, colName+"_hist.png"); err != nil {
        log.Fatal(err)
    }
}

本程序會為每個直方圖創建一個.png圖像:

image

觀察上圖以及計算出的匯總信息,下一步考慮是否符合線性回歸的假設條件。可以看到並不是所有的變量都是正態分佈的(鐘形的)。可以看到銷售額是鐘形的,而其他則不是正態的。

我們可以使用分位圖(quantile-quantile (q-q) p)統計工具來確定分佈與正態分佈的接近程度,甚至通過統計測試來確定變量是否服從正態分佈的概率。但大多數情況下,通過圖表就可以得出一個大致的結論。

下一步要做出決策,但至少有一部分數據在技術上並不會擬合到我們的線性回歸模型中,可以選擇如下一種方式進行處理:

  • 嘗試轉換變量,使其遵循正態分佈,並在線性回歸模型中使用這些轉換的變量。這種方式的好處是可以在模型的假設中進行操作,缺點是可能會讓模型難以理解,降低可解釋性
  • 使用不同的數據來解決問題
  • 在線性回歸假設中忽略該問題,並嘗試創建該模型

可能還有其他解決問題的方式,但我的建議是首先嘗試第三種選項。由於可以快速地訓練線性回歸模型,因此該選項並不會帶來多少壞處。如果最後得出滿意的模型,那麼就可以避免引入更多的複雜性。如果得到的模型不盡如人意,那麼此時再訴諸於其他選項。

選擇自變量

現在對我們的數據有了一些直覺上的了解,並且已經了解到數據是如何擬合線性回歸模型的假設的。那麼現在應該選擇哪個變量作為我們的自變量來預測因變量?

最簡單的方法是通過直觀地探索因變量和選擇的所有自變量之間的相關性,特別是可以通過繪製因變量與其他每個變量的散點圖(使用pkg.go.dev/gonum.org/v1/plot)來做決定:

// Open the advertising dataset file.
f, err := os.Open("Advertising.csv")
if err != nil {
    log.Fatal(err)
}
defer f.Close()
// Create a dataframe from the CSV file.
advertDF := dataframe.ReadCSV(f)
// Extract the target column.
yVals := advertDF.Col("Sales").Float()
// Create a scatter plot for each of the features in the dataset.
for _, colName := range advertDF.Names() {
    // pts will hold the values for plotting
    pts := make(plotter.XYs, advertDF.Nrow())
    // Fill pts with data.
    for i, floatVal := range advertDF.Col(colName).Float() {
        pts[i].X = floatVal
        pts[i].Y = yVals[i]
    }
    // Create the plot.
    p, err := plot.New()
    if err != nil {
        log.Fatal(err)
    }
        p.X.Label.Text = colName
    p.Y.Label.Text = "y"
    p.Add(plotter.NewGrid())
    s, err := plotter.NewScatter(pts)
    if err != nil {
        log.Fatal(err)
    }
    s.GlyphStyle.Radius = vg.Points(3)
    // Save the plot to a PNG file.
    p.Add(s)
    if err := p.Save(4*vg.Inch, 4*vg.Inch, colName+"_scatter.png"); err != nil {
        log.Fatal(err)
    }
}

如此可以創建如下散點圖:

image

通過這些散點圖,我們需要推斷出哪些屬性 (TV, Radio, 和/或 Newspaper)與我們的因變量(Sales)具有線性關係。是否可以在這些散點圖上畫一條線,以符合銷售趨勢和各自的屬性?這種方法並不總是行得通,且對於一個特定的問題,並不一定可以將其關聯到所有的屬性。

上述場景中,RadioTVSales呈線性關係,Newspaper可能與Sales有一定的關係,但相關性並不明顯。與TV的相關性是最明顯的,因此先選擇TV作為線性回歸模型的自變量,線性回歸公式如下:

\[Sales = m~TV+b
\]

這裡要注意的另一件事是,變量TV可能不是嚴格等方差的(在線性回歸的假設中討論過)。這一點需要注意(可能值得在項目中歸檔的),下面將繼續探究是否可以創建具有預測能力的線性回歸模型。當模型表現不佳時,需要重新審視這種假設。

創建訓練和測試集

為了避免過度擬合併保證模型的推廣,我們需要將數據集劃分為訓練集和測試集即評估和驗證(Evaluation and Validation)。這裡我們不會聚焦某個測試集,因為只需要進行一次模型訓練即可,而不會在訓練和測試之間來回迭代。但如果需要多個因變量進行驗證和/或需要迭代調整模型參數時,你可能希望創建一個保留集,保存到模型開發過程結束後進行驗證。

我們將使用github.com/go-gota/gota/blob/master/dataframe創建訓練和測試數據集,並將它們保存到各自的.csv文件中。該場景中,我們使用80/20的比例來劃分訓練和測試數據:

// Open the advertising dataset file.
f, err := os.Open("Advertising.csv")
if err != nil {
    log.Fatal(err)
}
defer f.Close()
// Create a dataframe from the CSV file.
// The types of the columns will be inferred.
advertDF := dataframe.ReadCSV(f)
// Calculate the number of elements in each set.
trainingNum := (4 * advertDF.Nrow()) / 5
testNum := advertDF.Nrow() / 5
if trainingNum+testNum < advertDF.Nrow() {
    trainingNum++
}
// Create the subset indices.
trainingIdx := make([]int, trainingNum)
testIdx := make([]int, testNum)
// Enumerate the training indices.
for i := 0; i < trainingNum; i++ {
    trainingIdx[i] = i
}
// Enumerate the test indices.
for i := 0; i < testNum; i++ {
    testIdx[i] = trainingNum + i
}
// Create the subset dataframes.
trainingDF := advertDF.Subset(trainingIdx)
testDF := advertDF.Subset(testIdx)
// Create a map that will be used in writing the data
// to files.
setMap := map[int]dataframe.DataFrame{
    0: trainingDF,
    1: testDF,
}
// Create the respective files.
for idx, setName := range []string{"training.csv", "test.csv"} {
    // Save the filtered dataset file.
    f, err := os.Create(setName)
    if err != nil {
        log.Fatal(err)
    }
    // Create a buffered writer.
    w := bufio.NewWriter(f)
    // Write the dataframe out as a CSV.
    if err := setMap[idx].WriteCSV(w); err != nil {
        log.Fatal(err)
    }
}

上述代碼會輸出如下訓練和測試集:

$ wc -l *.csv
    201 Advertising.csv
    41  test.csv
    161 training.csv
    403 total

這裡使用的數據並沒有經過排序。但如果需要按照響應、日期或其他方式處理數據,則最好隨機劃分訓練和測試集。如果不這麼做,訓練和測試集可能會包含特定範圍的響應,這樣響應可能會受到時間/日期等人為因素的影響。

訓練模型

下面將訓練(或擬合)我們的線性回歸模型。這也意味着需要找到誤差平方和最小的的斜率(m)和截距(b)。為了執行訓練,我們會使用來自Sajari的包:github.com/sajari/regressionSajari是一個重度依賴Go和機器學習的網站搜索公司,他們在生產中使用了github.com/sajari/regression

為了使用github.com/sajari/regression來訓練回歸模型,需要初始化一個regression.Regression值,並設置一對標籤,然後使用被標記的訓練數據來填充regression.Regression。之後就可以簡單地使用 調用Run()來對regression.Regression的值進行訓練,以此生成線性回歸模型。

// Open the training dataset file.
f, err := os.Open("training.csv")
    if err != nil {
    log.Fatal(err)
}
defer f.Close()
// Create a new CSV reader reading from the opened file.
reader := csv.NewReader(f)
// Read in all of the CSV records
reader.FieldsPerRecord = 4
trainingData, err := reader.ReadAll()
if err != nil {
    log.Fatal(err)
}
// In this case we are going to try and model our Sales (y)
// by the TV feature plus an intercept. As such, let's create
// the struct needed to train a model using github.com/sajari/regression.
var r regression.Regression
r.SetObserved("Sales")
r.SetVar(0, "TV")
// Loop of records in the CSV, adding the training data to the regression value.
for i, record := range trainingData {
    // Skip the header.
    if i == 0 {
        continue
    }
    // Parse the Sales regression measure, or "y".
    yVal, err := strconv.ParseFloat(record[3], 64)
    if err != nil {
        log.Fatal(err)
    }
    // Parse the TV value.
    tvVal, err := strconv.ParseFloat(record[0], 64)
    if err != nil {
        log.Fatal(err)
    }
    // Add these points to the regression value.
    r.Train(regression.DataPoint(yVal, []float64{tvVal}))
}
// Train/fit the regression model.
r.Run()
// Output the trained model parameters.
fmt.Printf("\nRegression Formula:\n%v\n\n", r.Formula)

編譯執行上述代碼,訓練的線性回歸公式會打印到標準輸出:

$ go build
$ ./myprogram

Regression Formula:
Predicted = 7.07 + TV*0.05

從結果可以看到,引用的包輸出的線性回歸的截距為7.07,斜率為0.5。這裡可以進行一些簡單的檢查,因為我們在散點圖中看到了TVSales之間的相關性是上升和向右的(即正相關),這也意味着公式的斜率應該是正數。

評估訓練模型

下面需要通過評估模型的表現來查看是否可以使用自變量TV來預測Sales。為此,需要加載測試集,使用訓練過的模型對每個測試例進行預測,然後計算第3章”評估和驗證”中討論的某個評估指標。

為此,我們使用平均絕對誤差(Mean Absolute Error (MAE))作為評估指標,這樣就可以直接對比結果和Sales,而不必太擔心異常值或極端值。

為了通過訓練的regression.Regression 值來預測Sales,只需解析測試集的值,並針對regression.Regression 的值調用Predict()。然後,計算預測值與觀測值的差值,得到差值的絕對值,然後將所有絕對值相加,得到MAE

// Open the test dataset file.
f, err = os.Open("test.csv")
if err != nil {
    log.Fatal(err)
}
defer f.Close()
// Create a CSV reader reading from the opened file.
reader = csv.NewReader(f)
// Read in all of the CSV records
reader.FieldsPerRecord = 4
testData, err := reader.ReadAll()
if err != nil {
    log.Fatal(err)
}
// Loop over the test data predicting y and evaluating the prediction
// with the mean absolute error.
var mAE float64
for i, record := range testData {
    // Skip the header.
    if i == 0 {
    continue
    }
    // Parse the observed Sales, or "y".
    yObserved, err := strconv.ParseFloat(record[3], 64)
    if err != nil {
        log.Fatal(err)
    }
    // Parse the TV value.
    tvVal, err := strconv.ParseFloat(record[0], 64)
    if err != nil {
        log.Fatal(err)
    }
    // Predict y with our trained model.
    yPredicted, err := r.Predict([]float64{tvVal})
    // Add the to the mean absolute error.
    mAE += math.Abs(yObserved-yPredicted) / float64(len(testData))
}
// Output the MAE to standard out.
fmt.Printf("MAE = %0.2f\n\n", mAE)

編譯並運行評估程序,得到如下結果:

$ go build
$ ./myprogram

Regression Formula:
Predicted = 7.07 + TV*0.05

MAE = 3.01

那麼如果MAE = 3.01,我們怎麼知道該值是好的還是壞的?這也是為什麼擁有一個良好的數據心智模型很重要的原因。我們已經計算了銷售額的平均值、範圍和標準差。平均銷售額為14.02,標準差為5.21。這樣我們的MAE小於銷售額數值的標準差,約為平均值的20%,說明我們的模型具有一定的預測能力。

恭喜,我們構建了第一個具有預測能力的機器學習模型。

為了更直觀地了解模型的運行狀況,可以藉助圖形來幫助可視化線性回歸線(利用gonum.org/v1/plot)。首先,創建一個可以執行預測的函數(而需要使用 github.com/sajari/regression包,相當於樁函數),通過這種方式可以提供輕量級內存訓練模型:

// predict uses our trained regression model to made a prediction.
func predict(tv float64) float64 {
    return 7.07 + tv*0.05
}

然後創建可視化的回歸線:

// Open the advertising dataset file.
f, err := os.Open("Advertising.csv")
if err != nil {
    log.Fatal(err)
}
defer f.Close()
// Create a dataframe from the CSV file.
advertDF := dataframe.ReadCSV(f)
// Extract the target column.
yVals := advertDF.Col("Sales").Float()
// pts will hold the values for plotting.
pts := make(plotter.XYs, advertDF.Nrow())
// ptsPred will hold the predicted values for plotting.
ptsPred := make(plotter.XYs, advertDF.Nrow())
// Fill pts with data.
for i, floatVal := range advertDF.Col("TV").Float() {
    pts[i].X = floatVal
    pts[i].Y = yVals[i]
    ptsPred[i].X = floatVal
    ptsPred[i].Y = predict(floatVal)
}
// Create the plot.
p, err := plot.New()
if err != nil {
    log.Fatal(err)
}
p.X.Label.Text = "TV"
p.Y.Label.Text = "Sales"
p.Add(plotter.NewGrid())
// Add the scatter plot points for the observations.
s, err := plotter.NewScatter(pts)
if err != nil {
    log.Fatal(err)
}
s.GlyphStyle.Radius = vg.Points(3)
// Add the line plot points for the predictions.
l, err := plotter.NewLine(ptsPred)
if err != nil {
    log.Fatal(err)
}
l.LineStyle.Width = vg.Points(1)
l.LineStyle.Dashes = []vg.Length{vg.Points(5), vg.Points(5)}
// Save the plot to a PNG file.
p.Add(s, l)
if err := p.Save(4*vg.Inch, 4*vg.Inch, "regression_line.png"); err != nil {
    log.Fatal(err)
}

編譯和運行後得到如下圖:

image

可以看到,我們訓練的回歸線與實際的數據趨勢相匹配。

多元線性回歸

線性回歸併不局限於依賴單個自變量的簡單線性公式。多元線性回歸與前面討論的類似,但具有多個自變量(x1,x2等)。這種場景下的直線方程如下:

\[y=m_1x_1+m_1x_2+…+m_Nx_N+b
\]

這裡x作為自變量,m作為與自變量相關的斜率,此外還有一個截距b

多元線性回歸相對比較難以可視化和思考,因為它不再是一條可以在二維中可視化的直線。而是一條在二維、三維或多維的線性曲面。但它使用了很多在一元線性回歸中用過的技術。

多元線性回歸具有與一元線性回歸相同的假設,但需要注意的是與之相關的陷阱:

  • 過擬合:通過為模型添加越來越多的自變量,會增加模型的複雜度,並存在過擬合的風險。可以使用之前推薦的技術:正則化(regularization)來解決這種問題。 正則化在模型中創建一個懲罰項,該懲罰項是一個與模型複雜性有關的函數,有助於控制這種影響。
  • 相對比例(Relative Scale):在某些場景下,如果其中某個自變量的比例與另一個自變量的比例相差幾個數量級,那麼較大的變量可能會抵消較小變量帶來的影響,因此可能需要考慮規範化變量。

記住以上兩點,下面嘗試將Sales模型從一元線性回歸模型擴展到多元線性回歸模型。回顧一下前面章節中的散點圖,可以看到Radio似乎也與Sales呈線性關係,因此可以嘗試創建一個多元線性回歸模型,如下:

\[Sales=m_1TV+m_2Radio+b
\]

使用github.com/sajari/regression時,需要在regression.Regression中標記其他變量,並確保這些值在訓練數據點中成對出現:

// Open the training dataset file.
f, err := os.Open("training.csv")
if err != nil {
    log.Fatal(err)
}
defer f.Close()
// Create a new CSV reader reading from the opened file.
reader := csv.NewReader(f)
// Read in all of the CSV records
reader.FieldsPerRecord = 4
trainingData, err := reader.ReadAll()
if err != nil {
    log.Fatal(err)
}
// In this case we are going to try and model our Sales
// by the TV and Radio features plus an intercept.
var r regression.Regression
r.SetObserved("Sales")
r.SetVar(0, "TV")
r.SetVar(1, "Radio")
// Loop over the CSV records adding the training data.
for i, record := range trainingData {
    // Skip the header.
    if i == 0 {
        continue
    }
    // Parse the Sales.
    yVal, err := strconv.ParseFloat(record[3], 64)
    if err != nil {
        log.Fatal(err)
    }
    // Parse the TV value.
    tvVal, err := strconv.ParseFloat(record[0], 64)
    if err != nil {
        log.Fatal(err)
    }
    // Parse the Radio value.
    radioVal, err := strconv.ParseFloat(record[1], 64)
    if err != nil {
        log.Fatal(err)
    }
    // Add these points to the regression value.
    r.Train(regression.DataPoint(yVal, []float64{tvVal, radioVal}))
}
// Train/fit the regression model.
r.Run()
// Output the trained model parameters.
fmt.Printf("\nRegression Formula:\n%v\n\n", r.Formula)

編譯並運行,得到如下回歸公式:

$ go build
$ ./myprogram

Regression Formula:
Predicted = 2.93 + TV*0.05 + Radio*0.18

可以看到,回歸公式增加了一個額外的自變量項。截距值也發生了變化。

可以用與一元線性回歸模型類似的方式,使用Predict()方法來測試該模型:

// Open the test dataset file.
f, err = os.Open("test.csv")
if err != nil {
    log.Fatal(err)
}
defer f.Close()
// Create a CSV reader reading from the opened file.
reader = csv.NewReader(f)
// Read in all of the CSV records
reader.FieldsPerRecord = 4
testData, err := reader.ReadAll()
if err != nil {
    log.Fatal(err)
}
// Loop over the test data predicting y and evaluating the prediction
// with the mean absolute error.
var mAE float64
for i, record := range testData {
    // Skip the header.
    if i == 0 {
        continue
    }
    // Parse the Sales.
    yObserved, err := strconv.ParseFloat(record[3], 64)
    if err != nil {
        log.Fatal(err)
    }
    // Parse the TV value.
    tvVal, err := strconv.ParseFloat(record[0], 64)
    if err != nil {
        log.Fatal(err)
    }
    // Parse the Radio value.
    radioVal, err := strconv.ParseFloat(record[1], 64)
    if err != nil {
        log.Fatal(err)
    }
    // Predict y with our trained model.
    yPredicted, err := r.Predict([]float64{tvVal, radioVal})
    // Add the to the mean absolute error.
    mAE += math.Abs(yObserved-yPredicted) / float64(len(testData))
}
// Output the MAE to standard out.
fmt.Printf("MAE = %0.2f\n\n", mAE)

運行該程序將為我們的新多元回歸模型得出如下MAE

$ go build
$ ./myprogram

Regression Formula:
Predicted = 2.93 + TV*0.05 + Radio*0.18
MAE = 1.26

新的多元回歸模型提高了MAE值。現在,我們可以根據廣告支出來預測Sales了。你還可以嘗試將Newspaper添加到模型。

注意,模型複雜性增加的同時,也會犧牲掉簡易性,並增加過擬合的風險,因此只考慮當添加的複雜性能夠提升模型的表現、並帶來更大的價值時。

非線性以及其他類型的回歸

雖然本章節主要關注線性回歸,但不會僅限於使用線性方程來執行回歸。你可以使用一個或多個非線性(如冪、指數或其他變換)自變量來為因變量建模。例如,我們可以通過一系列TV項來為Sales建模:

\[Sales=m_1TV+m_2TV^2+m_3TV^3+…+b
\]

注意,增加複雜性的同時也增加了過擬合的風險。

為了實現非線性回歸,不能使用github.com/sajari/regression(僅限於線性回歸),但可以使用github.com/go-hep/hep/tree/main/fit 來擬合或訓練特定的非線性模型。在Go社區中有很多人已經或正在開發非線性模型工具。

除了OLS外還有其他線性回歸技術,可以幫助克服最小二乘線性回歸中的一些假設和弱點。包括嶺回歸和套索回歸(lasso regression)。這兩種技術使用懲罰回歸係數來減輕自變量的多重共線性和非正態性帶來的影響。

github.com/berkmancenter/ridge中實現了Go語言的嶺回歸。與 github.com/sajari/regression不同,我們的自變量和因變量數據是通過gonum矩陣輸入github.com/berkmancenter/ridge的。為了說明該方法,我們首先構造一個包含廣告支出 (TV, Radio, 和Newspaper)的矩陣,以及包含Sales數據的矩陣。注意在github.com/berkmancenter/ridge中,如果想在模型中有一個截距,則需要為截距的輸入自變量矩陣顯式地添加一列,該列中的每個值僅為1.0

// Open the training dataset file.
f, err := os.Open("training.csv")
if err != nil {
    log.Fatal(err)
}
defer f.Close()
// Create a new CSV reader reading from the opened file.
reader := csv.NewReader(f)
reader.FieldsPerRecord = 4
// Read in all of the CSV records
rawCSVData, err := reader.ReadAll()
if err != nil {
    log.Fatal(err)
}
// featureData will hold all the float values that will eventually be
// used to form our matrix of features.
featureData := make([]float64, 4*len(rawCSVData))
yData := make([]float64, len(rawCSVData))
// featureIndex and yIndex will track the current index of the matrix values.
var featureIndex int
var yIndex int
// Sequentially move the rows into a slice of floats.
for idx, record := range rawCSVData {
    // Skip the header row.
    if idx == 0 {
        continue
    }
    // Loop over the float columns.
    for i, val := range record {
        // Convert the value to a float.
        valParsed, err := strconv.ParseFloat(val, 64)
        if err != nil {
            log.Fatal(err)
        }
        if i < 3 {
        // Add an intercept to the model.
        if i == 0 {
            featureData[featureIndex] = 1
            featureIndex++
        }
        // Add the float value to the slice of feature floats.
        featureData[featureIndex] = valParsed
        featureIndex++
        }
        if i == 3 {
            // Add the float value to the slice of y floats.
            yData[yIndex] = valParsed
            yIndex++
        }
    }
}
// Form the matrices that will be input to our regression.
features := mat64.NewDense(len(rawCSVData), 4, featureData)
y := mat64.NewVector(len(rawCSVData), yData)

下面使用自變量和因變量矩陣創建一個新的ridge.RidgeRegression值,然後調用Regress() 方法來訓練模型,最後打印訓練的回歸公式:

// Create a new RidgeRegression value, where 1.0 is the
// penalty value.
r := ridge.New(features, y, 1.0)
// Train our regression model.
r.Regress()
// Print our regression formula.
c1 := r.Coefficients.At(0, 0)
c2 := r.Coefficients.At(1, 0)
c3 := r.Coefficients.At(2, 0)
c4 := r.Coefficients.At(3, 0)
fmt.Printf("\nRegression formula:\n")
fmt.Printf("y = %0.3f + %0.3f TV + %0.3f Radio + %0.3f Newspaper\n\n", c1, c2, c3, c4)

編譯並運行程序,得出如下回歸公式:

$ go build
$ ./myprogram
Regression formula:
y = 3.038 + 0.047 TV + 0.177 Radio + 0.001 Newspaper

可以看到TVRadio的係數與最小二乘回歸得到的結果類似,但略微不同。另外可以看到添加了一個Newspaper項。

可以通過創建自己的預測函數來測試嶺回歸公式:

// predict uses our trained regression model to made a prediction based on a
// TV, Radio, and Newspaper value.
func predict(tv, radio, newspaper float64) float64 {
    return 3.038 + tv*0.047 + 0.177*radio + 0.001*newspaper
}

然後使用該predict函數和測試數據來測試嶺回歸公式:

// Open the test dataset file.
f, err := os.Open("test.csv")
if err != nil {
    log.Fatal(err)
}
defer f.Close()
// Create a new CSV reader reading from the opened file.
reader := csv.NewReader(f)
// Read in all of the CSV records
reader.FieldsPerRecord = 4
testData, err := reader.ReadAll()
if err != nil {
    log.Fatal(err)
}
// Loop over the holdout data predicting y and evaluating the prediction
// with the mean absolute error.
var mAE float64
for i, record := range testData {
    // Skip the header.
    if i == 0 {
        continue
    }
    // Parse the Sales.
    yObserved, err := strconv.ParseFloat(record[3], 64)
    if err != nil {
        log.Fatal(err)
    }
    // Parse the TV value.
    tvVal, err := strconv.ParseFloat(record[0], 64)
    if err != nil {
        log.Fatal(err)
    }
    // Parse the Radio value.
    radioVal, err := strconv.ParseFloat(record[1], 64)
    if err != nil {
        log.Fatal(err)
    }
    // Parse the Newspaper value.
    newspaperVal, err := strconv.ParseFloat(record[2], 64)
    if err != nil {
        log.Fatal(err)
    }
    // Predict y with our trained model.
    yPredicted := predict(tvVal, radioVal, newspaperVal)
    // Add the to the mean absolute error.
    mAE += math.Abs(yObserved-yPredicted) / float64(len(testData))
}
// Output the MAE to standard out.
fmt.Printf("\nMAE = %0.2f\n\n", mAE)

編譯並運行,等到新的MAE

$ go build
$ ./myprogram

MAE = 1.26

注意在模型中添加Newspaper並不會提高MAE,因此在這種場景下並不適合添加Newspaper項,因為此時提高了複雜度,但並沒有顯著影響到模型。

添加到模型中的任何複雜性或複雜度都應該有其可衡量的理由。使用一個複雜的模型,這看起來很有趣,但同時也會讓人頭疼。

總結