线性回归

#!/usr/bin/env python
#-*- coding:utf-8 -*-
from numpy import *
import matplotlib.pyplot as plt
"""
打开一个用tab键分隔的文本文件
parameters:
    fileName -文件名
return:
    dataMat -数据矩阵
    labelMat -目标值向量
"""
def loadDataSet(fileName):      
    numFeat = len(open(fileName).readline().split('\t')) - 1 #得到列数，不包括最后一列，默认最后一列值为目标值
    dataMat = []; labelMat = []
    fr = open(fileName)
    for line in fr.readlines():
        lineArr =[]
        curLine = line.strip().split('\t')
        for i in range(numFeat):
            lineArr.append(float(curLine[i]))
        dataMat.append(lineArr)
        labelMat.append(float(curLine[-1]))
    return dataMat,labelMat
"""
计算最佳拟合直线
parameters:
    xArr -给定的输入值
    yArr -给定的输出值
return:
    ws -回归系数
"""
def standRegres(xArr,yArr):
    xMat = mat(xArr); yMat = mat(yArr).T    #将数据保存到矩阵中
    #计算x.T *x
    xTx = xMat.T @ xMat
    #使用linalg.det()方法来判断它的行列式是否为零，即是否可逆
    if linalg.det(xTx) == 0.0:
        return
       
    #使用最小二乘法计算w值
    ws = linalg.inv(xTx) @ xMat.T @ yMat
    return ws
"""
计算局部加权线性回归系数
parameters:
    testPoint -待预测数据
    xArr -给定输入值
    yArr -给定输出值
    k -高斯核的k值，决定对附近的点赋予多大的权重
return:
    testPoint * ws -回归预测的估计值
"""
def lwlr(testPoint, xArr, yArr, k=1.0):
    xMat = mat(xArr); yMat = mat(yArr).T #读入数据到矩阵
    m = shape(xMat)[0]
    #创建对角权重矩阵，该矩阵为方阵，矩阵维数为样本点个数
    theta = eye(m, m)
    #遍历整个数据集
    for i, x in enumerate(xMat):
    #计算待预测数据与每个样本点的差值
        x_i = (testPoint - x) @ (testPoint - x).T
    #计算每个样本点对应的权重值，随着样本点与待预测点距离的递增，权重将以指数级衰减
        theta[i][i] = exp(x_i / (-2 * k ** 2))
    #计算x.T θ x
    xT_theta_x = xMat.T @ theta @ xMat
    #判断矩阵是否可逆
    if linalg.det(xT_theta_x) == 0:
        return
    
    #使用最小二乘法计算w值
    ws = linalg.inv(xT_theta_x) @ (xMat.T @ (theta @ yMat))
    
    return testPoint*ws
"""
测试函数
parameters:
    testArr -测试数据集
    xArr -给定输入值
    yArr -给定输出值
    k -高斯核的k值
return:
    yHat -预测值
"""
def lwlrTest(testArr, xArr, yArr,k=1.0):
    m = shape(xArr)[0]
    yHat = zeros(m)
    for i in range(m):
        yHat[i] = lwlr(testArr[i],xArr,yArr,k)
    return yHat
"""
计算预测误差的平方和
parameters:
    yArr -给定y值
    yHatArr -预测y值
return:
    ((yArr-yHatArr)**2).sum() -误差矩阵
"""
def rssError(yArr,yHatArr):
    return ((yArr-yHatArr)**2).sum()
if __name__=='__main__':
    abX,abY = loadDataSet('linear_regression_abalone/abalone.txt')
    yHat01 = lwlrTest(abX[0:99],abX[0:99],abY[0:99],0.1)
    yHat1 = lwlrTest(abX[0:99],abX[0:99],abY[0:99],1)
    yHat10 = lwlrTest(abX[0:99],abX[0:99],abY[0:99],10)
    print("使用局部加权线性回归预测误差：")
    print("核为0.1时：",rssError(abY[0:99],yHat01.T))
    print("核为1时：",rssError(abY[0:99],yHat1.T))
    print("核为10时：",rssError(abY[0:99],yHat10.T))
    yHat01 = lwlrTest(abX[100:199],abX[0:99],abY[0:99],0.1)
    yHat1 = lwlrTest(abX[100:199],abX[0:99],abY[0:99],1)
    yHat10 = lwlrTest(abX[100:199],abX[0:99],abY[0:99],10)
    print("使用局部加权线性回归预测误差在新数据上的表现：")
    print("核为0.1时：",rssError(abY[100:199],yHat01.T))
    print("核为1时：",rssError(abY[100:199],yHat1.T))
    print("核为10时：",rssError(abY[100:199],yHat10.T))
    ws = standRegres(abX[0:99],abY[0:99])
    yHat = mat(abX[100:199])*ws
    print("使用标准线性回归预测误差为：",rssError(abY[100:199],yHat.T.A))

