iT邦幫忙

1

【不用機器學習套件手刻機器學習code】Multiple Linear Regression

機器學習是當今程式熱門領域,
雖說當前已經有許多機器學習的套件可以使用,
然而實際將程式細節實作出來感覺比較能對機器學習算法的原理比較了解,
於是有點自虐的嘗試自己摸索程式細節

今天自學的主題是Multiple Linear Regression,
Linear Regression的核心概念是
Y的值跟m個特徵有關,
可以將Y的值用一個方程式來表達:
y = w1*f1(x)+ w2*f2(x)+...+ wM*fM(x)
其中x是一個向量,f是我們可以選擇的函數,
通常會選多項式函數比較簡單,
w1, w2, ..., wM是我們想要求解的參數

目標: 利用特徵x預測y值
應用: 例如運用股市資訊預測股價

關於linear regression的概念,
個人推薦可以看李宏毅老師的ML Lecture 1: Regression - Case Study
李老師講解的蠻生動清楚易懂的,
先去看老師的講解應該會比較有概念

直接用sklearn實作的版本

如果你想直接使用機器學習套件sklearn來做linear regression,
不想自己手刻程式,
可以參考sklearn.linear_model.LinearRegression
可以讓你快速實作出linear regression的程式

sklearn程式例子

這邊我們手動設計資料滿足關係式 y = 1 * x_0 + 2 * x_1 + 3
實際把資料餵給程式求解時,
會假設y = w0 * x_0 + w1 * x_1 + b
w0, w1, b 是我們要求解的參數

import numpy as np
from sklearn.linear_model import LinearRegression

"""
設計資料滿足關係式 y = 1 * x_0 + 2 * x_1 + 3,
實際把資料餵給程式求解時,是假設 y = w0 * x_0 + w1 * x_1 + b,
w0, w1, b 是我們要求解的參數
"""
X = np.array([[1, 1], 
              [1, 2], 
              [2, 2], 
              [2, 3]])
y = np.dot(X, np.array([1, 2])) + 3

# 用sklearn模組的函數直接求解參數
reg = LinearRegression().fit(X, y)
print(reg.coef_) # [1. 2.]
print(reg.intercept_) # 3.0000

# 給定x值,預測y值
print(reg.predict(np.array([[3, 5],
                            [6, 0]]))) #[16.  9.]

手刻程式碼之數學式子推導

以下講解手刻程式碼的一些細節:

假設Y=XB,Y是m*1 vector, X是m*n矩陣,B是n*1vector,
Y跟X都是已知,是從訓練資料來的,
B是我們要求解的參數

則 X_T*Y=X_T*X*B(用X_T表示X的轉置矩陣),
解得參數B= [(X_T*X)^-1]X_T*Y
若使用MAP方法公式解則為 B= [(X_TX)^-1+ lamb*I(單位矩陣)]X_TY
參考: Linear Regression

公式解B= [(X_TX)^-1+ lamb*I(單位矩陣)]X_TY

我們把公式解以程式的形式寫出

注意做Linear Regression時,
訓練資料的數量必須多於feature的數量,
否則參數公式計算的(X.T*X)會沒有反矩陣

def parameter(X, Y, lamb=0):
    return np.dot(np.linalg.inv(np.dot(X.T,X)+lamb*np.eye(X.shape[1])),np.dot(X.T,Y))

資料轉換

上述我們假設y與x 的關係可以寫成y = w1*x1+ w2*x2+...+ wM*xM
然而更general的情況可能會將x透過函數做轉換,
比如說: y = w1*f1(x)+ w2*f2(x)+...+ wM*fM(x)

舉例來說y的值跟x1, x2, x3有關,
關係式為y= w1*x1*x2+ w2*x2*x3 + w3*x3*x1+ b
x1, x2, x3, y是已知,w1, w2, w3, b是我們想求解的參數
這樣的關係式依然是一個線性的關係

至於函數的選取,簡單起見常會選多項式函數
這邊嘗試撰寫程式:

from itertools import combinations_with_replacement as CWR

"""
 計算一個陣列透過idx的指定乘積
 例:idx=(i1,i2,...in)
 return nums[i1]*nums[i2]*...*nums[in]
"""
def prodOfArr(nums,idxs):
    return reduce(lambda x,y: x*nums[y], idxs, 1)

"""
 生成x的多項式函數,x資料維度=x_dim,次方最高到polyDeg
 例x_dim=4, polyDeg=2
 則X的一筆data形式為 [x1,x2,x3,x4]
 會得到函數: [1,x1,x2,x3,x4, x1*x2, x1*x3,x1*x4,x2*x3,x2*x4,x3*x4]
"""
def designPolyFunc(x_dim, polyDeg=1):
    funcs = [lambda X:1]
    for i in range(1,polyDeg+1):
        funcs += [lambda X, idx=c: prodOfArr(X,idx) for c in CWR(range(x_dim),i)]
    return funcs


