對(duì)于大多數(shù)數(shù)據(jù)科學(xué)家而言，線性回歸方法是他們進(jìn)行統(tǒng)計(jì)學(xué)建模和預(yù)測(cè)分析任務(wù)的起點(diǎn)。這種方法已經(jīng)存在了 200 多年，并得到了廣泛研究，但仍然是一個(gè)積極的研究領(lǐng)域。由于良好的可解釋性，線性回歸在商業(yè)數(shù)據(jù)上的用途十分廣泛。當(dāng)然，在生物數(shù)據(jù)、工業(yè)數(shù)據(jù)等領(lǐng)域也不乏關(guān)于回歸分析的應(yīng)用。

另一方面，Python 已成為數(shù)據(jù)科學(xué)家首選的編程語(yǔ)言，能夠應(yīng)用多種方法利用線性模型擬合大型數(shù)據(jù)集顯得尤為重要。

如果你剛剛邁入機(jī)器學(xué)習(xí)的大門(mén)，那么使用 Python 從零開(kāi)始對(duì)整個(gè)線性回歸算法進(jìn)行編碼是一次很有意義的嘗試，讓我們來(lái)看看怎么做吧。

數(shù)據(jù)

機(jī)器學(xué)習(xí)問(wèn)題的第一步是獲取數(shù)據(jù)，沒(méi)有可以學(xué)習(xí)的數(shù)據(jù)就沒(méi)有機(jī)器學(xué)習(xí)。本文將使用非常常規(guī)的線性回歸數(shù)據(jù)集——房?jī)r(jià)預(yù)測(cè)數(shù)據(jù)集。

這是一個(gè)包含俄勒岡州波特蘭市房?jī)r(jià)的簡(jiǎn)單數(shù)據(jù)集。該數(shù)據(jù)集中第一列是房屋面積（以平方英尺為單位），第二列是臥室的數(shù)量，第三列是房屋價(jià)格。該數(shù)據(jù)集中有多個(gè)特征（例如，house_size 和房間數(shù)），因此我們將研究多元線性回歸，標(biāo)簽 (y) 是我們將要預(yù)測(cè)的房?jī)r(jià)。

首先定義用于加載數(shù)據(jù)集的函數(shù)：

def load_data(filename):
  df = pd.read_csv(filename, sep=",", index_col=False)
  df.columns = ["housesize", "rooms", "price"]
  data = np.array(df, dtype=float)
  plot_data(data[:,:2], data[:, -1])
  normalize(data)
    return data[:,:2], data[:, -1]

我們稍后將調(diào)用上述函數(shù)來(lái)加載數(shù)據(jù)集。此函數(shù)返回 x 和 y。

歸一化數(shù)據(jù)

上述代碼不僅加載數(shù)據(jù)，還對(duì)數(shù)據(jù)執(zhí)行歸一化處理并繪制數(shù)據(jù)點(diǎn)。在查看數(shù)據(jù)圖之前，我們首先了解上述代碼中的 normalize(data)。

查看原始數(shù)據(jù)集后，你會(huì)發(fā)現(xiàn)第二列數(shù)據(jù)的值（房間數(shù)量）比第一列（即房屋面積）小得多。該模型不會(huì)將此數(shù)據(jù)評(píng)估為房間數(shù)量或房屋面積，對(duì)于模型來(lái)說(shuō)，它們只是一些數(shù)字。機(jī)器學(xué)習(xí)模型中某些列（或特征）的數(shù)值比其他列高可能會(huì)造成不想要的偏差，還可能導(dǎo)致方差和數(shù)學(xué)均值的不平衡。出于這些原因，也為了簡(jiǎn)化工作，我們建議先對(duì)特征進(jìn)行縮放或歸一化，使其位于同一范圍內(nèi)（例如 [-1,1] 或 [0,1]），這會(huì)讓訓(xùn)練容易許多。因此我們將使用特征歸一化，其數(shù)學(xué)表達(dá)如下：

• Z = (x — μ) / σ
• μ : mean
• σ : standard deviation

其中 z 是歸一化特征，x 是非歸一化特征。有了歸一化公式，我們就可以為歸一化創(chuàng)建一個(gè)函數(shù)：

def normalize(data):
  for i in range(0,data.shape[1]-1):
        data[:,i] = ((data[:,i] - np.mean(data[:,i]))/np.std(data[:, i]))

