线性回归模型

线性回归

一、理论部分

线性回归：在有限的数据中，通过调整参数，来拟合出一条直线，并以此直线（模型）来对未知的数据进行预测。

直线的一般形式：
$$y=w\times x+b$$
其中整个直线的状态由$w和b$决定，$w$决定直线的斜率（也就是倾斜角度），$b$决定直线在Y轴上的截距（控制直线的上下平移，也称为偏置）。因此，我们只需要计算出$w$和$b$的值即可确定出一条具体的直线，接下来研究如何找到$w$和$b$。

w和b对直线形状的影响

专有名词

名词	释义
训练集	用于训练模型的数据集
测试集	用于检验模型的质量的数据集
训练样本/训练点（Training Sample）	训练集中的每一项数据
特征（Features）	输入到模型中的数据，可以是数值、类别标签、图像像素值等。
标签（Label）	模型需要预测的目标值。对于分类问题，标签通常是类别名称；对于回归问题，标签则是一个连续值。

1. 单变量线性回归

1）数据预处理

我们先研究单变量线性回归，所谓单变量线性回归是指只有一个自变量的一次函数，如：$y=w\cdot x+b$即为单变量直线，只具有一个输入变量$x$。该直线可以表示在二维平面上（横轴为X，纵轴为Y）。

在拿到一组未划分的数据时，我们通常会将数据分为训练集与测试集，一个简单的划分方式是：取样本的前80%作为训练集，剩余20%作为测试集。

2）定义代价函数

假定我们已经求出来了$w$和$b$，那么我们就确定了一条直线，就可以使用该直线进行预测了，为了方便判断我们所预测出来的值y’与真实的值y之间的误差是多少，我们要定义*“一把尺子”*，用于衡量预测值$y^{\prime }$与真实值$y$之间的误差。这里我们使用均方误差来定义代价函数：

$$J(w,b)=\frac{1}{2m}\sum _{i=1}^m(f_{w,b}(x^{(i)})-y^{(i)}{)}^2$$
公式拆解：

$f_{w,b}(x^{(i)})-y^{(i)}$：其中$f_{w,b}(x^{(i)})$代表了使用训练好的模型预测得到的值，而$y^{(i)}$代表了每一个训练样本的真实结果，这两者的差值表示了模型预测得到的值与真实的值之间的误差。

为什么要将误差进行平方？

在所有样本集中，每一个误差可能为正数也可能为负数，进行求和的过程中会有一定概率抵消，这样做会导致在每一项的误差都很大时（如：-100，+90，-25，+30），进行求和之后得到一个小的值（-5），最终得出一个错误的判断。

$\frac{1}{2m}$：代表了对所有数据误差之和取平均值（这个平均值从某种意义上来讲可以代表整个模型的误差），得到均方误差。

为什么要除以2？

因为在后面进行梯度下降的时候，求导会将指数2给划分到系数上，因为对于大量的数据而言，一个常数对模型的影响不大，为了求导公式的简介，这里提前除以2，以便之后的抵消。

知道了代价函数，只需要想办法将代价降低即可，代价越低，说明我们预测的值与真实的值越接近。

通过观察误差成本函数，可以看到误差成本函数是一个二次函数，即是一个凸函数，凸函数的一个性质是：极值点就是最值点，由于代价函数是一个开口向上的二次函数（可以将公式进行展开，能够直观的感受到平方项的系数是大于0的），因此该凸函数只具有最小值，我们只需要找到代价函数的极小值就得到了最小值。对于误差成本函数$J(w,b)$ , 其公式展开可以写为：
$$J(w,b)=\frac{1}{2m}\sum _{i=1}^m((w{x}^{(i)}+b)-y^{(i)}{)}^2$$
$J$的大小取决于参数$w$和$b$，这可以通过梯度下降的方法进行求解，代价函数的形状大致如下：

3）梯度下降