"""
 透過函數們funcs將原始資料X=
 x01 x02 ... x0n (= vector x0)
 x11 x12 ... x1n (= vector x1)
 ...
 xm1 xm2 ... xmn (= vector xm)
 其中(m是data數量,n是x的feature數量)
 轉換成矩陣形式:
 1 f1(x0) f2(x0) ... fk(x0)
 1 f1(x1) f2(x1) ... fk(x1)
 ...
 1 f1(xm) f2(xm) ... fk(xm)
"""
def preprocess(funcs, X):
    mat = [[f(x) for f in funcs] for x in X]
    return np.array(mat)

機器學習概念: training set與testing set

一般來說,在機器學習中我們不會直接把所有資料都拿去訓練,
而是會將資料分成訓練集(training data)和測試集(testing data),
就好比說讓學生看一系列的題庫(訓練集),
在給學生考試(測試集),檢測學生是否懂了

此部分並非Linear Regression的算法核心,
我想偷懶直接call sklearn的函數,
範例(x_data, y_data是原始資料):

from sklearn.model_selection import train_test_split
# from sklearn.cross_validation import train_test_split # for較舊的版本
x_train, x_test, y_train, y_test = train_test_split(x_data, y_data, test_size=0.2)

如果要自己實作的話,我是覺得pandas工具適合作資料處理:

def train_test_split(x_data, y_data, test_size=0.1):
    x_train = x_data.sample(frac=1-test_size)
    x_test = x_data.drop(x_train.index)
    
    y_train = y_data.iloc[x_train.index]
    y_test= y_data.iloc[x_test.index]

    # 將資料轉換成np可用的矩陣
    return (np.array(arr) for arr in [x_train, x_test, y_train, y_test])

建 Linear Regression Model做預測

有個基礎概念之後,我們要進入Linear Regression的算法核心,
這邊定義一個函數,
可以根據training data 去算去一組參數,
並用這組參數去對testing data 做預測,計算誤差

為什麼要計算誤差呢?
因為我們總會想知道做的預測準不準

# 計算預測的Y和真實Y之間的root mean square
def RMS(perdict_Y, real_Y):
    return math.sqrt(sum((real_Y-perdict_Y)**2)*2/len(real_Y))


# 把訓練資料、測試資料(預處理過的data)、lamb(若用MAP方法才要傳lamb),做一次訓練
# 回傳值: 訓練得到的最佳參數、訓練資料的RMS、測試資料的RMS
def training(x_train,x_test, y_train, y_test, lamb=0):
    beta = parameter(x_train,y_train, lamb)
    
    ## 以下用train_data做predict,計算train_data的RMS
    train_perdict_Y = np.dot(x_train,beta)
    train_error= RMS(train_perdict_Y,y_train)
    
    ## 以下用test_data做predict,計算test_data的RMS
    perdict_Y = np.dot(x_test,beta)
    test_error= RMS(perdict_Y,y_test)
    return beta, train_error, test_error

程式實作範例

有了前述基礎知識後,我們可以來嘗試手刻一個例子,
看看機器學習能否學出來,
這邊我們創建一個隨機的100筆資料,
每筆資料有3個特徵,
並且我們設計y = 0.5 * x_0 + 0.4 * x_1 + 0.7 * x_2 + 150
我們想看看用Linear Regression的方法能否把y與x的關係學出來

np.random.seed(931) # 這邊設固定的random seed,確保練習用的資料相同
X = np.random.randint(low=4,high=1000,size=(100,3))
Y = np.dot(X, np.array([0.5, 0.4, 0.7])) + 150
print(X)
print(Y)

生成資料結果:

[[390 927 691]
 [  6 670 949]
 [157 456 403]
 [692 677 328]
 ...]
 [1199.5 1085.3  693.   996.4 ...]

舉例來說,第一筆x資料是[390 927 691],
對應的y即為0.5 * 390 + 0.4 * 927 + 0.7 * 691 + 150 = 1199.5

完整程式示例

from functools import reduce
import numpy as np
import math
from sklearn.model_selection import train_test_split
from itertools import combinations_with_replacement as CWR

"""
 計算一個陣列透過idx的指定乘積
 例:idx=(i1,i2,...in)
 return nums[i1]*nums[i2]*...*nums[in]
"""
def prodOfArr(nums,idxs):
    return reduce(lambda x,y: x*nums[y], idxs, 1)

