机器学习(2-4)：线性回归之四：多元线性回归

2.4  多元线性回归

2.4.1 假设函数、样本数量、代价函数

假设函数

代数表达

比如对于多元线性回归，假设函数表示为

其中 为模型参数， 为每个样本的 个特征值。这个表示可以简化，我们增加一个特征 ，这样

矩阵表达

设向量 ， ，则有：

样本

代数表达

假设样本个数为m，具体为 。

矩阵表达

设输入特征向量 ，输出特征向量

代价函数

代数表达

矩阵表达

设矩阵 ，则有

2.4.2 多元梯度下降法

1、假设函数：  或

2、参数： 或

3、样本：样本数为m，具体为 或

输入特征向量 ，输出特征向量

4、代价函数： =12m(Xθ-y)T(Xθ-y)

5、梯度下降：

特征缩放(Fearture Scaling)——归一化

由于样本不同特征的取值范围不一样，可能导致迭代很慢，为了减少特征取值的影响，可以对特征数据归一化。

一般地而言，将每个输入的特征值限制在,n ，这就是归一化。实际中不是绝对的，经验告诉我们，当 或 ，就不需要归一化，否则就需要归一化。

均值归一化

就是对于每个特征 ，求出它的期望 和标准差 ，然后转化为：

这样特征的新期望为0，新方差为1，迭代速度可以大大加快。

几个概念

总体标准差：

样本标准差： 

标准误差： 

方差：方差=标准差的平方。

学习速率(Learning rate)——步长

如何选择合适的α ，使梯度下降的迭代步数更少。

步长太大，会导致迭代过快，甚至有可能错过最优解。步长太小，迭代速度太慢，很长时间算法都不能结束。所以算法的步长需要多次运行后才能得到一个较为优的值。

在前面的算法描述中，提到取步长为1，但是实际上取值取决于数据样本，可以多取一些值，从大到小，分别运行算法，看看迭代效果，如果损失函数在变小，说明取值有效，否则要增大步长。

特征和多项式回归

根据训练集的m 个样本特征数n ，我们可以得出假设函数为

我们可以创建新特征来取得更好的效果。

创造新特征降低特征数

对于一块土地的价值，

其中 表示临界宽度， 表示深度。

我们可以创建一个新特征“面积=临街宽度x深度”，则假设函数就可以变为：

其中 表示土地面积。

多项式回归

上面所提到的都是线性回归（一元，多元），当线性回归不能准确拟合样本曲线时，我们可以多项式回归。

创造新特征使多项式回归变成线性回归

对于多项式回归

在其中，我们有可以创建新特征，令 ，则多项式回归就变成了多元线性回归

创造新特征规避多项式回归中偶函数的问题

对于多项式回归

若上面的假设函数中 是房屋面积，由于 是偶函数，有可能会造成面积越大，总价越小的尴尬。我们可以对假设函数进行改造为：

也可以改造为

创建新特征 ，则多项式回归就变成了多元线性回归：

2.4.3 正规方程法（区别于迭代方法的直接解法）

之前我们都是用梯度下降法来求代价函数的最小值。这里我们用方程直接求解最小值。

通过高等数学我们知道：

极值点出现在偏导数=0的地方，求解导数方程或偏导数方程组，即可求解极值点坐标。

对于一维 ，若 ， ，则得

对于n+1维 ，

求解上述n+1个方程组成的方程组，则可解出 。

步骤详解

房屋价格，样本数m=4，具体情况见下表：

面积(feet2)	卧室数量	层数	房龄(years)	价格（$1000）
2104	5	1	45	460
1416	3	2	40	232
1534	3	2	30	315
852	2	1	36	178

第一步，在数据集中增加一列 ，其取值永远为1

1	2104	5	1	45	460
1	1416	3	2	40	232
1	1534	3	2	30	315
1	852	2	1	36	178

第二步，构建一个输入特征矩阵X和输出向量y

第三步，求解θ

梯度下降法和正规方程法的优缺点

梯度下降法需要选择学习速率，需要多次迭代；但可以工作在任意多的输入特征量n，尤其是n很大时。

