机器学习算法实现：从零开始的线性回归

来源：网络更新时间：2024-06-07 15:31:40

最近，我决定将本科学过的一些机器学习算法从零开始实现一下。虽然天天当调包侠也能愉快玩耍，但对于我这个强迫症患者来说，底层逻辑是必须要了解的。接下来，我会从最简单的多元线性回归开始，一步一步在jupyter里面实现。

【注1】：本文默认读者具有基本的机器学习基础、numpy基础，如果没有，建议先学习吴恩达老师的机器学习课程再来阅读本文。

【注2】：本文大部分代码会采用向量化加速编程，当然部分情况下也会用到循环遍历的情况。后面有时间会再出一篇循环遍历实现的文章。

【注3】：作者实力有限，有错在所难免，欢迎大家指出。

一、理论介绍

线性回归的预测函数可以写成向量的形式：\(y = \mathbf{w^Tx}+b\)，其中\(\mathbf{w} = (w_1,...,w_n)^T,\mathbf{x}=(x_1,...,x_n)^T\)。为了编程方便，将\(b\)收缩到\(w\)向量里面去，有\(\mathbf{w} = (w_0,w_1,...,w_n)^T,w_0=b\)。上面模型变为\(y = \mathbf{w^Tx}\)。损失函数和偏导公式如下。

有了以上的理论操作，下面可以开始操作了。上述理论推导本文不作说明，有需要可参考周志华老师的《西瓜书》。

二、代码实现

1、导入库

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

库就不用介绍了，很常见的机器学习会用到的库。

2、准备数据集

生成数据集的方式很多，random库，或者你有数据（csv等）都可以。本文重点不是这个，因此简单用numpy手打生成了。

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

注：我一般喜欢将向量是列向量还是行向量与数学公式中严格对应起来，否则，广播机制出现了问题真的很头疼的。

# 训练集第一列插入全为1的1列
new_column = np.ones((X_train.shape[0], 1))  # 创建一列全是 1 的数组
X_train = np.concatenate((new_column, X_train), axis=1)  # 在索引为 0 的位置插入新列
print(X_train)

3、定义预测函数（predict）

def predict(x,w):
    '''
    
    :param x: 要预测的样本点向量,x是行向量
    :param w: 训练出来的参数向量，w,是列向量
    :param b: 训练出来的参数 b ，b可以是标量，也可以是一个元素的数组
    :return: prediction,预测值
    '''
    return np.dot(x,w)

w_init = np.array([ [785.1811367994083],[0.39133535], [18.75376741], [-53.36032453], [-26.42131618]])
x_vec = X_train[0,:]
print(f"x_vec shape {x_vec.shape}, x_vec value: {x_vec}")

# make a prediction
f_wb = predict(x_vec,w_init)
print(f"f_wb shape {f_wb.shape}, prediction: {f_wb}")

4、代价（损失）函数

向量化加速：

def loss(X,y,w):
    m = X.shape[0]
    err = np.dot(X,w) -y
    return (1/(2*m))*np.dot(err.T,err)

cost = loss(X_train, y_train, w_init)
print(f'Cost at optimal w : {cost}')

5、计算参数梯度

向量化加速：

def compute_gradient(X, y, w):
    err = np.dot(X,w) - y
    return (1/len(X))*np.dot(X.T,err)

g = compute_gradient(X_train,y_train,w_init)
print(g)

6、批量梯度下降

def gradient_descent(X, y, w_in, loss, gradient_function, alpha, num_iters): 
    J_history = []
    # 用来存每进行一次梯度后的损失，方便查看损失变化。如果要画损失变化也是用这个
    # w = copy.deepcopy(w_in)  #avoid modifying global w within function
    w = w_in
    
    for i in range(num_iters):
        dj_dw = gradient_function(X, y, w)   
        w = w - alpha * dj_dw                              
        if i<100000:      
            J_history.append( loss(X, y, w))
        if i % math.ceil(num_iters / 10) == 0:
            # 控制间隔 打印出一次损失结果
            # 总迭代次数分成均分10组，在每组最后打印一次损失
            print(f"迭代次数（梯度下降次数） {i}: Cost {J_history[-1]}")
        
    return w,J_history

7、训练

initial_w = np.zeros_like(w_init)
iterations = 1000
alpha = 5.0e-7
# run gradient descent 
w_final,J_hist = gradient_descent(X_train, y_train, initial_w,loss, compute_gradient, 
                                                    alpha, iterations)
print(f"b,w found by gradient descent: {w_final} ")
m = X_train.shape[0]
for i in range(m):
    print(f"prediction: {np.dot(X_train[i], w_final)}, target value: {y_train[i]}")

8、可视化一下损失

J_hist = [i[0,0] for i in J_hist]
# 画一下损失变化图
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()

明显看到，随便取的数据集，loss迭代到后面就稳定600多了，线性回归不太适合该数据集。有时间在用多项式回归试试效果。