梯度下降的思想，主要运行了求偏导的方法，这与生物学上的控制变量的方法很相似，如：在控制$b$不变的情况下去更新$w$（可视$b$为常数），公式：$w^{\prime }=w-\alpha \frac{\partial J(w)}{\partial w}$ 表示依次更新了$w$，其中$\alpha$ 代表了学习率 用于表示步长，也可以理解为下降的速度，$\frac{\partial J(w)}{\partial w}$ 表示对$w$求偏导，得到的是 $W-J$(权重W与代价J的函数) 凸函数上的一个切线，用于表示函数值下降最快的方向，这两者的乘积表示向函数值下降最快的方向进行移动一个步长。这个距离需要依据数据集进行调整，如果$\alpha$过大（步长过大），会导致$w$直接越过最低点而到另一端的高点，以至于永远无法接近最小值，如果$\alpha$过小（步长过小），会导致$w$接近最小值时越来越慢，消耗计算成本。

---

学习率（$\alpha$）的调整方法：

1. 先设置一个较小的$\alpha$ 如：0.001 。

2. 然后每次增加10倍，最大为1。

3. 确定某一值后，如：0.01。

4. 再进行3倍的处理，如：$0.01\times 3=0.03,0.03\times 3=0.09$ （这样做的目的是使收敛的速度更快）。

---

求解偏导数的过程（寻找梯度下降的方向）：

• 求$\frac{\partial J(w)}{\partial w}$ :

  $\frac{\partial J(w)}{\partial w}$ = $\frac{\partial }{\partial w}\frac{1}{2m}\sum _{i=1}^m(f_{w,b}(x^{(i)})-y^{(i)}{)}^2$

  ​ = $\frac{\partial }{\partial w}\frac{1}{2m}\sum _{i=1}^{m}w{x}^{(i)}-y^{(i)}{)}^2$

  ​ =$\frac{1}{2m}\sum _{i=1}^m(f_{w,b}(x^{(i)})-y^{(i)})\cdot 2x^{(i)}$

  ​ =$\frac{1}{m}\sum _{i=1}^m(f_{w,b}(x^{(i)})-y^{(i)})x^{(i)}$

• 求$\frac{\partial J(w)}{\partial b}$ :

  $\frac{\partial J(w)}{\partial b}$ = $\frac{\partial }{\partial b}\frac{1}{2m}\sum _{i=1}^m(f_{w,b}(x^{(i)})-y^{(i)}{)}^2$

  ​ = $\frac{\partial }{\partial b}\frac{1}{2m}\sum _{i=1}^{m}w{x}^{(i)}-y^{(i)}{)}^2$

  = $\frac{1}{2m}\sum _{i=1}^m(f_{w,b}(x^{(i)})-y^{(i)})\cdot 2$

  ​ =$\frac{1}{m}\sum _{i=1}^m(f_{w,b}(x^{(i)})-y^{(i)})$

通过循环迭代的方法找到具体的$w$值与$b$值：

$while():$

​ $w^{\ast }=w-\alpha \frac{\partial }{\partial w}J(w,b)$

​ $b^{\ast }=b-\alpha \frac{\partial }{\partial b}J(w,b)$

​ $w=w^{\ast }$

​ $b=b^{\ast }$

初始的时候我们随机的对w和b进行取值，然后进行迭代，可以设定当误差小于某个阈值的时候退出梯度下降，或者自定义迭代的次数，通过有限步骤可以得到一个较为理想的$w$值和$b$值。

2.多变量线性回归

多变量线性回归，将维度扩充到了三维以至多维状态，如$y=w_1{x}_1+w_2{x}_2+b$ 可以理解为X轴和Y轴分别为$x_1和x_2$ Z轴为$y$，这样就是一个三维的立体状态，每一个训练样本都是三维空间的一个点，要找到一个合适直线拟合该三维空间的样本点。

做法：分别对每一个变量（$w_1,w_2,\dots ,w_n,b$）进行梯度下降处理

要点：对于多变量线性回归，其中不同的特征值取值范围不同，如：年龄特征取值范围：$0$ ~ $100$，占地面积：$0m^2$ ~ $10000m^2$ ，还可能会有奇异样本，奇异样本数据的存在会引起训练时间增大，同时也可能导致无法收敛，因此，当存在奇异样本数据时，在进行训练之前需要对预处理数据进行归一化；反之，不存在奇异样本数据时，则可以不进行归一化。针对此问题，我们需要对数据进行特征进行缩放（归一化处理）。