正规方程法不需要选择学习速率，不需要迭代；但需要计算 ，当n很大时，计算会非常缓慢。

如何选择呢？

对于线性回归，当n为10000以下，用正规方程法，否则用梯度下降法。

当算法越来越复杂时，正规方程法就不适用了，此时选择梯度下降法。

正规方程在矩阵不可逆情况下的解决办法

不可逆怎么办？

矩阵不可逆有两种情况：

1. 一是特征值线性相关
2. 二是样本数量m小于特质值数量n。

解决办法有：

1. 针对特征值线性相关，去掉线性相关的特征值
2. 针对样本数量m小于特质值数量n，要么删除一些特征值，要么正则化（regularization），在后面会讲正则化。

2.4.5 示例：【Python】多元线性回归（梯度下降法）

一个房屋价格数据集，其中有2个变量（房子的大小，卧室的数量）和目标（房子的价格）。 我们使用我们已经应用的技术来分析数据集。

样本数据data2.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

Python代码（梯度下降法）

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d

'''
  样本输入特征数为N (N=1 对于一元线性回归)，输出特征为y。
  假设函数：y=theta0+theta1*x1+theta2*x2+...+thetaN*xN
      令x0 =1,
        x = [x0,x1,x2,...,xN]
        theta = [theta0,theta1,theta2,...thetaN]
      则有：
        y = x * theta.T
  样本个数为M, 样本的输入特征矩阵X为 (M行N+1列)，输出特征y为向量(M行1列)
  代价函数：cost = sum((X*theta.T-y)^2)/(2*M)
  
'''

def computeCost(X, y, theta):
    '''
    作用：计算代价函数,向量化来计算参数
    :param X: 样本的输入特征，矩阵(M行N+1列)。
    :param y: 样本的输出特征，向量(M行1列)
    :param theta: parameters，一元线性回归的参数，向量(1行N+1列
    :return: 代价函数的值，浮点数
    '''
    inner = np.power(((X * theta.T) - y), 2)
    # print(inner)
    return np.sum(inner) / (2 * len(X))

def gradientDescent(X, y, theta, alpha, epoch):
    '''
    作用:获得最终梯度下降后的theta值以及cost
    :param X: 输入特性矩阵(M,N+1),二元线性回归N=2，矩阵形状(M,3)
    :param y: 输出特性矩阵(1,M)
    :param theta: 假设函数参数矩阵(1,N+1)。二元线性回归N=2，矩阵形状(1,3)
    :param alpha: 步长，又称学习速率,浮点数。
    :param epoch: 要执行的迭代次数，整数。
    :return:
    '''
    # 1 变量初始化，储存数据
    # 1.1 初始化参数矩阵为0
    np.matrix(np.zeros(theta.shape)) # 初始化一个临时矩阵(1, 3)
    # flatten()降维 即返回一个折叠成一维的数组。
    # 但是该函数只能适用于numpy对象，即array或者mat，普通的list列表是不行的。
    parameters = int(theta.flatten().shape[1]) # 参数theta的数量 3
    # print(parameters)
    # 1.2 初始化保存每一次迭代计算代价函数值的数组为0
    cost = np.zeros(epoch)  # 初始化一个ndarray, 包含每次训练后的cost #1000个0的矩阵
    # print(cost)
    # 1.3 初始化保存每一次迭代计算假设函数参数的数组(epoch,3)为0
    counterTheta = np.zeros((epoch, 3)) #1000 * 3的数组
    m = X.shape[0]  # 样本个数

    # 2 迭代计算
    for i in range(epoch):
        '''
        使用 vectorization同时更新所有的θ，可以大大提高效率,此处都是相对应的进行计算
        X.shape, theta.shape, y.shape, X.shape[0]
        ((97, 2), (1, 2), (97, 1), 97)
        '''
        # 相当于theta1 theta2不停做偏导并且更新 theta[theta1, theta2] temp是临时的theta
        temp = theta - (alpha / m) * (X * theta.T - y).T * X
        theta = temp
        counterTheta[i] = theta
        cost[i] = computeCost(X, y, theta)
        pass
    return counterTheta, theta, cost

