【邱希鵬】nndl-chap2-linear_regression

線性回歸

問題描述

有一個函數 ，使得。現 在不知道函數 \(f(\cdot)\)的具體形式，給定滿足函數關系的一組訓練樣本，請使用線性回歸模型擬合出函數\(y=f(x)\)。

(可嘗試一種或幾種不同的基函數，如多項式、高斯或sigmoid基函數）

數據集

根據某種函數關系生成的train 和test 數據。

題目要求：

• [ ] 按順序完成 exercise-linear_regression.ipynb中的填空

  1. 先完成最小二乘法的優化 (參考書中第二章 2.3節中的公式)
  2. 附加題：實現“多項式基函數”以及“高斯基函數”（可參考PRML）
  3. 附加題：完成梯度下降的優化 (參考書中第二章 2.3節中的公式)

• [ ] 參照lienar_regression-tf2.0.ipnb使用tensorflow2.0 使用梯度下降完成線性回歸

• [ ] 使用訓練集train.txt 進行訓練，使用測試集test.txt 進行評估（標准差），訓練模型時請不要使用測試集。

exercise-linear_regression

import numpy as np
import matplotlib.pyplot as plt

def load_data(filename):
    """載入數據。"""
    xys = []
    with open(filename, 'r') as f:
        for line in f:
            xys.append(map(float, line.strip().split()))
        xs, ys = zip(*xys)
        return np.asarray(xs), np.asarray(ys)

不同的基函數 (basis function)的實現

def identity_basis(x):
    ret = np.expand_dims(x, axis=1)
    return ret

def multinomial_basis(x, feature_num=10):
    x = np.expand_dims(x, axis=1) # shape(N, 1)
    feat = [x]
    for i in range(2, feature_num+1):
        feat.append(x**i)
    ret = np.concatenate(feat, axis=1)
    return ret

def gaussian_basis(x, feature_num=10):
    centers = np.linspace(0, 25, feature_num)
    width = 1.0 * (centers[1] - centers[0])
    x = np.expand_dims(x, axis=1)
    x = np.concatenate([x]*feature_num, axis=1)
    
    out = (x-centers)/width
    ret = np.exp(-0.5 * out ** 2)
    return ret

返回一個訓練好的模型 填空順序 1 用最小二乘法進行模型優化

先完成最小二乘法的優化 (參考書中第二章 2.3中的公式)

再完成梯度下降的優化 (參考書中第二章 2.3中的公式)

在main中利用訓練集訓練好模型的參數，並且返回一個訓練好的模型。

計算出一個優化后的w，請分別使用最小二乘法以及梯度下降兩種辦法優化w。

• 最小二乘法

def main(x_train, y_train):
    """
    訓練模型，並返回從x到y的映射。
    
    """
    basis_func = identity_basis                             # shape(N, 1)的函數
    phi0 = np.expand_dims(np.ones_like(x_train), axis=1)    # shape(N,1)大小的全1 array
    phi1 = basis_func(x_train)                              # 將x_train的shape轉換為(N, 1)
    phi = np.concatenate([phi0, phi1], axis=1)              # phi.shape=(300,2) phi是增廣特征向量的轉置
    
    
    #==========
    #todo '''計算出一個優化后的w，請分別使用最小二乘法以及梯度下降兩種辦法優化w'''
    #==========              
    w = np.dot(np.linalg.pinv(phi), y_train)      # np.linalg.pinv(phi)求phi的偽逆矩陣(phi不是列滿秩) w.shape=[2,1]
    
    def f(x):
        phi0 = np.expand_dims(np.ones_like(x), axis=1)
        phi1 = basis_func(x)
        phi = np.concatenate([phi0, phi1], axis=1).T
        
        y = np.dot(w, phi)
        return y

    return f

• 梯度下降方法:

def main(x_train, y_train):
    """
    訓練模型，並返回從x到y的映射。
    """
    basis_func = multinomial_basis
    phi0 = np.expand_dims(np.ones_like(x_train), axis=1)      #phi0.shape=(300,1)
    phi1 = basis_func(x_train)                                #phi1.shape=(300,1)

    phi = np.concatenate([phi0, phi1], axis=1)                #phi.shape=(300,2)phi是增廣特征向量的轉置

    ### START CODE HERE ###
    #梯度下降法
    def dJ(theta, phi, y):
        return phi.T.dot(phi.dot(theta)-y)*2.0/len(phi)

    def gradient(phi, y, initial_theta, eta=0.001, n_iters=10000):
        w = initial_theta

        for i in range(n_iters):
            gradient = dJ(w, phi, y)                #計算梯度
            w = w - eta *gradient                   #更新w

        return w

    initial_theta = np.zeros(phi.shape[1])
    w = gradient(phi, y_train, initial_theta)
    ### END CODE HERE ###

    def f(x):
        phi0 = np.expand_dims(np.ones_like(x), axis=1)
        phi1 = basis_func(x)
        phi = np.concatenate([phi0, phi1], axis=1)
        y = np.dot(phi, w)
        return y

    return f

def evaluate(ys, ys_pred):
    """評估模型。"""
    std = np.sqrt(np.mean(np.abs(ys - ys_pred) ** 2))
    return std

# 程序主入口（建議不要改動以下函數的接口）
if __name__ == '__main__':
    train_file = 'train.txt'
    test_file = 'test.txt'
    # 載入數據
    x_train, y_train = load_data(train_file)
    x_test, y_test = load_data(test_file)
    print(x_train.shape)
    print(x_test.shape)
    
    
    # 使用線性回歸訓練模型，返回一個函數f()使得y = f(x)
    f = main(x_train, y_train)

    y_train_pred = f(x_train)
    std = evaluate(y_train, y_train_pred)
    print('訓練集預測值與真實值的標准差：{:.1f}'.format(std))
    
    # 計算預測的輸出值
    y_test_pred = f(x_test)
    # 使用測試集評估模型
    std = evaluate(y_test, y_test_pred)
    print('預測值與真實值的標准差：{:.1f}'.format(std))

    #顯示結果
    plt.plot(x_train, y_train, 'ro', markersize=3)
#     plt.plot(x_test, y_test, 'k')
    plt.plot(x_test, y_test_pred, 'k')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('Linear Regression')
    plt.legend(['train', 'test', 'pred'])
    plt.show()

運行結果：（調用identity_basis）
(300,)
(200,)
訓練集預測值與真實值的標准差：2.0
預測值與真實值的標准差：2.1

運行結果：（調用multinomial_basis）
(300,)
(200,)
訓練集預測值與真實值的標准差：2.0
預測值與真實值的標准差：2.2

運行結果：（調用gaussian_basis）
(300,)
(200,)
訓練集預測值與真實值的標准差：0.4
預測值與真實值的標准差：0.4

linear_regression_tensorflow2.0

設計基函數(basis function) 以及數據讀取

import numpy as np
import matplotlib.pyplot as plt

def identity_basis(x):
    ret = np.expand_dims(x, axis=1)
    return ret

def multinomial_basis(x, feature_num=10):
    x = np.expand_dims(x, axis=1) # shape(N, 1)
    feat = [x]
    for i in range(2, feature_num+1):
        feat.append(x**i)
    ret = np.concatenate(feat, axis=1)
    return ret

def gaussian_basis(x, feature_num=10):
    centers = np.linspace(0, 25, feature_num)
    width = 1.0 * (centers[1] - centers[0])
    x = np.expand_dims(x, axis=1)
    x = np.concatenate([x]*feature_num, axis=1)
    
    out = (x-centers)/width
    ret = np.exp(-0.5 * out ** 2)
    return ret

def load_data(filename, basis_func=gaussian_basis):
    """載入數據。"""
    xys = []
    with open(filename, 'r') as f:
        for line in f:
            xys.append(map(float, line.strip().split()))
        xs, ys = zip(*xys)
        xs, ys = np.asarray(xs), np.asarray(ys)
        
        o_x, o_y = xs, ys
        phi0 = np.expand_dims(np.ones_like(xs), axis=1)
        phi1 = basis_func(xs)
        xs = np.concatenate([phi0, phi1], axis=1)
        return (np.float32(xs), np.float32(ys)), (o_x, o_y)

定義模型

import tensorflow as tf
from tensorflow.keras import optimizers, layers, Model

class linearModel(Model):
    def __init__(self, ndim):
        super(linearModel, self).__init__()
        self.w = tf.Variable(
            shape=[ndim, 1], 
            initial_value=tf.random.uniform(
                [ndim,1], minval=-0.1, maxval=0.1, dtype=tf.float32))
        
    @tf.function
    def call(self, x):
        y = tf.squeeze(tf.matmul(x, self.w), axis=1)
        return y

(xs, ys), (o_x, o_y) = load_data('train.txt')        
ndim = xs.shape[1]

model = linearModel(ndim=ndim)

訓練以及評估

optimizer = optimizers.Adam(0.1)
@tf.function
def train_one_step(model, xs, ys):
    with tf.GradientTape() as tape:
        y_preds = model(xs)
        loss = tf.reduce_mean(tf.sqrt(1e-12+(ys-y_preds)**2))
    grads = tape.gradient(loss, model.w)
    optimizer.apply_gradients([(grads, model.w)])
    return loss

@tf.function
def predict(model, xs):
    y_preds = model(xs)
    return y_preds

def evaluate(ys, ys_pred):
    """評估模型。"""
    std = np.sqrt(np.mean(np.abs(ys - ys_pred) ** 2))
    return std

for i in range(1000):
    loss = train_one_step(model, xs, ys)
    if i % 100 == 1:
        print(f'loss is {loss:.4}')
        
        
y_preds = predict(model, xs)
std = evaluate(ys, y_preds)
print('訓練集預測值與真實值的標准差：{:.1f}'.format(std))

(xs_test, ys_test), (o_x_test, o_y_test) = load_data('test.txt')

y_test_preds = predict(model, xs_test)
std = evaluate(ys_test, y_test_preds)
print('訓練集預測值與真實值的標准差：{:.1f}'.format(std))

plt.plot(o_x, o_y, 'ro', markersize=3)
plt.plot(o_x_test, y_test_preds, 'k')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Linear Regression')
plt.legend(['train', 'test', 'pred'])
plt.show()

loss is 11.67
loss is 1.655
loss is 1.608
loss is 1.572
loss is 1.533
loss is 1.495
loss is 1.454
loss is 1.411
loss is 1.366
loss is 1.32
訓練集預測值與真實值的標准差：1.5
訓練集預測值與真實值的標准差：1.8