1）数据预处理

数据归一化处理

实现数据归一化的方法有三种方法：

1. 同除以最大值：

   对每个特征中的所有值都除以该特征中的最大值。

2. 均值归一化：

   对每个特征中的值减去该特征的均值，然后再除以该特征最大值与最小值的之间的差值。

3. Z-score标准化：

   计算每个特征的标准差和均值

   对每个特征的所有值减去该特征的平均值，然后再除以该特征的均值

如果不进行归一化，那么由于特征向量中不同特征的取值相差较大，会导致代价函数变 “扁”。这样在进行梯度下降的时候，梯度的方向就会偏离最小值的方向，走很多弯路，会使训练时间过长。

进行归一化以后，目标函数会呈现比较 “圆”，这样训练速度大大加快，少走很多弯路。

数据归一化的好处：

1. 归一化后加快了梯度下降求最优解的速度，也即加快训练网络的收敛性。
2. 归一化有可能提高精度。

2）代价函数（与单变量线性回归的代价函数相同）

$$J(w,b)=\frac{1}{2m}\sum _{i=1}^m(f_{w,b}(x^{(i)})-y^{(i)}{)}^2$$

3）梯度下降

对于多元线性回归的梯度下降：

$w_1=W_1-\alpha \frac{1}{m}\sum _{i=1}^m(f_{\vec{w},b}{\vec{x}}^{(i)}-y^{(i)})x_1^{(i)}$

​ $⋮$

$w_n=W_n-\alpha \frac{1}{m}\sum _{i=1}^m(f_{\vec{w},b}{\vec{x}}^{(i)}-y^{(i)})x_n^{(i)}$

$b=b-\alpha \frac{1}{m}\sum _{i=1}^m(f_{\vec{w},b}{\vec{x}}^{(i)}-y^{(i)})$

解释：$w_1\cdots w_n$ 表示每一个变量的系数，b表示线性函数的常数项。

3.正规方程

1）数据预处理

省略……

2）代价函数

数学推导：

$J=\frac{1}{2m}\sum _{i=1}^m({\vec{\theta }}_i{\vec{x}}_i-y_i{)}^2$

​ =$\frac{1}{2m}||\vec{\theta }\vec{x}-y|{|}^2$

​ =$\frac{1}{2m}(\vec{\theta }\vec{x}-y{)}^T(\vec{\theta }\vec{x}-y)$

​ =$\frac{1}{2m}({\vec{\theta }}^T{\vec{x}}^T-y^T)(\vec{\theta }\vec{x}-y)$

​ =$\frac{1}{2m}({\vec{\theta }}^T{\vec{x}}^T\vec{x}\vec{\theta }-y^T\vec{x}\vec{\theta }-{\vec{\theta }}^T{\vec{x}}^{T}y+y^{T}y)$

3）梯度下降

对$\theta$求偏导得： $\Delta =\frac{\partial J}{\partial \theta }=\frac{1}{2m}(\frac{\partial {\vec{\theta }}^T{\vec{x}}^T\vec{x}\vec{\theta }}{\partial \theta }-\frac{\partial y^T\vec{x}\vec{\theta }}{\partial \theta }-\frac{\partial {\vec{\theta }}^T{\vec{x}}^{T}y}{\partial \theta }+\frac{y^{T}y}{\partial \theta })$

矩阵求导法则：

1. $\frac{\partial {\theta }^{T}A\theta }{\partial \theta }=(A+A^T)\theta$

2. $\frac{\partial X^{T}A}{\partial X}=A$

3. $\frac{\partial AX}{\partial X}=A^T$

4. $\frac{\partial A}{\partial X}=0$

