機器學習練習(一)——簡單線性迴歸
這篇文章是一系列 Andrew Ng 在 Coursera 上的機器學習課程的練習的一部分。這篇文章的原始程式碼,練習文字,資料檔案可從這裡獲得。
這些年來,我專業開發的一個關鍵時刻是當我發現 Coursera 的時候。我以前聽過慕課現象,但我從來沒有時間深入進去並上一門課。今年早些時候我終於“鋌而走險”(pull the trigger)並註冊了一門 Andrew Ng 的機器學習課程。我完成了從開始到結束所有的步驟,包括所有的程式設計練習。這次經歷為我打開了新世界的大門,讓我感到了這個教育平臺的力量,我從此沉迷其中,不可自拔。
這篇文章將是一系列包含程式設計練習的 Andrew 課程的第一篇。
這個課程我不太在乎的一方面是用 Octave 做作業。
儘管 Octave/Matlab 是一個好的平臺,但大多數現實世界的“資料科學家”不是用 R 就是用 Python(當然,也是其他的語言和工具正在被使用,但這兩個毫無疑問排在首位)。自從我打算開發我的 Python 技能,我就決定通過這些練習從零開始學習 Python 。完整的原始碼可以在我的 GitHub 上的
如果你有興趣,你也可以在根目錄的子資料夾下找到這些練習中用到的資料和原始練習 PDF 。
雖然隨著時間的推移我可以解釋一些練習中涉及的概念,但我不可能解釋所有的你可能需要完全理解的資訊。如果你真的被機器學習所吸引但還沒有深入瞭解,我推薦你去看這門課(它是完全免費的,也不需要做任何承諾)。正因為如此,讓我們開始吧!
檢查資料(Examining The Data)
在 練習 1 的第一部分,我們被要求用簡單線性迴歸實現預測食物卡車的利潤。假設你是一家連鎖餐館的 CEO 正在考慮在不同的城市開一家新店。連鎖店在不同的城市已經有卡車了並且你有不同城市人口和利潤的資料。你想找出一個新的食品卡車期望的利潤,通過只給出它將被放置的城市的人口數。
讓我們通過檢查資料開始,資料在我的 repository 上的 data 目錄下的 ex1data1.txt 檔案中。首先,我需要匯入一些庫。
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
現在,讓我們開始吧。我們可以使用 pandas 將資料載入到資料框(DataFrame)並且用 head 函式顯示前幾行。
path = os.getcwd() + '\data\ex1data1.txt'
data = pd.read_csv(path, header=None, names=['Population', 'Profit'])
data.head()
Population | Profit | |
---|---|---|
0 | 6.1101 | 17.5920 |
1 | 5.5277 | 9.1302 |
2 | 8.5186 | 13.6620 |
3 | 7.0032 | 11.8540 |
4 | 5.8598 | 6.8233 |
pandas 提供的另一個開箱即用的函式是 describe,它能夠對一個數據集進行基本的統計計算。這對於在一個專案的探索分析階段對資料的整體“感覺”有幫助。
data.describe()
Population | Profit | |
---|---|---|
count | 97.000000 | 97.000000 |
mean | 8.159800 | 5.839135 |
std | 3.869884 | 5.510262 |
min | 5.026900 | -2.680700 |
25% | 5.707700 | 1.986900 |
50% | 6.589400 | 4.562300 |
75% | 8.578100 | 7.046700 |
max | 22.203000 | 24.147000 |
檢查你資料的統計是有幫助的,但有時候你需要讓它視覺化。幸運的是這個資料集只有一個因變數,所以我們可以把它放在散點圖( scatter plot)來看看它到底是什麼樣子的。我們可以使用 pandas 提供的 plot 函式來完成,它真的只是 matplotlib 的包裝。
data.plot(kind='scatter', x='Population', y='Profit', figsize=(12,8))
它真的有助於實際上看看到底是什麼樣的,不是嗎?我們可以清楚的看到較少人口城市值的叢集,還有某種程度上增長的利潤和城市規模增長的線性趨勢。現在讓我們來到有趣的部分——從零開始用 python 實現一種線性迴歸演算法!
實現簡單線性迴歸
如果你不熟悉線性迴歸,它是一種建立一個因變數和一個或多個自變數(如果只有一個自變數那麼它被稱為簡單線性迴歸,如果有多個自變數那麼他被稱為多元線性迴歸)之間關係的方法。有很多不同型別的線性迴歸並且差異很大,這超過這裡討論的範圍,所以我不會去管,簡單地說——我們試圖建立一個關於資料
在這個實現中,我們將要使用一種被稱為梯度下降(gradient descent)的技術來找出引數
已經有足夠的理論了。讓我們來寫一些程式碼吧。首先我們需要一個代價(cost)函式。代價函式評估我們模型的質量,通過(模型引數和實際資料點)計算我們模型對資料點的預測的誤差(error)。例如,如果給定城市人口是 4 而我們預測值是 7,我們的誤差就是 (7-4)^2=3^2=9(假設一個L2約束或“最小二乘”損失函式)。我們對每個
def computeCost(X, y, theta):
inner = np.power(((X * theta.T) - y), 2)
return np.sum(inner) / (2 * len(X))
注意這裡沒有迴圈。我們利用 numpy 的線性代數能力計算一系列矩陣運算的結果。這遠比未優化的 for 迴圈計算效率高。
為了讓代價函式能夠無縫地應用在我們之前建立的 pandas 資料框上,我們需要做一些操作。首先,為了使矩陣運算正確(我不會詳細的說明為什麼這是必要的,但如果你感興趣它在練習的文字中——基本上,它解釋為線上性方程中的截距項( intercept term))我們需要在資料框的開頭插入一列 1 。其次,我們需要將我們的資料分為自變數
# append a ones column to the front of the data set
data.insert(0, 'Ones', 1)
# set X (training data) and y (target variable)
cols = data.shape[1]
X = data.iloc[:,0:cols-1]
y = data.iloc[:,cols-1:cols]
最後,我們需要將我們的資料框轉化為 numpy 矩陣並且例項化引數矩陣。
# convert from data frames to numpy matrices
X = np.matrix(X.values)
y = np.matrix(y.values)
theta = np.matrix(np.array([0,0]))
記住除錯矩陣運算的一個有用的技巧是看你處理矩陣的 shape(譯者注:形狀,行列數)。它也有助於記住當在你腦海中進行矩陣乘法看起來像 (i x j) * (j x k) = (i x k)時,i,j,k 是矩陣相應維的 shape。
X.shape, theta.shape, y.shape
((97L, 2L), (1L, 2L), (97L, 1L))
現在,我們測試我們的代價函數了。記住引數被初始化為
computeCost(X, y, theta)
32.072733877455676
到目前為止執行良好。現在我們需要定義一個函式用定義在練習文字中的更新規則在引數
def gradientDescent(X, y, theta, alpha, iters):
temp = np.matrix(np.zeros(theta.shape))
parameters = int(theta.ravel().shape[1])
cost = np.zeros(iters)
for i in range(iters):
error = (X * theta.T) - y
for j in range(parameters):
term = np.multiply(error, X[:,j])
temp[0,j] = theta[0,j] - ((alpha / len(X)) * np.sum(term))
theta = temp
cost[i] = computeCost(X, y, theta)
return theta, cost
梯度下降的思想是對於每次迭代,為了找出合適的方向來移動我們的引數向量我們計算出誤差項(error term )的梯度。換句話說,我們為了減少誤差正在計算引數需要的改變,從而使我們的結果更接近最佳的結果(即最合適)。
這是一個相當複雜的話題,我可以用一個整篇部落格來討論梯度下降。如果你對學習更多這方面的東西感興趣的話,我建議從這篇文章開始,並從那裡拓展出來。
我們又一次依靠 numpy 和線性代數獲得解。你可能注意到我的實現並不是 100% 最優的。特別的,有方法能消除內迴圈並一次性更新所有引數。暫時,我將把它留給讀者(我將在下一篇文章中討論它)。
現在,我們已經有一種方法評估結果也有了一種方法找出好的結果,是時候將它們應用到我們的資料集了。
# initialize variables for learning rate and iterations
alpha = 0.01
iters = 1000
# perform gradient descent to "fit" the model parameters
g, cost = gradientDescent(X, y, theta, alpha, iters)
g
matrix([[-3.24140214, 1.1272942 ]])
注意,在這裡我們已經初始化了一些新變數。如果你仔細看看梯度下降函式,它有引數叫alpha(
我們現在有一個引數向量描述對我們的資料集而言什麼是我們認為的最優線性模型。一個快速的方法來評估我們的迴歸模型有多好是——看看我們的新解對資料集的總誤差是多少:
computeCost(X, y, g)
4.5159555030789118
這確實比 32 好,但它並不是一個直觀的方法。幸運的是,有一些其他的技術任我們使用。
觀察結果
我們現在將使用 matplotlib 來使我們的結果視覺化。還記得之前的散點圖嗎?讓我們放一條代表我們模型的線在散點圖上來看看它有多適合。我們可以用 numpy 的 linspace 函式來建立均勻分佈的點在我們資料範圍內,然後用我們的模型評估這些點,來看看期望的利潤是多少。然後,我們可以把它做成一個線圖並繪製出來。
x = np.linspace(data.Population.min(), data.Population.max(), 100)
f = g[0, 0] + (g[0, 1] * x)
fig, ax = plt.subplots(figsize=(12,8))
ax.plot(x, f, 'r', label='Prediction')
ax.scatter(data.Population, data.Profit, label='Traning Data')
ax.legend(loc=2)
ax.set_xlabel('Population')
ax.set_ylabel('Profit')
ax.set_title('Predicted Profit vs. Population Size')
不錯!我們的解對於這個資料集看起來像是最優線性模型。由於梯度下降函式也輸出了一個向量和每次訓練迭代的代價(cost),我們也可以繪製出來。
fig, ax = plt.subplots(figsize=(12,8))
ax.plot(np.arange(iters), cost, 'r')
ax.set_xlabel('Iterations')
ax.set_ylabel('Cost')
ax.set_title('Error vs. Training Epoch')
注意,代價總是降低——這是一個被稱為凸優化問題(convex optimization problem)的一個例子。如果你對這個問題繪製了整個解空間(即對所有每可能的引數值繪製將代價當作一個模型引數的函式)你會看到它像一個“盆地”的“碗”形它代表最優解。
暫時就這些!在 第二部分 我們將通過拓展這個例子到多個變數來結束第一個練習。我也將展示上面的解如何用一個受歡迎的機器學習庫 scikit-learn 來實現。