'''
多元线性回归
'''

# 1.Prepare datasets
path = 'D:/dataAnalysis/MachineLearning/data2.txt'
# names添加列名，header用指定的行来作为标题，若原无标题且指定标题则设为None
data = pd.read_csv(path, header=None, names=['Size', 'Bedrooms', 'Price'])
print(data.head())#显示前五行
print(data.describe())

# 2.数据归一化
# 数据归一化后，最优解的寻优过程明显会变得平缓，更容易正确的收敛到最优解。简而言之，归一化的目的就是使得预处理的数据被限定在一定的范围内（比如[0,1]或者[-1,1]），从而消除奇异样本数据导致的不良影响。
temp = data  #复制一份
data = (data - data.mean()) / data.std() # data2.std()是标准差
print(data.head())#显示前五行
print(data.describe())

# 3 在数据集前增加1列，值为全1
data.insert(0, 'Ones', 1)       # 增加一条第一列，全部数值为1
# print(data)

# 4 从表中提取样本输入特征，输出特征
#   变量初始化:set X (training data) and y (target variable)
cols = data.shape[1]  # 列数
X = data.iloc[:, 0:cols - 1]        # 取前cols-1列，即输入向量
y = data.iloc[:, cols - 1:cols]     # 取最后一列作为目标向量
# print(X.head()) # 观察下 X (训练集) and y (目标变量)是否正确.
# print(y.head())

# 5 将样本输入特征，输出特征转化为矩阵
X = np.matrix(X.values)
y = np.matrix(y.values)
theta = np.matrix([0,0,0]) # theta 是一个(1,3)矩阵
computeCost(X, y, theta)

# 6 定义学习速率，迭代次数
alpha = 0.01 # 学习速率α
epoch = 3800 # 要执行的迭代次数

# 7 梯度下将法求解
counterTheta, final_theta, cost = gradientDescent(X, y, theta, alpha, epoch)
computeCost(X, y, final_theta)
print(computeCost(X, y, final_theta)) # 4.515955503078912
print(final_theta)

# 8 将归一化化后的数据还原
#data =  temp 

# 9 预测值
x1 = np.linspace(data.Size.min(), data.Size.max(), 100)  # xlabel start:返回样本数据开始点 stop:返回样本数据结束点 num:生成的样本数据量，默认为50
x2 = np.linspace(data.Bedrooms.min(), data.Bedrooms.max(), 100),# ylabel start:返回样本数据开始点 stop:返回样本数据结束点 num:生成的样本数据量，默认为50
x1 = np.matrix(x1) # 数组转换为矩阵
x2 = np.matrix(x2) # 数组转换为矩阵
f = final_theta[0, 0] + final_theta[0, 1] * x1.T + final_theta[0, 2] * x2.T  # zlabel price
# 矩阵转换为数组
#x1= x1.getA()
#x2= x2.getA()
#f= f.getA()
print(f)
print(f.shape)
print(final_theta)
print(final_theta.shape)

# 10 可视化
# 10.1 绘制样本离散图和预测曲线
fig1 = plt.subplots(figsize=(6, 4)) # 尺寸 
ax = plt.axes(projection='3d')
ax.scatter3D(data.Size, data.Bedrooms, data.Price, color="green") # 样本的离散值
ax.plot_surface(x1, x2, f,cmap='viridis', edgecolor='none') # 预测曲面：横坐标 纵坐标 颜色 标签
ax.legend(loc=2) # 2表示在左上角
ax.set_xlabel('Size')
ax.set_ylabel('Bedrooms')
ax.set_zlabel('Price')
ax.set_title('Room Price  vs. Room Size & Bedrooms')

# 10.2 绘制代价函数曲线
fig2, ax = plt.subplots(figsize=(8, 4))
ax.plot(np.arange(epoch), cost, 'r') # 横坐标 纵坐标 颜色
ax.set_xlabel('Iteration')
ax.set_ylabel('Cost')
ax.set_title('Error vs. Training Epoch')
plt.show()