可得 $\Delta =\frac{1}{2m}(\frac{\partial {\vec{\theta }}^T{\vec{x}}^T\vec{x}\vec{\theta }}{\partial \theta }-\frac{\partial y^T\vec{x}\vec{\theta }}{\partial \theta }-\frac{\partial {\vec{\theta }}^T{\vec{x}}^{T}y}{\partial \theta }+\frac{y^{T}y}{\partial \theta })=\frac{1}{2m}\cdot (2x^{T}x\theta -2x^{T}y)=\frac{1}{m}\cdot (x^{T}x\theta -x^{T}y)$
当$\Delta =0$ 时: $x^{T}x\theta =x^{T}y$ ,$\theta =(x^{T}x{)}^{-1}{x}^{T}y$ 可以计算出 $\theta$ 的值。

梯度下降与正规方程的比较：

• 梯度下降：需要选择学习率α，需要多次迭代，当特征数量n大时也能较好适用，适用于各种类型的模型

• 正规方程：不需要选择学习率α，一次计算得出，需要计算 $(x^{T}x{)}^{-1}$ ，如果特征数量$n$较大则运算代价大，因为矩阵逆的计算时间复杂度为 $O(n^3)$ ，通常来说当 $n$小于10000 时还是可以接受的，只适用于线性模型，不适合逻辑回归模型等其它模型。

4.多项式回归

由于在某些情况下，直线难以拟合全部的数据，需要曲线来对数据进行拟合，如：二次模型、三次模型等等。

一般情况下，数据的回归函数是未知的，即使已知也难以用一个简单的函数变换转化为线性模型，所以常用的做法是多项式回归(Polynomial Regression)，即使用多项式函数拟合数据。

如何选择多项式的次数：

多项式函数有多种，一般来说，需要先观察数据的形状，再去决定选用什么形式的多项式函数来处理问题。比如，从数据的散点图观察，如果有一个“弯 ”，就可以考虑用二次多项式（即对特征进行平方）；有两个“弯 ”，可以考虑用三次多项式（对特征取三次方）；有三个“弯 ”，则考虑用四次多项式（对特征取四次方），以此类推。

虽然真实的回归函数不一定是某个次数的多项式，但只要拟合的好，用适当的多项式来近似模拟真实的回归函数是可行的。

二、实验部分

文章最后附录部分有实验用到的所有原数据，ex1data1.txt是人口与利润之间的关系， ex1data2.txt是房子大小与卧室数量对房子价格的影响。

1. 单变量线性回归

1）加载数据

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
path = "ex1data1.txt"
data = pd.read_csv(path,header = None,names=['Population','Profit'])
data.head()   # 预览数据

2）查看数据

data.describe()    # 更加详细的数据描述

# 可视化训练数据
data.plot(kind = 'scatter',x = 'Population',y = 'Profit',figsize = (12,8))
plt.show()

3）定义代价函数

def computerCost(X,y,theta):    # 定义代价函数
    inner = np.power(((X*theta.T)-y),2)   # theta.T表示theta的转置
    return np.sum(inner)/(2*len(X))

data.insert(0,"One",1)  # 表示在第0列前面插入一列数，其表头为One，其值为1

在数据集第一列插入$1$的作用是为了便于矩阵的计算，在矩阵相乘的时候，涉及到权重$w$和偏置$b$，由于$b$没有与变量进行相乘，因此在变量前面插入一个$1$，用于与$b$进行相乘。

4）拆分数据

cols = data.shape[1]
X = data.iloc[:,0:cols - 1]  #“，”前只有“：”，表示所有的行，“，”后表示抽取数据中第[0列~第cols-1列)(左闭右开)，去掉最后一列，最后一列为预测值

y = data.iloc[:,cols - 1:cols]  #只取最后一列的值，表示预测值

X.head()

y.head()

X = np.matrix(X.values)
y = np.matrix(y.values)  #只将表格中的值装换为矩阵而不是包括序号与标题

#初始化theta
theta = np.matrix(np.array([0,0]))  #先是一个一维的数据，然后在转换为一个二维的矩阵

5）初始化参数

theta  
# => matrix([[0, 0]])

X.shape,theta.shape,y.shape  # 此时theta为一行列，需要进行转置
# => ((97, 2), (1, 2), (97, 1))

computerCost(X,y,theta)
# => 32.072733877455676