import numpy as np

import matplotlib.pyplot as plt

  

np.random.seed(42)

w = np.array([2, 1, 4, 5, 3])

d = len(w)

X = []

Y = []

for _ in range(1000000):

    x = np.random.randn(d)

    y = w.dot(x) + np.random.randn()

    X.append(x)

    Y.append(y)

X = np.array(X)

Y = np.array(Y)

  

def mse(y_true, y_test):

    return ((y_true - y_test) ** 2) / len(y_true)

  

def gradient(y_true, y_test):

    return 2 * (y_test - y_true) / len(y_true)

  

def batch_gradient_descent(w, alpha, x, y):

    y_pred = x.dot(w)

    error = mse(y, y_pred).mean()

    grad = np.dot(x.T, gradient(y, y_pred))

    w = w - alpha * grad

  

    return w, error

  

def stochastic_gradient_descent(w, alpha, x, y, epoch):

    alpha_update = alpha

    for i in range(len(x)):

        y_pred = x[i].dot(w)

        grad = np.dot(x[i].T, (y_pred - y[i])) * 2 / len(x)

        w = w- alpha_update * grad

        alpha_update = alpha_update / (epoch+1)

    error = mse(y, x.dot(w)).mean()

  

    return w, error

  
  

X_test = []

Y_test = []

for _ in range(10000):

    x = np.random.randn(d)

    y = w.dot(x) + np.random.randn()

    X_test.append(x)

    Y_test.append(y)

  

X_test = np.array(X_test)

Y_test = np.array(Y_test)

  
  
  

def l2_mse(y_true, y_test, l, w):

    return ((y_true - y_test) ** 2) / len(y_true) + l * np.sum(w ** 2)

  

def l2_gradient(y_true, y_test):

    return 2 * (y_test - y_true) / len(y_true)

  

def batch_gradient_descent_with_l2(w, alpha, x, y, l):

    y_pred = x.dot(w)

    error = l2_mse(y, y_pred, l, w).mean()

    grad = np.dot(x.T, l2_gradient(y, y_pred))

    w = w - alpha * grad - alpha * l * w *2

  

    return w, error

  
  
  

if __name__ == "__main__":

    train_loss = []

    test_loss = []

    print("Batch Gradient Descent")

    for epoch in range(1000):

        w, error = batch_gradient_descent(w, 0.01, X, Y) # train

        y_pred = X_test.dot(w)

        error_test = mse(Y_test, y_pred).mean() # test

        if epoch % 100 == 0:

            print("Epoch: {}, TrainError: {}, TestError: {}".format(epoch, error, error_test))

  

        train_loss.append(error)

        test_loss.append(error_test)

    plt.plot(train_loss, label="Train-No-L2")

    plt.legend()

    plt.xlabel("Epoch")

    plt.ylabel("Loss")

    plt.show()

  

    plt.plot(test_loss, label="Test-No-L2")

    plt.xlabel("Epoch")

    plt.ylabel("Loss")

    plt.legend()

    plt.show()

    plt.plot(train_loss, label="Train-No-L2")

    plt.plot(test_loss, label="Test-No-L2")

    plt.legend()

    plt.show()

    # ============================================

    train_loss = []

    test_loss = []

    print("Batch Gradient Descent with L2")

    l = 0.0001 # lambda

    for epoch in range(1000):

        w, error = batch_gradient_descent_with_l2(w, 0.01, X, Y, l) # train

        y_pred = X_test.dot(w)

        error_test = l2_mse(Y_test, y_pred, l, w).mean() # test

        if epoch % 100 == 0:

            print("Epoch: {}, TrainError: {}, TestError: {}".format(epoch, error, error_test))

  

        train_loss.append(error)

        test_loss.append(error_test)

    plt.plot(train_loss, label="Train-L2")

    plt.legend()

    plt.xlabel("Epoch")

    plt.ylabel("Loss")

    plt.show()

  

    plt.plot(test_loss, label="Test-L2")

    plt.xlabel("Epoch")

    plt.ylabel("Loss")

    plt.legend()

    plt.show()

    plt.plot(train_loss, label="Train-L2")

    plt.plot(test_loss, label="Test-L2")

    plt.legend()

    plt.show()