2.4.6 示例：【Python】多元线性回归（正规方程法）

样本数据Adertisting.xlsx

TV	radio	newspaper	sales
230.1	37.8	69.2	22.1
44.5	39.3	45.1	10.4
17.2	45.9	69.3	9.3
151.5	41.3	58.5	18.5
180.8	10.8	58.4	12.9
8.7	48.9	75	7.2
57.5	32.8	23.5	11.8
120.2	19.6	11.6	13.2
8.6	2.1	1	4.8
199.8	2.6	21.2	10.6
66.1	5.8	24.2	8.6
214.7	24	4	17.4
23.8	35.1	65.9	9.2
97.5	7.6	7.2	9.7
204.1	32.9	46	19
195.4	47.7	52.9	22.4
67.8	36.6	114	12.5
281.4	39.6	55.8	24.4
69.2	20.5	18.3	11.3
147.3	23.9	19.1	14.6
218.4	27.7	53.4	18
237.4	5.1	23.5	12.5
13.2	15.9	49.6	5.6
228.3	16.9	26.2	15.5
62.3	12.6	18.3	9.7
262.9	3.5	19.5	12
142.9	29.3	12.6	15
240.1	16.7	22.9	15.9
248.8	27.1	22.9	18.9
70.6	16	40.8	10.5
292.9	28.3	43.2	21.4
112.9	17.4	38.6	11.9
97.2	1.5	30	9.6
265.6	20	0.3	17.4
95.7	1.4	7.4	9.5
290.7	4.1	8.5	12.8
266.9	43.8	5	25.4
74.7	49.4	45.7	14.7
43.1	26.7	35.1	10.1
228	37.7	32	21.5
202.5	22.3	31.6	16.6
177	33.4	38.7	17.1
293.6	27.7	1.8	20.7
206.9	8.4	26.4	12.9
25.1	25.7	43.3	8.5
175.1	22.5	31.5	14.9
89.7	9.9	35.7	10.6
239.9	41.5	18.5	23.2
227.2	15.8	49.9	14.8
66.9	11.7	36.8	9.7
199.8	3.1	34.6	11.4
100.4	9.6	3.6	10.7
216.4	41.7	39.6	22.6
182.6	46.2	58.7	21.2
262.7	28.8	15.9	20.2
198.9	49.4	60	23.7
7.3	28.1	41.4	5.5
136.2	19.2	16.6	13.2
210.8	49.6	37.7	23.8
210.7	29.5	9.3	18.4
53.5	2	21.4	8.1
261.3	42.7	54.7	24.2
239.3	15.5	27.3	15.7
102.7	29.6	8.4	14
131.1	42.8	28.9	18
69	9.3	0.9	9.3
31.5	24.6	2.2	9.5
139.3	14.5	10.2	13.4
237.4	27.5	11	18.9
216.8	43.9	27.2	22.3
199.1	30.6	38.7	18.3
109.8	14.3	31.7	12.4
26.8	33	19.3	8.8
129.4	5.7	31.3	11
213.4	24.6	13.1	17
16.9	43.7	89.4	8.7
27.5	1.6	20.7	6.9
120.5	28.5	14.2	14.2
5.4	29.9	9.4	5.3
116	7.7	23.1	11
76.4	26.7	22.3	11.8
239.8	4.1	36.9	12.3
75.3	20.3	32.5	11.3
68.4	44.5	35.6	13.6
213.5	43	33.8	21.7
193.2	18.4	65.7	15.2
76.3	27.5	16	12
110.7	40.6	63.2	16
88.3	25.5	73.4	12.9
109.8	47.8	51.4	16.7
134.3	4.9	9.3	11.2
28.6	1.5	33	7.3
217.7	33.5	59	19.4
250.9	36.5	72.3	22.2
107.4	14	10.9	11.5
163.3	31.6	52.9	16.9
197.6	3.5	5.9	11.7
184.9	21	22	15.5
289.7	42.3	51.2	25.4
135.2	41.7	45.9	17.2
222.4	4.3	49.8	11.7
296.4	36.3	100.9	23.8
280.2	10.1	21.4	14.8
187.9	17.2	17.9	14.7
238.2	34.3	5.3	20.7
137.9	46.4	59	19.2
25	11	29.7	7.2
90.4	0.3	23.2	8.7
13.1	0.4	25.6	5.3
255.4	26.9	5.5	19.8
225.8	8.2	56.5	13.4
241.7	38	23.2	21.8
175.7	15.4	2.4	14.1
209.6	20.6	10.7	15.9
78.2	46.8	34.5	14.6
75.1	35	52.7	12.6
139.2	14.3	25.6	12.2
76.4	0.8	14.8	9.4
125.7	36.9	79.2	15.9
19.4	16	22.3	6.6
141.3	26.8	46.2	15.5
18.8	21.7	50.4	7
224	2.4	15.6	11.6
123.1	34.6	12.4	15.2
229.5	32.3	74.2	19.7
87.2	11.8	25.9	10.6
7.8	38.9	50.6	6.6
80.2	0	9.2	8.8
220.3	49	3.2	24.7
59.6	12	43.1	9.7
0.7	39.6	8.7	1.6
265.2	2.9	43	12.7
8.4	27.2	2.1	5.7
219.8	33.5	45.1	19.6
36.9	38.6	65.6	10.8
48.3	47	8.5	11.6
25.6	39	9.3	9.5
273.7	28.9	59.7	20.8
43	25.9	20.5	9.6
184.9	43.9	1.7	20.7
73.4	17	12.9	10.9
193.7	35.4	75.6	19.2
220.5	33.2	37.9	20.1
104.6	5.7	34.4	10.4
96.2	14.8	38.9	11.4
140.3	1.9	9	10.3
240.1	7.3	8.7	13.2
243.2	49	44.3	25.4
38	40.3	11.9	10.9
44.7	25.8	20.6	10.1
280.7	13.9	37	16.1
121	8.4	48.7	11.6
197.6	23.3	14.2	16.6
171.3	39.7	37.7	19
187.8	21.1	9.5	15.6
4.1	11.6	5.7	3.2
93.9	43.5	50.5	15.3
149.8	1.3	24.3	10.1
11.7	36.9	45.2	7.3
131.7	18.4	34.6	12.9
172.5	18.1	30.7	14.4
85.7	35.8	49.3	13.3
188.4	18.1	25.6	14.9
163.5	36.8	7.4	18
117.2	14.7	5.4	11.9
234.5	3.4	84.8	11.9
17.9	37.6	21.6	8
206.8	5.2	19.4	12.2
215.4	23.6	57.6	17.1
284.3	10.6	6.4	15
50	11.6	18.4	8.4
164.5	20.9	47.4	14.5
19.6	20.1	17	7.6
168.4	7.1	12.8	11.7
222.4	3.4	13.1	11.5
276.9	48.9	41.8	27
248.4	30.2	20.3	20.2
170.2	7.8	35.2	11.7
276.7	2.3	23.7	11.8
165.6	10	17.6	12.6
156.6	2.6	8.3	10.5
218.5	5.4	27.4	12.2
56.2	5.7	29.7	8.7
287.6	43	71.8	26.2
253.8	21.3	30	17.6
205	45.1	19.6	22.6
139.5	2.1	26.6	10.3
191.1	28.7	18.2	17.3
286	13.9	3.7	15.9
18.7	12.1	23.4	6.7
39.5	41.1	5.8	10.8
75.5	10.8	6	9.9
17.2	4.1	31.6	5.9
166.8	42	3.6	19.6
149.7	35.6	6	17.3
38.2	3.7	13.8	7.6
94.2	4.9	8.1	9.7
177	9.3	6.4	12.8
283.6	42	66.2	25.5
232.1	8.6	8.7	13.4

