机器学习基础学习-多元线性回归问题（梯度下降法实现）_多元线性回归梯度下降法

1、基本概念

在之前的博客当中描述了怎样模拟出了梯度下降的过程

如果是多维情况，theta其实是一个向量，那么对其求导的损失函数也是向量，梯度就是损失函数对每个方向的theta求偏导。
和之前的一维线性回归相比，我们对只是对w这个数字进行求导，而现在针对多系数theta，对theta整个向量进行求导。

有两个参数的梯度下降法进行可视化，一圈一圈代表等高线，圈上的值就是梯度。越外层J的取值越大，越里层J的取值越小，在图的中心位置达到损失函数的最小值


对每一个维度的theta进行求导（这里的Xb注意依然是加上了一列恒为1的一列，代表截距）

在上面的式子当中，可以看到求导之后的结果其实是一个求和的过程，也就是说在这里如果样本数量越大，梯度也就越大，这显然是不合理的，最后我们的梯度应该和样本数量m是无关的。所以下面改变目标函数（除以样本数量m，再提取公共的2），最后发现这个结果就是y的预测值和真值之间的MSE值（均方误差）

2、代码实现

这里为了可视化方便，采取的数据都是一维的数据，但是实际上是对矩阵进行处理。

（1）损失函数

np.sum((y - X_b.dot(theta)) ** 2 / len(X_b)
1

（2）求梯度

# 梯度(比较笨的方法)
def dJ(theta, X_b, y):
  res = np.empty(len(theta)) # 开一个和theta一样大的空间(因为要对theta的每一个元素求偏导)
  res[0] = np.sum(X_b.dot(theta) - y) # 第一行单独写，其他的要写一个循环
  for i in range(1, len(theta)):
    res[i] = (X_b.dot(theta) - y).dot(X_b[:,i]) # X_b[:,i]相当于取出第i列，也就是第i个特征值
  return res * 2 / len(X_b)
1
2
3
4
5
6
7

（3）X_b

# X添加一列样本个数行、1列，每个取值都是1的一列
X_b = np.hstack([np.ones((len(X), 1)), X])
1
2

（4）总体代码

# 多维线性回归问题（梯度下降法）

import numpy as np
import matplotlib.pyplot as plt

# 设置随机种子
np.random.seed(666)
# 用一维数据方便可视化，但实际是对矩阵进行操作
x = np.random.random(size = 100)
y = x * 3. + 4. + np.random.normal(size = 100) # 噪音：均值为0，方差为1

X = x.reshape(-1, 1) # 数组array中的方法，作用是将数据重新组织成100行1列的数据

# print('X数组', X.shape)
# print('y数组', y.shape)

plt.scatter(x, y)
plt.show()

# 损失函数
def J(theta, X_b, y):
  try:
    return np.sum((y - X_b.dot(theta)) ** 2 / len(X_b))
  except:
    return float('inf') # 返回float最大值

# 梯度(比较笨的方法)
def dJ(theta, X_b, y):
  res = np.empty(len(theta)) # 开一个和theta一样大的空间(因为要对theta的每一个元素求偏导)
  res[0] = np.sum(X_b.dot(theta) - y) # 第一行单独写，其他的要写一个循环
  for i in range(1, len(theta)):
    res[i] = (X_b.dot(theta) - y).dot(X_b[:,i]) # X_b[:,i]相当于取出第i列，也就是第i个特征值
  return res * 2 / len(X_b)

# 直接复制之前模拟梯度下降过程的代码,然后进行更改
def gradient_descent(X_b, y, initial_theta, eta, n_iters = 1e4, epsilon = 1e-8):
  theta = initial_theta
  i_iters = 0

  while i_iters < n_iters:
      gradient = dJ(theta, X_b, y) # 求梯度
      last_theta = theta # theta重新赋值前，记录上一场的值
      theta = theta - eta * gradient # 通过一定的eta学习率取得下一个点的theta
      # 最近两点的损失函数差值小于一定精度，退出循环
      if(abs(J(theta, X_b, y) - J(last_theta, X_b, y)) < epsilon):
          break
      i_iters += 1
  return theta


''' 
  初始化数值
'''
# X添加一列样本个数行、1列，每个取值都是1的一列
X_b = np.hstack([np.ones((len(X), 1)), X])
initial_theta = np.zeros(X_b.shape[1]) # 设置n+1维的向量,X_b.shape[1]:第一行的维数
eta = 0.1

theta = gradient_descent(X_b, y, initial_theta, eta)

print('theta', theta)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

运行结果

这里的theta[0]就是截距，theta[1]就是斜率，和我们一开始设置的一致，说明梯度下降法成功训练了我们的模型

y = x * 3. + 4. + np.random.normal(size = 100)
1

最后

所有的思想都在上面，最后我把多元线性回归的两种方法总结在了一起

# 多元线性回归问题；解析解+梯度下降法（也能解决一元线性回归问题）
import numpy as np
from sklearn import datasets
from sklearn.metrics import r2_score
from sklearn.model_selection import train_test_split

class LinearRegression:
  def __init__(self):
      # 初始化Linear Regression模型
      self.coef_ = None # 系数，对应theta1-n，对应的向量
      self.interception_ = None # 截距，对应theta0
      self._theta = None # 定义私有变量，整体计算的theta



  '''
    梯度下降
  '''
  def fit_gd(self, X_train, y_train, eta = 0.01, n_iters = 1e4): 
    # 根据训练数据集X_train, y_ .train训练Linear Regression模型
    # X_train的样本数量和y_train的标记数量应该是一致的
    # 使用shape[0]读取矩阵第一维度的长度，在这里就是列数
    assert X_train.shape[0] == y_train.shape[0], \
    "the size of x_ .train must be equal to the size of y_ train"
    # 损失函数
    def J(theta, X_b, y):
      try:
        return np.sum((y - X_b.dot(theta)) ** 2 / len(y))
      except:
        return float('inf') # 返回float最大值

    # 梯度(比较笨的方法)
    def dJ(theta, X_b, y):
      res = np.empty(len(theta)) # 开一个和theta一样大的空间(因为要对theta的每一个元素求偏导)
      res[0] = np.sum(X_b.dot(theta) - y) # 第一行单独写，其他的要写一个循环
      for i in range(1, len(theta)):
        res[i] = (X_b.dot(theta) - y).dot(X_b[:,i]) # X_b[:,i]相当于取出第i列，也就是第i个特征值
      return res * 2 / len(X_b)


    # 求解theta矩阵
    def gradient_descent(X_b, y, initial_theta, eta, n_iters = 1e4, epsilon = 1e-8):
      theta = initial_theta
      cur_iters = 0

      while cur_iters < n_iters:
          gradient = dJ(theta, X_b, y) # 求梯度
          last_theta = theta # theta重新赋值前，记录上一场的值
          theta = theta - eta * gradient # 通过一定的eta学习率取得下一个点的theta
          # 最近两点的损失函数差值小于一定精度，退出循环
          if(abs(J(theta, X_b, y) - J(last_theta, X_b, y)) < epsilon):
              break
          cur_iters += 1
      return theta

    # 得到X_b
    X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
    initial_theta = np.zeros(X_b.shape[1]) # 设置n+1维的向量,X_b.shape[1]:第一行的维数
    # X_b.T是X_b的转置,.dot是点乘,np.linalg.inv是求逆
    # 获取theta
    self._theta = gradient_descent(X_b, y_train, initial_theta, eta, n_iters)

    self.interception_ = self._theta[0] # 截距
    self.coef_ = self._theta[1:] # 系数 
    return self




  '''
    解析解
  '''
  # 正规化方程
  # 训练过程（解析解）
  def fit_normal(self, X_train, y_train): # X_train和y_train都是矩阵
    # 根据训练数据集X_train, y_ .train训练Linear Regression模型
    # X_train的样本数量和y_train的标记数量应该是一致的
    # 使用shape[0]读取矩阵第一维度的长度，在这里就是列数
    assert X_train.shape[0] == y_train.shape[0], \
    "the size of x_ .train must be equal to the size of y_ train"
    # np.hstack():在水平方向上平铺，就是在横向上多加一列
    # np.ones(矩阵大小, 列数)是增加一列恒为1的一列
    # 得到X_b
    X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
    # X_b.T是X_b的转置,.dot是点乘,np.linalg.inv是求逆
    self._theta = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y_train)

    self.interception_ = self._theta[0] # 截距
    self.coef_ = self._theta[1:] # 系数
    return self
  
  '''
    预测过程
  '''
  def predict(self, X_predict):
    # 给定待预测数据集X_predict,返回表示X_predict的结果向量
    X_b = np.hstack([np.ones((len(X_predict), 1)), X_predict])
    return X_b.dot(self._theta)

  '''
    确定精度，评价多元线性回归的结果
  '''
  def score(self, X_test, y_test):
    # 根据测试数据集X_test 和y_test 确定当前模型的准确度
    y_predict = self.predict(X_test)
    return r2_score(y_test, y_predict) # r2_score求真值y_test和预测值y_predict的r方

  '''
    显示属性
  '''
  def __repr__(self):
      return "LinearRegression()"


