引言

扩展数据结构和之前开发的例程，以支持多个特征。有几个例程被更新，使得实验看起来有些冗长，但实际上只是对之前的例程进行了小的调整，因此快速回顾是可行的

一、多变量线性回归

1.1 目标

1. 扩展我们的回归模型例程以支持多个特征
2. 扩展数据结构以支持多个特征
3. 重写预测、成本和梯度例程以支持多个特征
4. 利用NumPy的np.dot来向量化它们的实现，以提高速度和简化性

1.2 工具

在本实验中，我们将使用以下工具：

• NumPy，一个流行的科学计算库
• Matplotlib，一个流行的数据绘图库

import copy, math
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('./deeplearning.mplstyle')
np.set_printoptions(precision=2)  # 减少NumPy数组的显示精度

二、问题陈述

使用房屋价格预测的启发示例。训练数据集包含三个具有四个特征（面积、卧室数量、楼层数和房屋年龄）的示例，如下表所示。请注意，与之前的实验不同，面积是以平方英尺（sqft）为单位，而不是1000平方英尺

面积（平方英尺）	卧室数量	楼层数	房屋年龄	价格（千元美元）
2104	5	1	45	460
1416	3	2	40	232
852	2	1	35	178

使用这些值构建一个线性回归模型，以便然后预测其他房屋的价格。例如，一个面积为1200平方英尺、3个卧室、1个楼层、40年历史的房屋。

#创建X_train和y_train变量
X_train = np.array([[2104, 5, 1, 45], [1416, 3, 2, 40], [852, 2, 1, 35]])
y_train = np.array([460, 232, 178])

2.1 包含示例的矩阵X

类似于上面的表格，示例存储在NumPy矩阵X_train中。矩阵的每一行代表一个示例。当有𝑚个训练示例（在我们的例子中𝑚是三个），并且有𝑛个特征（在我们的例子中有四个），矩阵𝐗的维度为(𝑚, 𝑛)（m行，n列）。

矩阵解释：
𝐱(𝑖)是包含示例i的向量。𝐱(𝑖)=(𝑥(𝑖)0,𝑥(𝑖)1,⋯,𝑥(𝑖)𝑛−1)
𝑥(𝑖)𝑗是示例i中的元素j。方括号内的上标表示示例编号，而下标表示元素。
显示输入数据。

# 数据存储在NumPy数组/矩阵中
print(f"X形状: {X_train.shape}, X类型: {type(X_train)})")
print(X_train)
print(f"y形状: {y_train.shape}, y类型: {type(y_train)})")
print(y_train)

输出结果：

2.2 参数向量w, b

𝐰是一个具有𝑛个元素的向量。

• 每个元素都包含与一个特征相关的参数。
• 在我们的数据集中，n是4。
• 理论上，我们将其绘制为列向量
  

𝑏是一个标量参数。
为了演示，𝐰 和 𝑏 将加载一些初始选定的值，这些值接近最优值。𝐰 是一个1-D NumPy向量