Python代码

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.model_selection import KFold
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
 
# 1、创建数据集
# 1.1 通过read_excel来读取我们的目的数据集
data = pd.read_excel("D:/dataAnalysis/MachineLearning/Advertising.xlsx")
# 1.2 得到我们所需要的数据集且查看其前几列以及数据形状
print('data head:\n',data.head(),'\nShape:',data.shape)
# 1.3 数据描述
print('data describe:\n',data.describe())
# 1.4 缺失值检验
print('data isnull:\n',data[data.isnull()==True].count())
# 1.5 绘图 
data.boxplot()
plt.savefig("boxplot.jpg")
plt.show()
# 1.6 计算数据集的相关系数矩阵
#     相关系数矩阵 r(相关系数) = x和y的协方差/(x的标准差*y的标准差) == cov（x,y）/σx*σy
#     相关系数0~0.3弱相关0.3~0.6中等程度相关0.6~1强相关
print(data.corr())
 
# 1.7 建立散点图来查看数据集里的数据分布
#     seaborn的pairplot函数绘制X的每一维度和对应Y的散点图。通过设置size和aspect参数来调节显示的大小和比例。
#     可以从图中看出，TV特征和销量是有比较强的线性关系的，而Radio和Sales线性关系弱一些，Newspaper和Sales线性关系更弱。
#     通过加入一个参数kind='reg'，seaborn可以添加一条最佳拟合直线和95%的置信带。
sns.pairplot(data, x_vars=['TV','radio','newspaper'], y_vars='sales', size=7, aspect=0.8,kind = 'reg')
plt.savefig("pairplot.jpg")
plt.show()

