机器学习线性回归递归下降算法实现

代码：https://github.com/Jasonssssss9/Machine-Learning/tree/main/LinearRegression/Implementation_1

一、基本原理

1.线性回归原理

线性回归是一种常见的有监督机器学习算法，用于预测一个连续值的输出

线性回归预测函数：$$h_\theta(x^i)=\theta_1x_1^i+\theta_2x_2^i+...+\theta_nx_n^i+\theta_0=\sum_{j=0}^{n}\theta_jx_j^i=\theta^Tx^i$$ 其中x^i为第i个样本，x_j为x的第j维分量(第j个特征值)，θj为参数θ的第j维分量

写为矩阵形式：$$h_\theta(x)=X\theta$$（这里的X矩阵增加了一列全1） θ为一个n维列向量，n为特征值个数；X为一个m×n矩阵，m为元素个数

• 模型训练：根据训练集数据特征值X和目标值y，借助最大似然估计法+梯度下降优化方法(也可以采用其他优化方法)得到参数θ值

• 模型预测：将测试集数据特征值X输入预测函数，得到预测值，并且与真实值进行比较

2.线性回归模型误差项

每个预测值都有误差，记第i个数据误差为ε^i，则对每个样本，均有：$$y^i=\theta^Tx^i+\varepsilon^i$$ ① θ 和x^i为n维列向量

对于误差需要进行如下假设：所有误差均为独立同分布且满足正态分布(高斯分布)，均值为0，方差为σ^2

从而可得σ^i的概率密度函数：$$P(\varepsilon^i)=\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(\varepsilon^i)^2}{2\sigma^2})$$ ②

联立①②式，可得：$$P(y^i|x^i;\theta)=\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(y^i-\theta^Tx^i)^2}{2\sigma^2})$$ 理解：关键是x^i和θ是给定的，这样y^i才有取值的可能，因此$$\theta^Tx^i$$实际上是一个定值 均值$$E(\varepsilon ^i)=0 => E(\varepsilon^i+C)=C$$，并且加上常数方差不变 即$$y^i=\theta^Tx^i+\varepsilon^i$$ 也满足正态分布，均值为$$\theta^Tx^i$$，方差为σ^2，从而可得上式

3.最大似然估计求目标函数

带入似然函数公式，得似然函数：$$L(\theta )=\prod_{i=1}^{m}P(y^i|x^i;\theta )$$

对数似然函数：$$ln(L(\theta))=\sum_{i=1}^{m}ln(\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(y^i-\theta^Tx^i)^2}{2\sigma^2}))$$

展开得：$$ln(L(\theta))=mln\frac{1}{\sqrt{2\pi}\sigma}-\frac{1}{\sigma^2} \frac{1}{2}\sum_{i=1}^{m} (y^i-\theta^Tx^i)^2$$

求L(θ)的最大值，实际上就是求ln(L(θ))的最大值，也实际上就是求$$\frac{1}{2}\sum_{i=1}^{m} (y^i-\theta^Tx^i)^2$$的最小值

因此目标函数为$$J(\theta)=\frac{1}{2}\sum_{i=1}^{m}(h_\theta(x^i)-y^i)^2=\frac{1}{2}(X\theta-y)^T(X\theta-y)$$

重要处理思路：将求和\sum转换为向量乘法

平方按上式思路转化为一个列向量转置后自乘自身。如果是一次方，可以引入一个全1的向量

这里可以采用求导为0的方法直接求θ

 

令上式为0，得$$\theta=(X^TX)^{-1}X^Ty$$

但是当数据量较大时，求X转置，求逆，复杂度非常高，效率非常低，并且还可能出现X^TX的逆矩阵不存在的情况，因此最好使用优化方法求解

4.梯度下降法求解

目标函数：$$J(\theta)=\frac{1}{2m}\sum_{i=1}^{m}(h_\theta(x^i)-y^i)^2$$ 加了一个m，表示平均值，不影响结果θ的取值

梯度：$$\frac{\partial J(\theta)}{\partial \theta_j} =-\frac{1}{m}\sum_{i=1}^{m}(y^i-h_\theta(x^i))x_j^i$$（其中j表示第j维）