b_init = 785.1811367994083
w_init = np.array([ 0.39133535, 18.75376741, -53.36032453, -26.42131618])
print(f"w_init形状: {w_init.shape},

输出结果：

三、多变量模型预测

具有多个变量的模型的预测由线性模型给出：
𝑓𝐰,𝑏(𝐱)=𝑤0𝑥0+𝑤1𝑥1+…+𝑤𝑛−1𝑥𝑛−1+𝑏(1)或者用向量表示法：𝑓𝐰,𝑏(𝐱)=𝐰⋅𝐱+𝑏(2)
其中 ⋅ 是向量点积
为了演示点积，我们将使用(1)和(2)来实现预测。

3.1 逐个元素进行单个预测

我们之前的预测是将一个特征值乘以一个参数，然后加上偏置参数。将我们之前预测的实现扩展到多个特征，可以通过遍历每个元素，执行其参数的乘法，然后在最后加上偏置参数来实现(1)上述。

def predict_single_loop(x, w, b): 
    """
    使用线性回归进行单个预测
    参数:
      x (ndarray): 形状为 (n,) 的多个特征示例
      w (ndarray): 形状为 (n,) 的模型参数    
      b (标量):  模型参数     
      
    返回:
      p (标量):  预测
    """
    n = x.shape[0]
    p = 0
    for i in range(n):
        p_i = x[i] * w[i]  
        p = p + p_i         
    p = p + b                
    return p

# 从我们的训练数据中获取一行
x_vec = X_train[0,:]
print(f"x_vec 形状 {x_vec.shape}, x_vec 值: {x_vec}")

# 进行预测
f_wb = predict_single_loop(x_vec, w_init, b_init)
print(f"f_wb 形状 {f_wb.shape}, 预测: {f_wb}")

输出结果：

注意$x_{v}ec$的形状。它是一个具有4个元素的1-D NumPy向量，形状为(4,)。结果$f_{w}b$是一个标量。

3.2 单个预测向量

注意方程(1)可以使用向量点积(2)来实现。我们可以利用向量运算来加速预测。
回顾Python/NumPy实验室，NumPy的np.dot()可以用来执行向量点积。

def predict(x, w, b): 
    """
    使用线性回归进行单个预测
    参数:
      x (ndarray): 形状为 (n,) 的多个特征示例
      w (ndarray): 形状为 (n,) 的模型参数   
      b (标量): 模型参数 
      
    返回:
      p (标量): 预测
    """
    p = np.dot(x, w) + b     
    return p    
# 从我们的训练数据中获取一行
x_vec = X_train[0,:]
print(f"x_vec 形状 {x_vec.shape}, x_vec 值: {x_vec}")
# 进行预测
f_wb = predict(x_vec,w_init, b_init)
print(f"f_wb 形状 {f_wb.shape}, 预测: {f_wb}")

输出结果：

结果和形状与之前使用循环的版本相同。从现在开始，np.dot将被用于这些操作。预测现在是一个单一的语句。大多数例程将直接实现它，而不是调用一个单独的预测例程

四、使用多变量计算成本

多变量成本函数的方程是：
$$J(\mathbf{w},b)=\frac{1}{2m}\sum _{i=0}^{m-1}(f_{\mathbf{w},b}({\mathbf{x}}^{(i)})-y^{(i)}{)}^2\text{\hspace{1em}(3)}$$

其中：
$$f_{\mathbf{w},b}({\mathbf{x}}^{(i)})=\mathbf{w}\cdot {\mathbf{x}}^{(i)}+b\text{\hspace{1em}(4)}$$
与之前的实验室相比，𝐰和${\mathbf{x}}^{(i)}$是向量而不是标量，支持多个特征。
下面是方程(3)和(4)的实现。请注意，这使用了本课程中的一种标准模式，其中使用一个for循环遍历所有m个示例。

def compute_cost(X, y, w, b): 
    """
    计算成本
    参数:
      X (ndarray (m,n)): 数据，m个示例具有n个特征
      y (ndarray (m,)): 目标值
      w (ndarray (n,)): 模型参数  
      b (标量): 模型参数
      
    返回:
      成本 (标量): 成本
    """
    m = X.shape[0]
    成本 = 0.0
    for i in range(m):                                
        f_wb_i = np.dot(X[i], w) + b           #(n,)(n,) = 标量 (见 np.dot)
        成本 = 成本 + (f_wb_i - y[i])**2       #标量
    成本 = 成本 / (2 * m)                      #标量    
    return 成本

# 使用我们预先选择的优化参数计算并显示成本。
成本 = compute_cost(X_train, y_train, w_init, b_init)
print(f'在最佳w处的成本 : {成本}')

输出结果：

预期结果：最佳w处的成本 : 1.5578904045996674e-12

五、带有多个变量的梯度下降

对于多个变量的梯度下降：
重复，直到收敛：
repeat  until convergence:    {    w j = w j − α ∂ J ( w , b ) ∂ w j    for j = 0..n-1 b    = b − α ∂ J ( w , b ) ∂ b } \begin{align*} \text{repeat}&\text{ until convergence:} \; \lbrace \newline\; & w_j = w_j - \alpha \frac{\partial J(\mathbf{w},b)}{\partial w_j} \tag{5} \; & \text{for j = 0..n-1}\newline &b\ \ = b - \alpha \frac{\partial J(\mathbf{w},b)}{\partial b} \newline \rbrace \end{align*} repeat}​ until convergence:{wj​=wj​−α∂wj​∂J(w,b)​b  =b−α∂b∂J(w,b)​​for j = 0..n-1(5)​

其中，n 是特征的数量，参数 𝑤𝑗 和 𝑏 同时更新，其中
$$\begin{array}{rl}\frac{\partial J(\mathbf{w},b)}{\partial w_j}& =\frac{1}{m}\sum _{i=0}^{m-1}(f_{\mathbf{w},b}({\mathbf{x}}^{(i)})-y^{(i)})x_j^{(i)}\\ \frac{\partial J(\mathbf{w},b)}{\partial b}& =\frac{1}{m}\sum _{i=0}^{m-1}(f_{\mathbf{w},b}({\mathbf{x}}^{(i)})-y^{(i)})\end{array}\text{\hspace{1em}(6)(7)}$$

• 𝑚 是数据集中的训练示例数量
• $f_{\mathbf{w},b}({\mathbf{x}}^{(i)})$是模型的预测，而 $y^{(i)}$是目标值

5.1 使用多个变量计算梯度

下面是计算方程（6）和（7）的实现。有多种实现方式。在这个版本中，有一个
外循环遍历所有m个示例。
$\frac{\partial J(\mathbf{w},b)}{\partial b}$对于示例可以直接计算并累加在第二个循环中，遍历所有n个特征：
$\frac{\partial J(\mathbf{w},b)}{\partial w_j}$为每个$w_j$计算。

def compute_gradient(X, y, w, b): 
    """
    计算线性回归的梯度
    参数:
      X (ndarray (m,n)): 数据，m个示例具有n个特征
      y (ndarray (m,)): 目标值
      w (ndarray (n,)): 模型参数  
      b (标量): 模型参数
      
    返回:
      dj_dw (ndarray (n,)): 成本关于参数w的梯度。 
      dj_db (标量): 成本关于参数b的梯度。 
    """
    m,n = X.shape           #(示例数量, 特征数量)
    dj_dw = np.zeros((n,))
    dj_db = 0.
    for i in range(m):                             
        err = (np.dot(X[i], w) + b) - y[i]   
        for j in range(n):                         
            dj_dw[j] = dj_dw[j] + err * X[i, j]    
        dj_db = dj_db + err                        
    dj_dw = dj_dw / m                                
    dj_db = dj_db / m                                
        
    return dj_db, dj_dw

# 计算和显示梯度
tmp_dj_db, tmp_dj_dw = compute_gradient(X_train, y_train, w_init, b_init)
print(f'初始w,b处的dj_db: {tmp_dj_db}')
print(f'初始w,b处的dj_dw: \n {tmp_dj_dw}')

输出结果：

预期结果：
初始w,b处的dj_db: -1.6739251122999121e-06
初始w,b处的dj_dw: [-2.73e-03 -6.27e-06 -2.22e-06 -6.92e-05]

5.2 带有多个变量的梯度下降

下面的例行程序实现了上面的方程（5）。

def gradient_descent(X, y, w_in, b_in, cost_function, gradient_function, alpha, num_iters): 
    """
    执行批量梯度下降以学习theta。通过采取num_iters个梯度步骤并使用学习率alpha来更新theta。
    
    参数:
      X (ndarray (m,n))   : 数据，m个示例具有n个特征
      y (ndarray (m,))    : 目标值
      w_in (ndarray (n,)) : 初始模型参数  
      b_in (标量)       : 初始模型参数
      cost_function       : 计算成本的函数
      gradient_function   : 计算梯度的函数
      alpha (float)       : 学习率
      num_iters (int)     : 运行梯度下降的迭代次数
      
    返回:
      w (ndarray (n,)) : 参数的更新值
      b (标量)       : 参数的更新值
      """
    
    # 一个数组来存储每次迭代的成本J和w，主要用于后续的绘图
    J_history = []
    w = copy.deepcopy(w_in)  # 在函数内部避免修改全局w
    b = b_in
    
    for i in range(num_iters):
        # 计算梯度并更新参数
        dj_db,dj_dw = gradient_function(X, y, w, b)   ##None
        
        # 使用w, b, alpha和梯度更新参数
        w = w - alpha * dj_dw               ##None
        b = b - alpha * dj_db               ##None
        
        # 每迭代一次就保存成本J
        if i<100000:      # 防止资源耗尽
            J_history.append(cost_function(X, y, w, b))
        
        # 每迭代10次打印一次成本，或者如果迭代次数小于10，则打印所有迭代
        if i% math.ceil(num_iters / 10) == 0:
            print(f"迭代 {i:4d}: 成本 {J_history[-1]:8.2f}   ")
        
    return w, b, J_history #返回最终的w,b和J历史，用于绘图

在下一个单元格中，您将测试该实现：

# 初始化参数
initial_w = np.zeros_like(w_init)
initial_b = 0.
# 一些梯度下降设置
iterations = 1000
alpha = 5.0e-7
# 运行梯度下降
w_final, b_final, J_hist = gradient_descent(X_train, y_train, initial_w, initial_b,
                                                    compute_cost, compute_gradient, 
                                                    alpha, iterations)
print(f"通过梯度下降找到的b,w: {b_final:0.2f},{w_final} ")
m, _ = X_train.shape
for i in range(m):
    print(f"预测: {np.dot(X_train[i], w_final) + b_final:0.2f}, 目标值: {y_train[i]}")

输出结果：

预期结果：
通过梯度下降找到的b,w: -0.00, [ 0.2 0. -0.01 -0.07]
预测: 426.19, 目标值: 460
预测: 286.17, 目标值: 232
预测: 171.47, 目标值: 178

# 绘制成本与迭代次数的关系图
fig, (ax1, ax2) = plt.subplots(1, 2, constrained_layout=True, figsize=(12, 4))
ax1.plot(J_hist)
ax2.plot(100 + np.arange(len(J_hist[100:])), J_hist[100:])
ax1.set_title("Cost vs. iteration");  ax2.set_title("Cost vs. iteration (tail)")
ax1.set_ylabel('Cost')             ;  ax2.set_ylabel('Cost') 
ax1.set_xlabel('iteration step')   ;  ax2.set_xlabel('iteration step') 
plt.show()

输出结果：

成本仍在下降，我们的预测并不十分准确，后续继续优化

总结

1. 重新开发了线性回归的例程，现在支持多个变量
2. 利用NumPy的np.dot来向量化实现