"""
 生成x的多項式函數,x資料維度=x_dim,次方最高到polyDeg
 例x_dim=4, polyDeg=2
 則X的一筆data形式為 [x1,x2,x3,x4]
 會得到函數: [1,x1,x2,x3,x4, x1*x2, x1*x3,x1*x4,x2*x3,x2*x4,x3*x4]
"""
def designPolyFunc(x_dim, polyDeg=1):
    funcs = [lambda X:1]
    for i in range(1,polyDeg+1):
        funcs += [lambda X, idx=c: prodOfArr(X,idx) for c in CWR(range(x_dim),i)]
    return funcs


"""
 透過函數們funcs將原始資料X=
 x01 x02 ... x0n (= vector x0)
 x11 x12 ... x1n (= vector x1)
 ...
 xm1 xm2 ... xmn (= vector xm)
 其中(m是data數量,n是x的feature數量)
 轉換成矩陣形式:
 1 f1(x0) f2(x0) ... fk(x0)
 1 f1(x1) f2(x1) ... fk(x1)
 ...
 1 f1(xm) f2(xm) ... fk(xm)
"""
def preprocess(funcs, X):
    mat = [[f(x) for f in funcs] for x in X]
    return np.array(mat)

# 公式解B= [(X_TX)^-1+ lamb\*I(單位矩陣)]X_TY
def parameter(X, Y, lamb=0):
    return np.dot(np.linalg.inv(np.dot(X.T,X)+lamb*np.eye(X.shape[1])),np.dot(X.T,Y))

# 計算預測的Y和真實Y之間的root mean square
def RMS(perdict_Y, real_Y):
    return math.sqrt(sum((real_Y-perdict_Y)**2)*2/len(real_Y))


# 把訓練資料、測試資料(預處理過的data)、lamb(若用MAP方法才要傳lamb),做一次訓練
# 回傳值: 訓練得到的最佳參數、訓練資料的RMS、測試資料的RMS
def training(x_train,x_test, y_train, y_test, lamb=0):
    beta = parameter(x_train,y_train, lamb)
    
    ## 以下用train_data做predict,計算train_data的RMS
    train_perdict_Y = np.dot(x_train,beta)
    train_error= RMS(train_perdict_Y,y_train)
    
    ## 以下用test_data做predict,計算test_data的RMS
    perdict_Y = np.dot(x_test,beta)
    test_error= RMS(perdict_Y,y_test)
    return beta, train_error, test_error



if __name__ == '__main__':
    np.random.seed(931) # 這邊設固定的random seed,確保練習用的資料相同
    X = np.random.randint(low=4,high=1000,size=(100,3))
    Y = np.dot(X, np.array([0.5, 0.4, 0.7])) + 150
    
    # 資料預處理
    funcs = designPolyFunc(3, polyDeg=1)
    X_pro = preprocess(funcs, X)
    print(X_pro)
    
    # 資料切割成訓練集、測試集
    x_train, x_test, y_train, y_test = train_test_split(X_pro, Y, test_size=0.2)
    beta, train_error, test_error = training(x_train, x_test, y_train, y_test)
    print(beta)
    print(train_error)
    print(test_error)

產生結果:

[150.    0.5   0.4   0.7]
1.4244919937188647e-12
1.5027533715667123e-12

參數beta = [150. 0.5 0.4 0.7]
也就是說這個方法有完全猜出來我們設的y = 0.5 * x_0 + 0.4 * x_1 + 0.7 * x_2 + 150的關係式,
並且得到的誤差幾乎是0

用高次多項式誤差會爆炸?

然而,這邊較特別的是,
若我把程式中的funcs = designPolyFunc(3, polyDeg=1)
改成funcs = designPolyFunc(3, polyDeg=2)
則結果變為:

[ 1.12895586e+03  1.07765747e-01 -6.58666127e-01  1.44931910e+00
 -3.56896923e-01  1.81659452e-01 -7.14622243e-01  1.05329766e+00
 -1.08003516e+00  3.67905658e-01]
475110.06876586087
639828.3908878103

整個誤差直接爆炸

這個是overfitting嗎?
應該不算,overfitting是說資料在training data 上表現很好,
但是在testing data上表現不佳,
可是這個是training data和testing data的誤差都炸掉了,
不知是否觀念有想錯,
原因待釐清...

(2020/9/21補充)
誤差爆炸的原因已找到,
如果把X經過資料預處理後的矩陣印出來看的話,
會發現二次方項的多項式數值變的非常大…
可能導致計算上的誤差

[     1    390    927    691 152100 361530 269490 859329 640557 477481]
 [     1      6    670    949     36   4020   5694 448900 635830 900601]
 [     1    157    456    403  24649  71592  63271 207936 183768 162409]
 ...

