數學推導+Python實現機器學習演算法:線性迴歸

步入量化學習艾莉絲發表於2018-09-27

很多同學在學習機器學習的時候,理論粗略看一遍之後就直接上手程式設計了,非常值得表揚。但是他不是真正的上手寫演算法,而是去直接呼叫 sklearn 這樣的 package,這就不大妥當了。筆者不是說調包不好,在實際工作和研究中,封裝好的簡單易用的 package 給我們的工作帶來了莫大的便利,大大提高了我們機器學習模型和演算法的實現效率。但這僅限於使用過程中。


筆者相信很多有企圖心的同學肯定不滿足於僅僅去使用這些 package 而不知模型和演算法的細節。所以,如果你是一名機器學習演算法的學習者,在學習過程中最好不要一上來就使用這些封裝好的包,而是根據自己對演算法的理解,在手推過模型和演算法的數學公式後,僅依靠 numpy 和 pandas 等基礎包的情況下手寫機器學習演算法。如此一遍過程之後,再去學習如何呼叫 sklearn 等機器學習庫,相信各位更能體會到調包的便利和樂趣。之後再去找資料實戰和打比賽做專案,相信你一定會成為一名優秀的機器學習演算法工程師。


本機器學習系列文章的兩個主題就是數學推導+純 numpy 實現。第一講我們從最基礎的線性迴歸模型開始。相信大家對迴歸演算法一定是相當熟悉了,特別是我們們有統計背景的同學。所以,筆者直接上數學推導。


迴歸分析的數學推導

本來想著上筆者的手推草稿的,但字跡過於張揚,在 word 裡或者用 markdown 寫公式又太耗費時間,這裡就直接借用周志華老師的機器學習教材上的推導過程:

數學推導+Python實現機器學習演算法:線性迴歸
數學推導+Python實現機器學習演算法:線性迴歸
數學推導+Python實現機器學習演算法:線性迴歸
數學推導+Python實現機器學習演算法:線性迴歸
數學推導+Python實現機器學習演算法:線性迴歸


推廣到矩陣形式就是:

數學推導+Python實現機器學習演算法:線性迴歸
數學推導+Python實現機器學習演算法:線性迴歸


以上便是線性迴歸模型中引數估計的推導過程。


迴歸分析的 numpy 實現

按照慣例,動手寫演算法之前我們需要理一下編寫思路。迴歸模型主體部分較為簡單,關鍵在於如何在給出 mse 損失函式之後基於梯度下降的引數更新過程。首先我們需要寫出模型的主體和損失函式以及基於損失函式的引數求導結果,然後對引數進行初始化,最後寫出基於梯度下降法的引數更新過程。當然,我們也可以寫出交叉驗證來得到更加穩健的引數估計值。話不多說,直接上程式碼。


迴歸模型主體:

import numpy as np

def linear_loss(X, y, w, b):
    num_train = X.shape[0]
    num_feature = X.shape[1]    
    # 模型公式
    y_hat = np.dot(X, w) + b    
    # 損失函式
    loss = np.sum((y_hat-y)**2)/num_train    
    # 引數的偏導
    dw = np.dot(X.T, (y_hat-y)) /num_train
    db = np.sum((y_hat-y)) /num_train    
    return y_hat, loss, dw, db
複製程式碼

引數初始化:

def initialize_params(dims):
    w = np.zeros((dims, 1))
    b = 0
    return w, b
複製程式碼

基於梯度下降的模型訓練過程:

def linar_train(X, y, learning_rate, epochs):
    w, b = initialize(X.shape[1])  
    loss_list = []  
    for i in range(1, epochs):        
    # 計算當前預測值、損失和引數偏導
        y_hat, loss, dw, db = linar_loss(X, y, w, b)  
        loss_list.append(loss)      
        # 基於梯度下降的引數更新過程
        w += -learning_rate * dw
        b += -learning_rate * db        
        # 列印迭代次數和損失

        if i % 10000 == 0:
            print('epoch %d loss %f' % (i, loss)) 
               
        # 儲存引數
        params = {            
            'w': w,            
            'b': b
        }        
        
        # 儲存梯度
        grads = {            
            'dw': dw,            
            'db': db
        }    
            
    return loss_list, loss, params, grads
複製程式碼

以上便是線性迴歸模型的基本實現過程。下面以 sklearn 中的 diabetes 資料集為例進行簡單的訓練。


資料準備:

from sklearn.datasets import load_diabetes
from sklearn.utils import shuffle

diabetes = load_diabetes()
data = diabetes.data
target = diabetes.target 

# 打亂資料
X, y = shuffle(data, target, random_state=13)
X = X.astype(np.float32)

# 訓練集與測試集的簡單劃分
offset = int(X.shape[0] * 0.9)
X_train, y_train = X[:offset], y[:offset]
X_test, y_test = X[offset:], y[offset:]
y_train = y_train.reshape((-1,1))
y_test = y_test.reshape((-1,1))