# 加载波士顿房价数据集，并划分为X_train,y_train
# 波士顿房价
boston = datasets.load_boston()

# X_train = boston.data
# y_train = boston.target

X = boston.data
y = boston.target

X = X[y < 50.0]
y = y[y < 50.0]
# 分出测试集合训练集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25)

'''
  实例化(数学解解法)
'''
reg = LinearRegression()
# 进行训练
reg.fit_normal(X_train, y_train) # LinearRegression()
# 预测r方值
print('预测精确度', reg.score(X_test, y_test))
# 截距
print('数学解系数reg.coef_', reg.coef_)
print('数学解截距reg.interception_', reg.interception_)
print('数学解预测值', reg.predict(X_test))
print('y_test测试值', y_test)



'''
  实例化(梯度下降解法)
'''
lin_reg = LinearRegression()
# 设置随机种子
np.random.seed(666)
# 用一维数据方便可视化，但实际是对矩阵进行操作
x = np.random.random(size = 100)
y = x * 3. + 4. + np.random.normal(size = 100) # 噪音：均值为0，方差为1

X = x.reshape(-1, 1) # 数组array中的方法，作用是将数据重新组织成100行1列的数据
# 梯度下降法进行训练
lin_reg.fit_gd(X, y) # LinearRegression()

print('梯度下降系数reg.coef_', lin_reg.coef_)
print('梯度下降截距reg.interception_', lin_reg.interception_)




1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165

运行结果