上述代碼遍歷每一列，并使用每一列中所有數(shù)據(jù)元素的均值和標(biāo)準(zhǔn)差對(duì)其執(zhí)行歸一化。

繪制數(shù)據(jù)

在對(duì)線性回歸模型進(jìn)行編碼之前，我們需要先問(wèn)「為什么」。

為什么要使用線性回歸解決這個(gè)問(wèn)題？這是一個(gè)非常有用的問(wèn)題，在寫(xiě)任何具體代碼之前，你都應(yīng)該非常清楚要使用哪種算法，以及在給定數(shù)據(jù)集和待解決問(wèn)題的情況下，這是否真的是最佳選擇。

我們可以通過(guò)繪制圖像來(lái)證明對(duì)當(dāng)前數(shù)據(jù)集使用線性回歸有效的原因。為此，我們?cè)谏厦娴?load_data 中調(diào)用了 plot_data 函數(shù)，現(xiàn)在我們來(lái)定義一下 plot_data 函數(shù)：

def plot_data(x, y):
  plt.xlabel('house size')
  plt.ylabel('price')
  plt.plot(x[:,0], y, 'bo')
    plt.show()

調(diào)用該函數(shù)，將生成下圖：

房屋面積與房屋價(jià)格關(guān)系圖。

如上圖所示，我們可以粗略地?cái)M合一條線。這意味著使用線性近似能夠做出較為準(zhǔn)確的預(yù)測(cè)，因此可以采用線性回歸。

準(zhǔn)備好數(shù)據(jù)之后就要進(jìn)行下一步，給算法編寫(xiě)代碼。

假設(shè)

首先我們需要定義假設(shè)函數(shù)，稍后我們將使用它來(lái)計(jì)算代價(jià)。對(duì)于線性回歸，假設(shè)是：

但數(shù)據(jù)集中只有 2 個(gè)特征，因此對(duì)于當(dāng)前問(wèn)題，假設(shè)是：

其中 x1 和 x2 是兩個(gè)特征（即房屋面積和房間數(shù)量）。然后編寫(xiě)一個(gè)返回該假設(shè)的簡(jiǎn)單 Python 函數(shù)：

def h(x,theta):
    return np.matmul(x, theta)

接下來(lái)我們來(lái)看代價(jià)函數(shù)。

代價(jià)函數(shù)

使用代價(jià)函數(shù)的目的是評(píng)估模型質(zhì)量。

代價(jià)函數(shù)的等式為：

代價(jià)函數(shù)的代碼如下：

def cost_function(x, y, theta):
    return ((h(x, theta)-y).T@(h(x, theta)-y))/(2*y.shape[0])

到目前為止，我們定義的所有 Python 函數(shù)都與上述線性回歸的數(shù)學(xué)意義完全相同。接下來(lái)我們需要將代價(jià)最小化，這就要用到梯度下降。

梯度下降

梯度下降是一種優(yōu)化算法，旨在調(diào)整參數(shù)以最小化代價(jià)函數(shù)。

梯度下降的主要更新步是：

因此，我們將代價(jià)函數(shù)的導(dǎo)數(shù)乘以學(xué)習(xí)率（α），然后用參數(shù)（θ）的當(dāng)前值減去它，獲得新的更新參數(shù)（θ）。

def gradient_descent(x, y, theta, learning_rate=0.1, num_epochs=10):
  m = x.shape[0]
  J_all = []
  
  for _ in range(num_epochs):
    h_x = h(x, theta)
    cost_ = (1/m)*(x.T@(h_x - y))
    theta = theta - (learning_rate)*cost_
    J_all.append(cost_function(x, y, theta))

    return theta, J_all

gradient_descent 函數(shù)返回 theta 和 J_all。theta 顯然是參數(shù)向量，其中包含假設(shè)的θs 值，J_all 是一個(gè)列表，包含每個(gè) epoch 后的代價(jià)函數(shù)。J_all 變量并非必不可少，但它有助于更好地分析模型。

整合到一起

接下來(lái)要做的就是以正確的順序調(diào)用函數(shù)

