1

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

Linear Regression的核心概念是
Y的值跟m個特徵有關，

`y = w1*f1(x)+ w2*f2(x)+...+ wM*fM(x)`

w1, w2, ..., wM是我們想要求解的參數

# 直接用sklearn實作的版本

## sklearn程式例子

w0, w1, b 是我們要求解的參數

``````import numpy as np
from sklearn.linear_model import LinearRegression

"""

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跟X都是已知，是從訓練資料來的，
B是我們要求解的參數

## 公式解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))
``````

## 資料轉換

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

``````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)
``````

``````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做預測

``````# 計算預測的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
``````

# 程式實作範例

``````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 ...]
``````

## 完整程式示例

``````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
``````

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

``````[ 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
``````

(2020/9/21補充)

``````[     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]
...
``````

## 改寫後的程式

``````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結果相符。