从而得到迭代公式：$${\theta_j}'=\theta_j-\alpha\frac{1}{num}\sum_{k=i}^{i+num}(h_{\theta_j}(x^k)-y^k)x_j^k$$

• num=m，则选取所有样本进行梯度下降，虽然准确，但如果样本数量过大，则会太慢

  这种方式被称为批量梯度下降或全梯度下降，数据较少时选择此法

• num=1，随机选取一个数据进行梯度下降，被称为随机梯度下降，虽然快但不准，不会选择此法

• 选择一部分数据进行计算，num为这部分数据的个数，这部分数据从i到i+num

  这种方式为小批量梯度下降法，数据量较大时选取此方法

θ'在θ的基础上进行迭代，不仅需要梯度，还需要定义一个步长(学习率)参数α，表示一次走多少

实际上是一个凸优化模型：

• $$\frac{\partial J(\theta)}{\partial \theta_j}$$的正负决定迭代方向

  • 小于0时，函数图像下降，θ' > θ

  • 大于0时，函数图像上升，θ' < θ

• $$|\frac{\partial J(\theta)}{\partial \theta_j}|$$和α决定一次走多远

  • 绝对值较大，说明离目标远，一次走的多

  • 绝对值较小，说明离目标近，一次走的少

 

二、numpy代码实现

本实验以python numpy库为基础，借助sklearn.preprocessing进行数据预处理

1.训练和预测模型

在LinearRegression.py文件中定义LinearRegression类，实现梯度下降算法

• 初始化

  class LinearRegression:
      '''
      一个线性回归机器学习模型
      主要分为两个功能：(1)传入训练集进行模型训练, 得到θ参数; (2)传入测试集进行预测
      '''
  ​
      def  __init__(self, features, targets):
          # 预处理特征值:
          features  = self.__preprocessing(features)
  ​
          self.features = features #特征值
          self.targets = targets   #目标值
          self.theta = np.zeros(self.features.shape[1]) #初始化θ值

  主要设置features(特征值)，targets(目标值)和theta(参数)为成员变量，首先需要调用__proprocessing函数对features对特征值进行预处理

• 预处理

  def __preprocessing(self, data):
      '''
      预处理数据
      '''
      scaler = StandardScaler()
      data = scaler.fit_transform(data)
      # 给features增加一列1
      num_rows = data.shape[0]
      ones_col = np.ones((num_rows,1))
      data = np.hstack((data, ones_col))
      return data

  使用sklearn提供的StandardScaler进行数据标准化，即features先减去均值再除以标准差，让数据分布在原点周围，并且具有单位标准差

  之后再给features增加一列1 因为最终的线性回归方程为y=aX+b，增加一列1保证可以获得常数项

  注意，这里很容易忘记常数项

• 训练模型

  • train函数

    def train(self, alpha = 0.01, num_iterations = 500):
        '''
        传入学习率和迭代次数参数，进行模型训练
        返回最终得到的theta，以及迭代过程中每一步的均方误差作为测试信息
        '''
        mse_history = self.__gradient_descent(alpha, num_iterations)
        return self.theta, mse_history

    train为暴露给外部的接口，实现模型训练，需要传入学习率(步长)α和最大迭代次数

    返回参数theta，以及一个mse_history列表

    这个mse_history列表存储迭代中每一步的均方误差

  • __gradient_descent函数

    def __gradient_descent(self, alpha, num_iterations):
        '''
        实际进行梯度下降算法
        生成最终的θ
        返回迭代中每一步的均方误差，用来在测试中显示迭代过程
        '''
        mse_history = []
        for _ in range(num_iterations):
            self.__gradient_step(alpha) #更新θ
            mse_history.append(self.__current_MSE())
        return mse_history

    进行梯度下降算法，迭代num_iterations次，每一步调用__gradient_step函数更新theta，再调用__current_MSE函数获取当前均方误差，加入到mes_history中

  • __gradient_step和__curr_prediction函数

    def __curr_prediction(self, data, theta):
        '''
        计算预测结果
        '''
        return np.dot(data, theta)
    ​
    def __gradient_step(self, alpha):
        '''
        实现每一步计算θ的公式
        '''
        curr_prediction = self.__curr_prediction(self.features, self.theta)
        delta = curr_prediction - self.targets
        self.theta = self.theta - alpha*(1/self.features.shape[0])*(np.dot(delta.T, self.features)).T

    以上即矩阵实现梯度下降法的公式

  • __current_MSE函数

    def __current_MSE(self):
        '''
        计算当前模型的均方误差(Mean Squared Error)
        '''
        curr_prediction = self.__curr_prediction(self.features, self.theta)
        delta = curr_prediction - self.targets
        mse = (1/2)*np.dot(delta.T, delta)/self.features.shape[0]
        return mse

    计算当前均方误差