print('X_train=', X_train.shape)
print('X_test=', X_test.shape)
print('y_train=', y_train.shape)
print('y_test=', y_test.shape)
複製程式碼

數學推導+Python實現機器學習演算法:線性迴歸


執行訓練:

loss_list, loss, params, grads = linar_train(X_train, y_train, 0.001, 100000)
複製程式碼

數學推導+Python實現機器學習演算法:線性迴歸


檢視訓練得到的迴歸模型引數值:

print(params)
複製程式碼

數學推導+Python實現機器學習演算法:線性迴歸


下面定義一個預測函式對測試集結果進行預測:

def predict(X, params):
    w = params['w']
    b = params['b']

    y_pred = np.dot(X, w) + b    
    return y_pred

y_pred = predict(X_test, params)
y_pred[:5]
複製程式碼

數學推導+Python實現機器學習演算法:線性迴歸


利用 matplotlib 對預測結果和真值進行展示:

import matplotlib.pyplot as plt
f = X_test.dot(params['w']) + params['b']

plt.scatter(range(X_test.shape[0]), y_test)
plt.plot(f, color = 'darkorange')
plt.xlabel('X')
plt.ylabel('y')
plt.show()
複製程式碼

數學推導+Python實現機器學習演算法:線性迴歸

可見全變數的資料對於線性迴歸模型的擬合併不好,一來資料本身的分佈問題,二來簡單的線性模型對於該資料擬合效果差。當然,我們只是為了演示線性迴歸模型的基本過程,不要在意效果。


訓練過程中損失的下降:

plt.plot(loss_list, color = 'blue')
plt.xlabel('epochs')
plt.ylabel('loss')
plt.show()
複製程式碼

數學推導+Python實現機器學習演算法:線性迴歸


封裝一個線性迴歸類

筆者對上述過程進行一個簡單的 class 封裝,其中加入了自定義的交叉驗證過程進行訓練:

import numpy as np
from sklearn.utils import shuffle
from sklearn.datasets import load_diabetes

class lr_model():    
    def __init__(self):        
        pass

    def prepare_data(self):
        data = load_diabetes().data
        target = load_diabetes().target
        X, y = shuffle(data, target, random_state=42)
        X = X.astype(np.float32)
        y = y.reshape((-1, 1))
        data = np.concatenate((X, y), axis=1)        
        return data   
         
    def initialize_params(self, dims):
        w = np.zeros((dims, 1))
        b = 0
        return w, b    
    
    def linear_loss(self, X, y, w, b):
        num_train = X.shape[0]
        num_feature = X.shape[1]

        y_hat = np.dot(X, w) + b
        loss = np.sum((y_hat-y)**2) / num_train
        dw = np.dot(X.T, (y_hat - y)) / num_train
        db = np.sum((y_hat - y)) / num_train        
        return y_hat, loss, dw, db    
    
    def linear_train(self, X, y, learning_rate, epochs):
        w, b = self.initialize_params(X.shape[1])        
        for i in range(1, epochs):
            y_hat, loss, dw, db = self.linear_loss(X, y, w, b)
            w += -learning_rate * dw
            b += -learning_rate * db            
            if i % 10000 == 0:
                print('epoch %d loss %f' % (i, loss))
            
            params = {                
                'w': w,                
                'b': b
            }
            grads = {                
                'dw': dw,                
                'db': db
            }        
         return loss, params, grads    
    
     def predict(self, X, params):
        w = params['w']
        b = params['b']
        y_pred = np.dot(X, w) + b        
        return y_pred    
        
     def linear_cross_validation(self, data, k, randomize=True):        
        if randomize:
            data = list(data)
            shuffle(data)

        slices = [data[i::k] for i in range(k)]        
        for i in range(k):
            validation = slices[i]
            train = [data                        
            for s in slices if s is not validation for data in s]
            train = np.array(train)
            validation = np.array(validation)            
            yield train, validation
            
            
if __name__ == '__main__':
    lr = lr_model()
    data = lr.prepare_data()   
 
    for train, validation in lr.linear_cross_validation(data, 5):
        X_train = train[:, :10]
        y_train = train[:, -1].reshape((-1, 1))
        X_valid = validation[:, :10]
        y_valid = validation[:, -1].reshape((-1, 1))

        loss5 = []
        loss, params, grads = lr.linear_train(X_train, y_train, 0.001, 100000)
        loss5.append(loss)
        score = np.mean(loss5)
        print('five kold cross validation score is', score)
        y_pred = lr.predict(X_valid, params)
        valid_score = np.sum(((y_pred - y_valid) ** 2)) / len(X_valid)
        print('valid score is', valid_score)
複製程式碼

以上便是本節的內容,基於 numpy 手動實現一個簡單的線性迴歸模型。


參考資料:

周志華 機器學習

相關文章