邱錫鵬nndl-chap2_1-線性回歸

1.問題描述

有一個函數$f:R\rightarrow R$ ，現在不知道函數 $f(x)$的具體形式，給定滿足函數關系的一組訓練樣本$\left \{ (x_{1},y_{1}),...,(x_{N},y_{N}) \right \}$，N=300，請使用線性回歸模型擬合出函數$y=f(x)$。(可嘗試一種或幾種不同的基函數，如多項式、高斯或sigmoid基函數）

2.數據集

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

3.代碼實現

3.1載入數據

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)                   #nsarray將結構數據轉換為ndarray類型

3.2不同基函數的實現

對於函數空間中的連續函數都可以表示成一系列基函數的線性組合，就像是在向量空間中每個向量都可以表示成基向量的線性組合一樣。

舉例：多項式基：{1，t, t²}是實系數二次多項式集合的基，每一個形如a+bt+ct²的二次多項式都可以寫成由基函數1、t、t²組成的線性組合。

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


def multinomial_basis(x, feature_num=10):
    '''多項式基函數'''
    x = np.expand_dims(x, axis=1)                
    
    ### START CODE HERE ###
    feat = [x]
    for i in range(2, feature_num+1):
        feat.append(x**i)
    ret = np.concatenate(feat, axis=1) 
    ### END CODE HERE ###
    
    return ret


def gaussian_basis(x, feature_num=10):
    '''高斯基函數'''
    ### START CODE HERE ###
    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)
    ### END CODE HERE ###

    return ret

3.3訓練模型

phi是增廣矩陣向量的轉置，不同的基函數會形成不同的phi，以identity_basic為例，$phi=\begin{bmatrix}1. & X^{(1)}\\1. & X^{(2)}\\... & \\ 1. & X^{(300)}\end{bmatrix}$

解法1：最小二乘法求解線性回歸參數：$w=(XX^{T})^{-1}Xy$       $(XX^{T})^{-1}X$也稱為X^T的偽逆矩陣

$y=phi*w$

def main(x_train, y_train):
    """
    訓練模型，並返回從x到y的映射。
    """
    basis_func = identity_basic
    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 ###
    #最小二乘法計算w
    w = np.dot(np.linalg.pinv(phi), y_train)                  #np.linalg.pinv(phi)求phi的偽逆矩陣(phi不是列滿秩) w.shape=[2,1]
    ### 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

解法2：梯度下降求解線性回歸參數：$w\leftarrow w+\alpha X(y-X^{T}w)$

詳細推導見https://www.cnblogs.com/cxq1126/p/13051607.html

def main(x_train, y_train):
    """
    訓練模型，並返回從x到y的映射。
    """
    basis_func = identity_basic
    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

3.4評估模型

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

完整代碼如下（使用identity_basis基函數和最小二乘法為例）：

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)                   #nsarray將結構數據轉換為ndarray類型


#不同的基函數的實現
def identity_basis(x):
    ret = np.expand_dims(x, axis=1)                             #shape(N, 1)
    return ret


def multinomial_basis(x, feature_num=10):
    '''多項式基函數'''
    x = np.expand_dims(x, axis=1)                
    
    ### START CODE HERE ###
    feat = [x]
    for i in range(2, feature_num+1):
        feat.append(x**i)
    ret = np.concatenate(feat, axis=1) 
    ### END CODE HERE ###
    
    return ret


def gaussian_basis(x, feature_num=10):
    '''高斯基函數'''
    ### START CODE HERE ###
    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)
    ### END CODE HERE ###

    return ret


def main(x_train, y_train):
    """
    訓練模型，並返回從x到y的映射。
    """
    basis_func = identity_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 ###
    '''使用最小二乘法計算w'''
    w = np.dot(np.linalg.pinv(phi), y_train)                  #np.linalg.pinv(phi)求phi的偽逆矩陣(phi不是列滿秩)
    ### 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 = 'data/train.txt'
    test_file = 'data/test.txt'
    # 載入數據
    x_train, y_train = load_data(train_file)
    x_test, y_test = load_data(test_file)
   
    print(x_train.shape)                     #(300,)
    print(x_test.shape)                      #(200,)

    # 使用線性回歸訓練模型，返回一個函數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()

4.運行結果

以下三張圖是分別調用identity_basis、multinomial_basis、gaussian_basis三種不同基函數（最小二乘法計算的w）運行而來。

訓練集預測值與真實值的標准差：2.0                                         訓練集預測值與真實值的標准差：1.5 
預測值與真實值的標准差：2.2                                                    預測值與真實值的標准差：1.6     

訓練集預測值與真實值的標准差：0.4

預測值與真實值的標准差：0.4 

以下二張圖是分別調用identity_basis、gaussian_basis二種不同基函數（梯度下降法計算的w）運行而來。

訓練集預測值與真實值的標准差：2.0                                                   訓練集預測值與真實值的標准差：1.9
預測值與真實值的標准差：2.2                                                              預測值與真實值的標准差：2.1