【机器学习】线性回归

一、问题概述

• 线性：变量之间的关系是一次函数关系，例如直线、平面等。相反，若变量之间的关系不是一次函数关系，称之为非线性。
• 回归：通过有限的条件计算所得的往往是测量值，为了得到真实值，需要进行无限次测量，通过测试数据接近到真实值，这就是回归。

综上，线性回归问题 即通过一个线性模型（一次函数）尽可能准确地预测新样本的输出，通过已知数据预测未知结果。例如：房价预测、电影票房预估等。

1、线性模型

例如：现在有两个变量$x_1、x_2$，假设结果为$h_{\theta }(x)$，则有$$h_{\theta }(x)={\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2$$。
其中 ${\theta }_0$ 称为偏置项，${\theta }_1、{\theta }_2$称为权重。

若有$n$个变量，令偏置项 ${\theta }_0$ 的变量值$x_0=1$，则预测值为$$h_{\theta }(x)=\sum _{i=0}^n{\theta }_i{x}_i={\theta }^{T}x$$。

其中 $\theta =[{\theta }_0,{\theta }_1,{\theta }_2]$ 为一个行矩阵。

2、误差损失

假设一个线性模型：

可以看出预测值与真实值之间并不是完全对等的，因此引入误差 $\epsilon$。

对于每个样本：$y^{(i)}={\theta }^T{x}^{(i)}+{\epsilon }^{(i)}$。${\epsilon }^{(i)}$称为误差项，不同样本的误差项也不一定相同。

误差 ${\epsilon }^{(i)}$ 符合独立同分布并且服从均值为$0$、方差为${\theta }^2$的高斯分布（正态分布）。

3、似然函数

由于误差服从高斯分布，则有：
$$p({\epsilon }^{(i)})=\frac{1}{\sqrt{2\pi \sigma }}\exp (-\frac{({\epsilon }^{(i)}{)}^2}{2{\sigma }^2})$$
将 $y^{(i)}={\theta }^T{x}^{(i)}+{\epsilon }^{(i)}$ 代入得到：
$$p(y^{(i)}|x^{(i)};\theta )=\frac{1}{\sqrt{2\pi \sigma }}\exp (-\frac{(y^{(i)}-{\theta }^T{x}^{(i)}{)}^2}{2{\sigma }^2})$$
因此 似然函数 为：
$$L(\theta )=\prod _{i=1}^{m}p(y^{(i)}|x^{(i)};\theta )=\prod _{i=1}^m\frac{1}{\sqrt{2\pi \sigma }}\exp (-\frac{(y^{(i)}-{\theta }^T{x}^{(i)}{)}^2}{2{\sigma }^2})$$
转化为对数似然：
$$\begin{array}{rl}\ln L(\theta )& =\sum _{i=1}^m\ln \frac{1}{\sqrt{2\pi \sigma }}\exp (-\frac{(y^{(i)}-{\theta }^T{x}^{(i)}{)}^2}{2{\sigma }^2})\\ & =m\ln \frac{1}{\sqrt{2\pi \sigma }}-\frac{1}{{\sigma }^2}\cdot \frac{1}{2}\sum _{i=1}^m(y^{(i)}-{\theta }^T{x}^{(i)}{)}^2\end{array}$$
令
$$J(\theta )=\frac{1}{2}\sum _{i=1}^m(y^{(i)}-{\theta }^T{x}^{(i)}{)}^2（最小二乘法）$$
$J(\theta )$即损失函数 ，由于似然函数来源于高斯分布，则应使似然函数越大越好，即 $J(\theta )$越小越好。

怎么让$J(\theta )$取得最小值呢？——求导：

注：这一部分不需要掌握，当时需要理解。

4、梯度下降

对于最小值的求解问题，我们可以将其近似看做一个小山，那么局部最小值就变成了下山问题：

由于 梯度（切线）方向总是下降的，因此我们每次沿着梯度方向前进一小步逼近局部最优解的过程，即称为 梯度下降。

那步长怎么取呢？
  对于最优解的求解，我们的惯性思维自然是越快越好，即将步长设置成一个较大值，一般情况下可以很快的逼近局部最优解。但是，如果在一个逼近局部最优解的小区间，由于步长较大，求解过程可能在最优解附近反复横跳，将永远不会收敛，即永远取不到局部最优解。因此，我们的步长不应取得很大，也不宜取得很小，因此如果步长过小，逼近局部最优解得过程将会很慢，时间利用率不高。

由于权值及偏置项之间是独立的，因此我们分开求解：
对于一个 ${\theta }_j$，我们关于$J(\theta )$对其求偏导并取平均:
$$\frac{\partial J(\theta )}{\partial {\theta }_j}=-\frac{1}{m}\sum _{i=1}^m(y^i-h_{\theta }(x^i))x_j^i$$
而梯度是下降的，因此我们取步长 $\frac{1}{m}{\sum }_{i=1}^m(y^i-h_{\theta }(x^i))x_j^i$，则：
$${\theta }_j^{\prime }={\theta }_j+\frac{1}{m}\sum _{i=1}^m(y^i-h_{\theta }(x^i))x_j^i$$
这样虽然能比较容易地取得最优解，但是对于每个权重地每次更新，我们都考虑对所有样本进行计算，在权重数多或样本量大的情况下，梯度下降会非常缓慢。因此，我们考虑每次随机地选择一个样本进行更新：$${\theta }_j^{\prime }={\theta }_j+\frac{1}{batch}\sum _{i=1}^{batch}(y^i-h_{\theta }(x^i))x_j^i$$
这就是 随机梯度下降 。

但是新的问题随之而来，如果存在恶意样本，那么选择这个样本的时候可能会破坏原本的梯度下降，产生副作用。因此，我们考虑将全部样本集划分为小批量样本集，每次选择一个小批量数据集进行更新，即避免了样本量大导致的缓慢，又减少了恶意样本带来的破坏。

如果每次都这样进行更新，则样本量确定的情况下，步长也确定了，这似乎不太美妙，我们在训练的过程中更多的希望能自行控制步长，而不是完全依赖样本量，因此，我们在计算步长时引入 学习率 ，控制训练过程中的步长：$${\theta }_j^{\prime }={\theta }_j+\alpha \cdot \frac{1}{batch}\sum _{i=1}^{batch}(y^i-h_{\theta }(x^i))x_j^i$$

这就是 梯度下降更新参数 的方法。

二、算法实现

  上一阶段我们已经了解线性模型，损失函数以及梯度下降的步骤，这一阶段将使用python实现线性回归算法。

首先我们需要一些对数据集进行预处理的方法，主要作用是对数据集进行标准化和归一化，ML_algorithm包已上传百度网盘

链接: https://pan.baidu.com/s/12CsQM8CILOrHnS0rQQOYYw?pwd=6666 
提取码: 6666

1、 定义 LinearRegression 类，并在__init__函数中进行初始化，$self.theta$ 即为参数矩阵，初始化为一个值为0的列向量。

2、 在 train 方法中对线性模型进行训练，执行梯度下降更新参数矩阵。

3、 在 gradient_descent 函数中进行迭代，在 gradient_step 中进行梯度计算，即使用上一阶段得到的梯度下降更新参数的方法 更新参数矩阵 。

4、 在 loss_function 中进行损失计算，即 利用 $J(\theta )$ 计算损失 ，在计算中，由于我们使用矩阵进行计算，求得的是所有样本的损失和。因此，需要除以样本量，取平均值 。

🍁 线性回归算法完整代码：

import numpy as np
from ML_algorithm.utils.features.prepare_for_training import prepare_for_training

class LinearRegression:
    def __init__(self,data,labels,polynomial_degree = 0,sinusoid_degree = 0,normalize_data = True):
        """
        1.对数据进行预处理
        2.先得到所有的特征个数
        3.初始化参数矩阵
        """
        (data_processed,
         features_mean,
         features_deviation) = prepare_for_training(data,polynomial_degree,sinusoid_degree,normalize_data)
        
        self.data = data_processed
        self.labels = labels
        self.features_mean = features_mean
        self.features_deviation = features_deviation
        self.polynomial_degree = polynomial_degree
        self.sinusoid_degree = sinusoid_degree
        self.normalize_data = normalize_data

        num_features = self.data.shape[1] # 数据的特征个数
        self.theta = np.zeros((num_features,1)) # 初始化参数矩阵为一个初始值为0的num_features行1列的矩阵

    def train(self,alpha,num_iterations = 500):
        loss_list = self.gradient_descent(alpha,num_iterations)
        return self.theta,loss_list

    def gradient_descent(self,alpha,num_iterations):
        """
        梯度下降
        返回一个损失矩阵
        alpha : 学习率
        num_iterations : 迭代次数
        """
        loss_history = []
        for _ in range(num_iterations):
            self.gradient_step(alpha)
            loss_history.append(self.loss_function(self.data,self.labels))
        return loss_history
            
    def gradient_step(self,alpha):
        """
        梯度下降参数更新计算方法，注意是矩阵运算
        """
        num_examples = self.data.shape[0]
        prediction = LinearRegression.hypothesis(self.data,self.theta)
        delta = self.labels - prediction
        theta = self.theta
        theta = theta + alpha * (1 / num_examples) * (np.dot(delta.T,self.data)).T # 将循环累加转化为矩阵点乘，更新theta
        self.theta = theta

    def loss_function(self,data,labels):
        """
        损失计算
        """
        num_examples = data.shape[0]
        delta = labels - LinearRegression.hypothesis(self.data,self.theta)
        loss = (1 / 2) * np.dot(delta.T,delta) / num_examples # 取平均值
        return loss[0][0]

    @staticmethod
    def hypothesis(data,theta):
        """
        预测值计算
        """
        predictions = np.dot(data,theta)
        return predictions
    
    def get_loss(self,data,labels):
        data_processed = prepare_for_training(data,
                                              self.polynomial_degree,
                                              self.sinusoid_degree,
                                              self.normalize_data)[0]
        
        return self.loss_function(data_processed,labels)
    
    def get_prediction(self,data):
        """
        用训练的参数模型，预测得到回归值结果
        """
        data_processed = prepare_for_training(data,
                                              self.polynomial_degree,
                                              self.sinusoid_degree,
                                              self.normalize_data)[0]
        
        predictions = LinearRegression.hypothesis(data_processed,self.theta)

三、应用实例

  这一阶段，我们要将理论部分和实现的算法进行实战。

在ML_algorithm包中提供了数据集world-happiness-report-2017.csv，我们将利用这一数据集进行实战。

1、 首先，使用 read_csv 方法读取数据集，再用 .sample 方法对数据集进行切割，分成训练集和测试集。

2、 然后，对训练集和测试集进行分离，获得变量和结果对应的 向量。

绘制数据分布散点图：

3、 接着，导入 LinearRegression 类， 使用训练集对线性模型进行训练 。

绘制损失函数折线图：

4、 最后，对线性模型进行测试。

线性方程折线图：

可以看出：得到的线性模型能够较准确地吻合数据分布。

🍁 实战完整代码：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from linear_regression import LinearRegression

# 数据集处理部分
dataset = pd.read_csv('ML_algorithm/data/world-happiness-report-2017.csv')

train_data = dataset.sample(frac = 0.8) # 训练集
test_data = dataset.drop(train_data.index) # 测试集

input_param_name = 'Economy..GDP.per.Capita.'
output_param_name = 'Happiness.Score'

x_train = train_data[[input_param_name]].values
y_train = train_data[[output_param_name]].values

x_test = test_data[input_param_name].values
y_test = test_data[output_param_name].values

plt.scatter(x_train,y_train,label = 'Train data') # scatter绘制散点图，数据1
plt.scatter(x_test,y_test,label = 'Test data') # 数据2
plt.xlabel(input_param_name) # x轴
plt.ylabel(output_param_name) # y轴
plt.title('Happy') # 图名
plt.legend()
plt.show()

# 训练部分
num_iterations = 500 # 迭代轮次
learning_rate = 0.01 # 学习率

linear_regression = LinearRegression(x_train,y_train)
(theta,loss) = linear_regression.train(learning_rate,num_iterations)

print('开始时的损失',loss[0])
print('训练后的损失',loss[-1])

plt.plot(range(num_iterations),loss) # plot绘制折线图
plt.xlabel('Iter')
plt.ylabel('loss')
plt.title('GD')
plt.show()

# 测试部分
predictions_num = 100
x_precitions = np.linspace(x_train.min(),x_train.max(),predictions_num).reshape(predictions_num,1)
y_precitions = linear_regression.get_prediction(x_precitions)

plt.scatter(x_train,y_train,label = 'Train data') # scatter绘制散点图，数据1
plt.scatter(x_test,y_test,label = 'Test data') # 数据2
plt.plot(x_precitions,y_precitions,'r',label = 'Prediction')
plt.xlabel(input_param_name) # x轴
plt.ylabel(output_param_name) # y轴
plt.title('Happy') # 图名
plt.legend()
plt.show()