6）定义梯度下降函数

def gradientDecent(X,y,theta,alpha,iters):   #iters为迭代次数
    temp = np.matrix(np.zeros(theta.shape))   #构造一个与theta大小一样的零矩阵，用于存储更新后的theta
    parmaters = int (theta.ravel().shape[1])    #.ravel()的功能是将多维数组降至一维，用于求需要求的参数个数
    cost = np.zeros(iters)   #构建iters个0的数组，相当于对每次迭代的cost进行记录
    
    for i in range(iters):
        error = (X * theta.T - y)     #记录误差值，结果为一个数组
        for j in range(parmaters):    #对每一个参数进行更新，j用于表示每一个参数
            term = np.multiply(error,X[:,j])   #.multiply 是对矩阵当中的数对应相乘，这里表示与X矩阵的第j列相乘。
            temp[0,j] = theta[0,j] - ((alpha/len(X))*np.sum(term))  #存储更行后的theta的值，.sum()表示将矩阵中的数进行求和
        
        theta = temp      #更新theta
        cost[i] = computerCost(X,y,theta)  #计算此时的代价，并记录在cost中。
    return theta,cost

7）初始化超参数

alpha  = 0.01		# 学习率
iters = 1000		# 迭代次数

8）梯度下降

g,cost = gradientDecent(X,y,theta,alpha,iters)
g
# => matrix([[-3.24140214,  1.1272942 ]])

9）计算代价

computerCost(X,y,g)
# => 4.515955503078914

10）绘制线性模型

x = np.linspace(data.Population.min(),data.Population.max(),100) #抽取100个样本  (从data数据集中的最小值到最大值之间抽取100个样本)
f = g[0,0] + (g[0,1] * x)  #f = ax + b

fig,ax = plt.subplots(figsize = (12,8))    #figsize表示图的大小
ax.plot(x,f,'r',label = "Prediction")    #绘制直线，横坐标，纵坐标，直线名称
ax.scatter(data.Population,data.Profit,label = 'Training data')   #绘制点，横坐标，纵坐标，点的名称
ax.legend(loc = 4)  #显示图例位置
ax.set_xlabel('Population')  #设置x轴的名称
ax.set_ylabel('Profit')   #设置y轴的名称
ax.set_title('Predicted Profit vs. Population Size')  #设置标题的名称
plt.show()

11）绘制代价变化曲线

fig,ax = plt.subplots(figsize = (12,8))
ax.plot(np.arange(iters),cost,'r')
ax.set_xlabel('Interations')
ax.set_ylabel('Cost')
ax.set_title("Error vs. Training Epoc")
plt.show()

2. 多变量线性回归

1）加载数据

path = "ex1data2.txt"
data2 = pd.read_csv(path,header = None,names=["Size","Bedroom","Price"])
data2.head()

2）数据归一化处理

data2 = (data2 - data2.mean())/data2.std()
data2.head()

3）拆分数据

data2.insert(0,'Ones',1)  #在x的第一列插入1

clos = data2.shape[1]   #存储第二维（列）的数据量
X2 = data2.iloc[:,0:clos-1]  #对X2进行赋值
y2 = data2.iloc[:,clos-1:clos]  #对y2进行赋值

X2 = np.matrix(X2.values)  #将X2转为矩阵
y2 = np.matrix(y2.values)  #将y2转为矩阵
theta2 = np.matrix(np.array([0,0,0]))  #初始化theta2为0矩阵
computerCost(X2, y2, theta2)
# => 0.48936170212765967

4）梯度下降

g2,cost2 = gradientDecent(X2,y2,theta2,alpha,iters)   #记录放回值g2（theta2）和cost2
g2
# => matrix([[-1.10868761e-16,  8.78503652e-01, -4.69166570e-02]])

5）计算代价

computerCost(X2,y2,g2)
# => 0.13070336960771892

6）绘制代价变化曲线

fig,ax = plt.subplots(figsize = (12,8))
ax.plot(np.arange(iters),cost2,'x')
ax.set_xlabel('Iterations')
ax.set_ylabel('Cost')
ax.set_title('Error vs. Training Epoch')
plt.show()

