多元线性回归的公式与矢量化代码实现方法（笔记）¶

以房价预测为例

面积（sqft）	卧室数	楼层数	房龄	价格（千美元）
2104	5	1	45	460
1416	3	2	40	232
852	2	1	35	178

一共有三组样本，m=3；每组样本的x有4个特征，分别是面积，卧室数，楼层数，房龄。价格为目标值y。
\[ f_{\mathbf{w},b}(\mathbf{x}) = \mathbf{w} \cdot \mathbf{x} + b \tag{1} \]

1.创建样本矩阵与参数变量初始化¶

\[\mathbf{X} = \begin{pmatrix} x^{(0)}_0 & x^{(0)}_1 & \cdots & x^{(0)}_{n-1} \\ x^{(1)}_0 & x^{(1)}_1 & \cdots & x^{(1)}_{n-1} \\ \cdots \\ x^{(m-1)}_0 & x^{(m-1)}_1 & \cdots & x^{(m-1)}_{n-1} \end{pmatrix}\]

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

\[ \mathbf{w} = \begin{pmatrix} w_0 \\ w_1 \\ \cdots\\ w_{n-1} \end{pmatrix} \]

b为标量

#此初始化值接近最优
b_init = 785.1811367994083
w_init = np.array([ 0.39133535, 18.75376741, -53.36032453, -26.42131618])

2.计算成本函数¶

\[ J(\mathbf{w},b) = \frac{1}{2m} \sum\limits_{i = 0}^{m-1} \big(f_{\mathbf{w},b}(\mathbf{x}^{(i)}) - y^{(i)}\big)^2 \tag{2} \]

\[ f_{\mathbf{w},b}(\mathbf{x}^{(i)}) = \mathbf{w} \cdot \mathbf{x}^{(i)} + b \tag{3} \]

（2）为多变量下的成本函数，（3）为多变量下的回归方程。

def compute_cost(X, y, w, b):  #定义一个函数用于计算成本
    m = X.shape[0]             #将样本数赋给m
    cost = 0.0
    for i in range(m):         #用一个for循环计算每组样本预测值与目标值差的平方并求和
        f_wb_i = np.dot(X[i], w) + b         #计算样本预测值  
        cost = cost + (f_wb_i - y[i])**2       
    cost = cost / (2 * m)                      
    return cost

3.计算多变量的梯度¶

\[ \begin{align} \frac{\partial J(\mathbf{w},b)}{\partial w_j} &= \frac{1}{m} \sum\limits_{i=0}^{m-1} \big(f_{\mathbf{w},b}(\mathbf{x}^{(i)}) - y^{(i)}\big) x_{j}^{(i)} \tag{4} \\ \frac{\partial J(\mathbf{w},b)}{\partial b} &= \frac{1}{m} \sum\limits_{i=0}^{m-1} \big(f_{\mathbf{w},b}(\mathbf{x}^{(i)}) - y^{(i)}\big) \tag{5} \end{align} \]

def compute_gradient(X, y, w, b):   #定义一个函数用于计算梯度
    m,n = X.shape                           #样本数赋给m，特征数赋给n
    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]    #将目标与预测之差乘特征x得到特征w梯度的中间量
        dj_db = dj_db + err           #得到J对于b梯度的中间量             
    dj_dw = dj_dw / m                 #得到J对于w梯度的向量               
    dj_db = dj_db / m                 #得到J对于b的梯度              

    return dj_db, dj_dw

4.计算多变量的梯度下降¶

\[ \begin{align*} \text{repeat}&\text{ until convergence:}\;\{ \\[-2mm] & w_j = w_j - \alpha \frac{\partial J(\mathbf{w},b)}{\partial w_j} \tag{6} \quad \text{for } j=0..n-1 \\[1mm] & b\ \ = b - \alpha \frac{\partial J(\mathbf{w},b)}{\partial b} \\ \} \end{align*} \]

def gradient_descent(X, y, w_in, b_in, cost_function, gradient_function, alpha, num_iters):       #定义一个函数来迭代W与b，num_iters表示迭代次数
    # 用于在每次迭代中存储代价 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)   # 就是上一步中计算梯度的公式

        # 使用 w、b、学习率 alpha 和梯度更新参数
        w = w - alpha * dj_dw               # 向量
        b = b - alpha * dj_db               # 标量

        # 在每次迭代时保存成本 J
        if i<100000:      # 防止资源耗尽
            J_history.append( cost_function(X, y, w, b))

        # 每隔一定步数打印一次代价，步数为迭代次数除以10（向上取整），最多打印 10 次；若迭代次数少于 10，则按实际次数打印
        if i% math.ceil(num_iters / 10) == 0:
            print(f"Iteration {i:4d}: Cost {J_history[-1]:8.2f}   ")

    return w, b, J_history # 返回最终的 w、b 以及用于绘图的 J 历史

编辑于9月22日