# 2 将原数据集拆分训练集和测试集
#   利用sklearn里面的包来对数据集进行划分，以此来创建训练集和测试集
#   train_size表示训练集所占总数据集的比例
X_train,X_test,Y_train,Y_test = train_test_split(data[['TV','radio','newspaper']],data.sales,train_size=.80)
#X_train为训练数据标签,X_test为测试数据标签,exam_X为样本特征,exam_y为样本标签，train_size 训练数据占比
print(f'原始数据特征: {data.shape} ,训练数据特征: {X_train.shape} ,测试数据特征: {X_test.shape};')
print(f'原始数据标签: {data.sales.shape} ,训练数据标签: {Y_train.shape} ,测试数据标签: {Y_test.shape}.')

# 3 多元线性回归 
#   回归方程 y = a0 + a1*x1 + a2*x2 + a3*x3(模型建立最佳拟合线)
#   点误差 = 实际值 - 拟合值
#   误差平方和（Sum of square error） SSE = Σ（实际值-预测值）^2
#   最小二乘法 ： 使得误差平方和最小（最佳拟合）
#   3.1 调用线性回归API
model = LinearRegression()
#   3.2 求优化参数
model.fit(X_train,Y_train)
a  = model.intercept_#截距
b = model.coef_#回归系数
print("最佳拟合线:截距",a,",回归系数：",b)
 
#   3.3 R方检测
#   决定系数r平方
#   对于评估模型的精确度
#   y误差平方和 = Σ(y实际值 - y预测值)^2
#   y的总波动 = Σ(y实际值 - y平均值)^2
#   有多少百分比的y波动没有被回归拟合线所描述 = SSE/总波动
#   有多少百分比的y波动被回归线描述 = 1 - SSE/总波动 = 决定系数R平方
#   对于决定系数R平方来说
#   1） 回归线拟合程度：有多少百分比的y波动刻印有回归线来描述(x的波动变化)
#   2）值大小：R平方越高，回归模型越精确(取值范围0~1)，1无误差，0无法完成拟合
score = model.score(X_test,Y_test)
print(score)
 
#   3.4对线性回归进行预测
Y_pred = model.predict(X_test)
print(Y_pred)
plt.plot(range(len(Y_pred)),Y_pred,'b',label="predict")
#   显示图像
#   plt.savefig("predict.jpg")
plt.show()
 
plt.figure()
plt.plot(range(len(Y_pred)),Y_pred,'b',label="predict")
plt.plot(range(len(Y_pred)),Y_test,'r',label="test")
plt.legend(loc="upper right") #显示图中的标签
plt.xlabel("the number of sales")
plt.ylabel('value of sales')
plt.savefig("ROC.jpg")
plt.show()