x,y = load_data("house_price_data.txt")
y = np.reshape(y, (46,1))
x = np.hstack((np.ones((x.shape[0],1)), x))
theta = np.zeros((x.shape[1], 1))
learning_rate = 0.1
num_epochs = 50
theta, J_all = gradient_descent(x, y, theta, learning_rate, num_epochs)
J = cost_function(x, y, theta)
print("Cost: ", J)
print("Parameters: ", theta)

#for testing and plotting cost 
n_epochs = []
jplot = []
count = 0
for i in J_all:
  jplot.append(i[0][0])
  n_epochs.append(count)
  count += 1
jplot = np.array(jplot)
n_epochs = np.array(n_epochs)
plot_cost(jplot, n_epochs)

test(theta, [1600, 2])

首先調(diào)用 load_data 函數(shù)載入 x 和 y 值。x 值包含訓(xùn)練樣本，y 值包含標(biāo)簽（在這里就是房屋的價(jià)格）。

你肯定注意到了，在整個(gè)代碼中，我們一直使用矩陣乘法的方式來(lái)表達(dá)所需。例如為了得到假設(shè)，我們必須將每個(gè)參數(shù)（θ）與每個(gè)特征向量（x）相乘。我們可以使用 for 循環(huán)，遍歷每個(gè)樣本，每次都執(zhí)行一次乘法，但如果訓(xùn)練的樣本過(guò)多，這可能不是最高效的方法。

在這里更有效的方式是使用矩陣乘法。本文所用的數(shù)據(jù)集具備兩個(gè)特征：房屋面積和房間數(shù)，即我們有（2+1）三個(gè)參數(shù)。將假設(shè)看作圖形意義上的一條線，用這種方式來(lái)思考額外參數(shù)θ0，最終額外的θ0 也要使這條線符合要求。

有利的假設(shè)函數(shù)圖示。

現(xiàn)在我們有了三個(gè)參數(shù)和兩個(gè)特征。這意味著θ或參數(shù)向量（1 維矩陣）的維數(shù)是 (3,1)，但特征向量的維度是 (46,2)。你肯定會(huì)注意到將這樣兩個(gè)矩陣相乘在數(shù)學(xué)上是不可能的。再看一遍我們的假設(shè)：

如果你仔細(xì)觀察的話，實(shí)際上這很直觀：如果在特征向量 (x) {維度為 (46, 3)} 的開(kāi)頭添加額外的一列，并且對(duì) x 和 theta 執(zhí)行矩陣乘法，將得出 hθ(x) 的方程。

記住，在實(shí)際運(yùn)行代碼來(lái)實(shí)現(xiàn)此功能時(shí)，不會(huì)像 hθ(x) 那樣返回表達(dá)式，而是返回該表達(dá)式求得的數(shù)學(xué)值。在上面的代碼中，x = np.hstack((np.ones((x.shape[0],1)), x)) 這一行在 x 開(kāi)頭加入了額外一列，以備矩陣乘法需要。

在這之后，我們用零初始化 theta 向量，當(dāng)然你也可以用一些小隨機(jī)值來(lái)進(jìn)行初始化。我們還指定了訓(xùn)練學(xué)習(xí)率和 epoch 數(shù)。

定義完所有超參數(shù)之后，我們就可以調(diào)用梯度下降函數(shù)，以返回所有代價(jià)函數(shù)的歷史記錄以及參數(shù) theta 的最終向量。在這里 theta 向量定義了最終的假設(shè)。你可能注意到，由梯度下降函數(shù)返回的 theta 向量的維度為 (3,1)。

還記得函數(shù)的假設(shè)嗎？

所以我們需要三個(gè)θ，theta 向量的維度為 (3,1)，因此 theta [0]、theta [1] 和 theta [2] 實(shí)際上分別為θ0、θ1 和 θ2。J_all 變量是所有代價(jià)函數(shù)的歷史記錄。你可以打印出 J_all 數(shù)組，來(lái)查看代價(jià)函數(shù)在梯度下降的每個(gè) epoch 中逐漸減小的過(guò)程。

代價(jià)和 epoch 數(shù)量的關(guān)系圖。

我們可以通過(guò)定義和調(diào)用 plot_cost 函數(shù)來(lái)繪制此圖，如下所示：

def plot_cost(J_all, num_epochs):
  plt.xlabel('Epochs')
  plt.ylabel('Cost')
  plt.plot(num_epochs, J_all, 'm', linewidth = "5")
    plt.show()