3. 正规方程

#正规方程
def normalEqn(X,y):
    theta = np.linalg.inv(X.T@X)@X.T@y   #.linalg中包含线性代数中的函数，求矩阵的逆（inv）、特征值等。@表示矩阵相乘
    return theta

final_theta2 = normalEqn(X,y)
final_theta2
# => matrix([[-3.89578088], [ 1.19303364]])

三、总结

训练模型的一般步骤：

1. 数据预处理。
2. 根据具体问题选择模型。
3. 设置代价函数。
4. 利用梯度下降算法寻找最优的参数。
5. 进行模型的评估，调整超参数。
6. 使用模型进行预测。

四、附录

1. ex1data1.txt

6.1101,17.592
5.5277,9.1302
8.5186,13.662
7.0032,11.854
5.8598,6.8233
8.3829,11.886
7.4764,4.3483
8.5781,12
6.4862,6.5987
5.0546,3.8166
5.7107,3.2522
14.164,15.505
5.734,3.1551
8.4084,7.2258
5.6407,0.71618
5.3794,3.5129
6.3654,5.3048
5.1301,0.56077
6.4296,3.6518
7.0708,5.3893
6.1891,3.1386
20.27,21.767
5.4901,4.263
6.3261,5.1875
5.5649,3.0825
18.945,22.638
12.828,13.501
10.957,7.0467
13.176,14.692
22.203,24.147
5.2524,-1.22
6.5894,5.9966
9.2482,12.134
5.8918,1.8495
8.2111,6.5426
7.9334,4.5623
8.0959,4.1164
5.6063,3.3928
12.836,10.117
6.3534,5.4974
5.4069,0.55657
6.8825,3.9115
11.708,5.3854
5.7737,2.4406
7.8247,6.7318
7.0931,1.0463
5.0702,5.1337
5.8014,1.844
11.7,8.0043
5.5416,1.0179
7.5402,6.7504
5.3077,1.8396
7.4239,4.2885
7.6031,4.9981
6.3328,1.4233
6.3589,-1.4211
6.2742,2.4756
5.6397,4.6042
9.3102,3.9624
9.4536,5.4141
8.8254,5.1694
5.1793,-0.74279
21.279,17.929
14.908,12.054
18.959,17.054
7.2182,4.8852
8.2951,5.7442
10.236,7.7754
5.4994,1.0173
20.341,20.992
10.136,6.6799
7.3345,4.0259
6.0062,1.2784
7.2259,3.3411
5.0269,-2.6807
6.5479,0.29678
7.5386,3.8845
5.0365,5.7014
10.274,6.7526
5.1077,2.0576
5.7292,0.47953
5.1884,0.20421
6.3557,0.67861
9.7687,7.5435
6.5159,5.3436
8.5172,4.2415
9.1802,6.7981
6.002,0.92695
5.5204,0.152
5.0594,2.8214
5.7077,1.8451
7.6366,4.2959
5.8707,7.2029
5.3054,1.9869
8.2934,0.14454
13.394,9.0551
5.4369,0.61705

2. ex1data2.txt

2104,3,399900
1600,3,329900
2400,3,369000
1416,2,232000
3000,4,539900
1985,4,299900
1534,3,314900
1427,3,198999
1380,3,212000
1494,3,242500
1940,4,239999
2000,3,347000
1890,3,329999
4478,5,699900
1268,3,259900
2300,4,449900
1320,2,299900
1236,3,199900
2609,4,499998
3031,4,599000
1767,3,252900
1888,2,255000
1604,3,242900
1962,4,259900
3890,3,573900
1100,3,249900
1458,3,464500
2526,3,469000
2200,3,475000
2637,3,299900
1839,2,349900
1000,1,169900
2040,4,314900
3137,3,579900
1811,4,285900
1437,3,249900
1239,3,229900
2132,4,345000
4215,4,549000
2162,4,287000
1664,2,368500
2238,3,329900
2567,4,314000
1200,3,299000
852,2,179900
1852,4,299900
1203,3,239500