比較好的方法應該是先將原始資料X做Normalization(見[資料分析&機器學習] 第2.4講:資料前處理(Missing data, One-hot encoding, Feature Scaling)),
將資料縮放至[0,1]區間內,
消除大矩陣運算的誤差

改寫後的程式

先將X做Normalization後,
成功解決誤差炸掉的問題了

from functools import reduce
import numpy as np
import math
from sklearn.model_selection import train_test_split
from itertools import combinations_with_replacement as CWR

"""
 計算一個陣列透過idx的指定乘積
 例:idx=(i1,i2,...in)
 return nums[i1]*nums[i2]*...*nums[in]
"""
def prodOfArr(nums,idxs):
    return reduce(lambda x,y: x*nums[y], idxs, 1)

"""
 生成x的多項式函數,x資料維度=x_dim,次方最高到polyDeg
 例x_dim=4, polyDeg=2
 則X的一筆data形式為 [x1,x2,x3,x4]
 會得到函數: [1,x1,x2,x3,x4, x1*x2, x1*x3,x1*x4,x2*x3,x2*x4,x3*x4]
"""
def designPolyFunc(x_dim, polyDeg=1):
    funcs = [lambda X:1]
    for i in range(1,polyDeg+1):
        funcs += [lambda X, idx=c: prodOfArr(X,idx) for c in CWR(range(x_dim),i)]
    return funcs


"""
 透過函數們funcs將原始資料X=
 x01 x02 ... x0n (= vector x0)
 x11 x12 ... x1n (= vector x1)
 ...
 xm1 xm2 ... xmn (= vector xm)
 其中(m是data數量,n是x的feature數量)
 轉換成矩陣形式:
 1 f1(x0) f2(x0) ... fk(x0)
 1 f1(x1) f2(x1) ... fk(x1)
 ...
 1 f1(xm) f2(xm) ... fk(xm)
"""
def preprocess(funcs, X):
    mat = [[f(x) for f in funcs] for x in X]
    return np.array(mat)

# 公式解B= [(X_TX)^-1+ lamb\*I(單位矩陣)]X_TY
def parameter(X, Y, lamb=0):
    return np.dot(np.linalg.inv(np.dot(X.T,X)+lamb*np.eye(X.shape[1])),np.dot(X.T,Y))

# 計算預測的Y和真實Y之間的root mean square
def RMS(perdict_Y, real_Y):
    return math.sqrt(sum((real_Y-perdict_Y)**2)*2/len(real_Y))


# 把訓練資料、測試資料(預處理過的data)、lamb(若用MAP方法才要傳lamb),做一次訓練
# 回傳值: 訓練得到的最佳參數、訓練資料的RMS、測試資料的RMS
def training(x_train,x_test, y_train, y_test, lamb=0):
    beta = parameter(x_train,y_train, lamb)
    
    ## 以下用train_data做predict,計算train_data的RMS
    train_perdict_Y = np.dot(x_train,beta)
    train_error= RMS(train_perdict_Y,y_train)
    
    ## 以下用test_data做predict,計算test_data的RMS
    perdict_Y = np.dot(x_test,beta)
    test_error= RMS(perdict_Y,y_test)
    return beta, train_error, test_error

def normalization(arr):
    return np.apply_along_axis(lambda x: (x-min(x)) / (max(x)-min(x)), axis=0, arr=arr)
    
if __name__ == '__main__':
    np.random.seed(931) # 這邊設固定的random seed,確保練習用的資料相同
    X = np.random.randint(low=4,high=1000,size=(100,3))
    Y = np.dot(X, np.array([0.5, 0.4, 0.7])) + 150
    
    # 資料預處理
    funcs = designPolyFunc(3, polyDeg=2)
    X = normalization(X)
    X_pro = preprocess(funcs, X)
    
    # 資料切割成訓練集、測試集
    x_train, x_test, y_train, y_test = train_test_split(X_pro, Y, test_size=0.2)
    beta, train_error, test_error = training(x_train, x_test, y_train, y_test)
    print(beta)
    print(train_error)
    print(test_error)


1 則留言

0
z7218758
iT邦新手 5 級 ‧ 2021-06-16 11:21:27

站長您好,
使用您上述的

# 公式解B= [(X_TX)^-1+ lamb\*I(單位矩陣)]X_TY
def parameter(X, Y, lamb=0):
    return np.dot(np.linalg.inv(np.dot(X.T,X)+lamb*np.eye(X.shape[1])),np.dot(X.T,Y))

lamb = 0,得出結果與sklearn linear regression結果相符。
但是lamb =1,得出結果與sklearn ridge regression結果不相同,能否請您幫忙解惑呢?

我要留言

立即登入留言