現(xiàn)在我們可以使用這些參數(shù)來(lái)找到標(biāo)簽，例如給定房屋面積和房間數(shù)量時(shí)的房屋價(jià)格。

測(cè)試

現(xiàn)在你可以測(cè)試調(diào)用測(cè)試函數(shù)的代碼，該函數(shù)會(huì)將房屋面積、房間數(shù)量和 logistic 回歸模型返回的最終 theta 向量作為輸入，并輸出房屋價(jià)格。

def test(theta, x):
  x[0] = (x[0] - mu[0])/std[0]
  x[1] = (x[1] - mu[1])/std[1]

  y = theta[0] + theta[1]*x[0] + theta[2]*x[1]
    print("Price of house: ", y)

完整代碼

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

#variables to store mean and standard deviation for each feature
mu = []
std = []

def load_data(filename):
  df = pd.read_csv(filename, sep=",", index_col=False)
  df.columns = ["housesize", "rooms", "price"]
  data = np.array(df, dtype=float)
  plot_data(data[:,:2], data[:, -1])
  normalize(data)
  return data[:,:2], data[:, -1]

def plot_data(x, y):
  plt.xlabel('house size')
  plt.ylabel('price')
  plt.plot(x[:,0], y, 'bo')
  plt.show()

def normalize(data):
  for i in range(0,data.shape[1]-1):
    data[:,i] = ((data[:,i] - np.mean(data[:,i]))/np.std(data[:, i]))
    mu.append(np.mean(data[:,i]))
    std.append(np.std(data[:, i]))


def h(x,theta):
  return np.matmul(x, theta)

def cost_function(x, y, theta):
  return ((h(x, theta)-y).T@(h(x, theta)-y))/(2*y.shape[0])

def gradient_descent(x, y, theta, learning_rate=0.1, num_epochs=10):
  m = x.shape[0]
  J_all = []
  
  for _ in range(num_epochs):
    h_x = h(x, theta)
    cost_ = (1/m)*(x.T@(h_x - y))
    theta = theta - (learning_rate)*cost_
    J_all.append(cost_function(x, y, theta))

  return theta, J_all 

def plot_cost(J_all, num_epochs):
  plt.xlabel('Epochs')
  plt.ylabel('Cost')
  plt.plot(num_epochs, J_all, 'm', linewidth = "5")
  plt.show()

def test(theta, x):
  x[0] = (x[0] - mu[0])/std[0]
  x[1] = (x[1] - mu[1])/std[1]

  y = theta[0] + theta[1]*x[0] + theta[2]*x[1]
  print("Price of house: ", y)

x,y = load_data("house_price_data.txt")
y = np.reshape(y, (46,1))
x = np.hstack((np.ones((x.shape[0],1)), x))
theta = np.zeros((x.shape[1], 1))
learning_rate = 0.1
num_epochs = 50
theta, J_all = gradient_descent(x, y, theta, learning_rate, num_epochs)
J = cost_function(x, y, theta)
print("Cost: ", J)
print("Parameters: ", theta)

#for testing and plotting cost 
n_epochs = []
jplot = []
count = 0
for i in J_all:
  jplot.append(i[0][0])
  n_epochs.append(count)
  count += 1
jplot = np.array(jplot)
n_epochs = np.array(n_epochs)
plot_cost(jplot, n_epochs)

test(theta, [1600, 3])

總結(jié)

這就是線性回歸的全部代碼了。

現(xiàn)在你已經(jīng)學(xué)會(huì)了從零開(kāi)始成功編寫(xiě)線性回歸模型。能夠理解和編寫(xiě)整個(gè)算法并不是一件容易的事，你或許需要時(shí)不時(shí)地回看才能完全理解。但這些努力是值得的，線性回歸通常是人們學(xué)習(xí)機(jī)器學(xué)習(xí)算法的第一步，在這之后你可以選擇另一個(gè)適用于線性回歸處理的數(shù)據(jù)集，并嘗試剛寫(xiě)好的算法。

原文鏈接：

https://towardsdatascience.com/coding-linear-regression-from-scratch-c42ec079902

以上就是如何用Python徒手寫(xiě)線性回歸的詳細(xì)內(nèi)容，更多關(guān)于python 手寫(xiě)線性回歸的資料請(qǐng)關(guān)注腳本之家其它相關(guān)文章！