• 使用模型进行预测

  def predict(self, data):
      '''
      根据训练得到的模型进行数据预测
      '''
      #预处理数据：
      data = self.__preprocessing(data)
  ​
      return self.__curr_prediction(data, self.theta)

  传入需要预测的数据(测试集数据)，返回预测结果

2.测试1

# test_1.py
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
​
from LinearRegression import LinearRegression
​
# 获取世界幸福指数数据
data = pd.read_csv('world-happiness-report-2017.csv')
​
train_data = data.sample(frac = 0.8)
test_data = data.drop(train_data.index)
​
# 为了方便打印结果，这里特征值只选取一列数据
features_name = 'Economy..GDP.per.Capita.'
target_name = 'Happiness.Score'
​
# 获取特征值和目标值
X_train = train_data[[features_name]].values
# 注意这里要用两个[]，因为即使只有一列数据，X也应该是一个二维数组
y_train = train_data[target_name].values
X_test = test_data[[features_name]].values
y_test = test_data[target_name].values
​
# 进行模型训练
linear_regression = LinearRegression(X_train,y_train)
(theta,mse_history) = linear_regression.train(0.01, 1000)
​
# 打印初始和最终均方误差
print ('My initial mse:',mse_history[0])
print ('My final mse:',mse_history[-1])
​
# 打印均方误差和迭代次数图像
plt.plot(mse_history)
plt.xlabel('Iteration')
plt.ylabel('MSE')
plt.title('Iteration-MSE Diagram')
plt.show()
​
predictions_num = 100
x_predictions = np.linspace(X_train.min(),X_train.max(),predictions_num).reshape(predictions_num,1)
y_predictions = linear_regression.predict(x_predictions)
​
plt.scatter(X_train,y_train,label='Train data')
plt.scatter(X_test,y_test,label='Test data')
plt.plot(x_predictions,y_predictions,'r',label = 'Prediction')
plt.xlabel(features_name)
plt.ylabel(target_name)
plt.title('Linear Regression Diagram')
plt.legend()
plt.show()

步骤：

①获取数据；②对数据进行划分；③进行模型训练；④打印均方误差查看训练效果；⑤进行数据预测；⑥打印数据预测结果

注意：这里实际上没有用到测试集，直接通过X_train.min()和X_train_max()生成一个新的100个测试数据

结果：

MSE迭代减小图像：

拟合结果

打印结果

3.测试2

测试1的数据较为简单，并且只选取了特征值的一列。测试2借助sklearn提供的糖尿病数据，将本代码和sklearn提供的线性回归递归下降算法的预测结果进行比较

# test_2.py
​
from sklearn.linear_model import SGDRegressor
from sklearn.datasets import load_diabetes
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
​
from LinearRegression import LinearRegression
​
​
# 加载糖尿病数据集
diabetes = load_diabetes()
X = diabetes.data
y = diabetes.target
​
# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
​
​
# 使用sklearn中的递归下降算法进行模型训练
model = SGDRegressor(eta0=0.01, max_iter=10000, tol=1e-3)
model.fit(X_train, y_train)
predictions_1 = model.predict(X_test)
​
# 计算均方误差
mse_1 = mean_squared_error(y_test, predictions_1)
print("Sklearn mse:", mse_1)
​
​
# 使用LinearRegression的方法进行模型训练
linear_regression = LinearRegression(X_train, y_train)
(theta, mse_history) = linear_regression.train(0.01, 1000)
​
# 打印一个mes_history
plt.plot(mse_history)
plt.xlabel('Iteration')
plt.ylabel('MSE')
plt.title('Iteration-MSE diagram')
plt.show()
​
predictions_2 = linear_regression.predict(X_test)
mse_2 = mean_squared_error(y_test, predictions_2)
print("My mse:", mse_2)

 

可以看到最终结果差异不大