使用Python实现线性回归：从基础到scikit-learn

线性回归是机器学习中最基础也是最重要的算法之一。本文将带领读者从基础的numpy实现，到使用成熟的scikit-learn库，全面了解线性回归的实现过程。我们将通过实际的代码示例和可视化来深入理解这个算法。

1. 环境准备

首先，让我们导入所需的库并设置环境：

from __future__ import division, print_function, unicode_literals
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import warnings
np.random.seed(42)
%matplotlib inline
mpl.rc('axes', labelsize=14)
mpl.rc('xtick', labelsize=12)
mpl.rc('ytick', labelsize=12)
warnings.filterwarnings(action="ignore", message="^internal gelsd")

这段代码导入了必要的库，设置了随机种子以确保结果可重现，并配置了matplotlib的一些参数。

2. 数据准备和可视化

假设我们已经有了训练数据X和y。让我们先来可视化这些数据：

plt.plot(X, y, "b.")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.axis([0, 2, 0, 15])
plt.show()

这将绘制一个散点图，展示我们的数据分布。

3. 使用numpy实现线性回归

现在，让我们使用numpy来手动实现线性回归：

X_b = np.c_[np.ones((100, 1)), X]  # 添加x0 = 1到每个实例
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)

这里，我们首先添加了一列1到X矩阵，然后使用正规方程计算最优的theta值。

4. 使用模型进行预测

有了theta_best，我们就可以进行预测了：

X_new = np.array([[0], [2]])
X_new_b = np.c_[np.ones((2, 1)), X_new]
y_predict = X_new_b.dot(theta_best)

5. 可视化预测结果

让我们把原始数据和预测结果可视化：

plt.plot(X_new, y_predict, "r-", linewidth=2, label="Predictions")
plt.plot(X, y, "b.")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.legend(loc="upper left", fontsize=14)
plt.axis([0, 2, 0, 15])
plt.show()

这将绘制一个图，显示原始数据点和我们的预测线。

6. 使用scikit-learn实现线性回归

最后，让我们看看如何使用scikit-learn来实现相同的功能：

from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(X, y)
print("截距:", lin_reg.intercept_)
print("系数:", lin_reg.coef_)
# 预测
print("预测结果:", lin_reg.predict(X_new))

使用scikit-learn，我们只需要几行代码就可以完成模型的训练和预测。

7. 梯度下降法

除了使用正规方程，我们还可以使用梯度下降法来训练线性回归模型。以下是批量梯度下降的实现：

eta = 0.1  # 学习率
n_iterations = 1000
m = 100

theta = np.random.randn(2,1)  # 随机初始化
for iteration in range(n_iterations):
    gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)
    theta = theta - eta * gradients

print("梯度下降法得到的theta:", theta)

我们还可以可视化梯度下降的过程：

theta_path_bgd = []

def plot_gradient_descent(theta, eta, theta_path=None):
    m = len(X_b)
    plt.plot(X, y, "b.")
    n_iterations = 1000
    for iteration in range(n_iterations):
        if iteration < 10:
            y_predict = X_new_b.dot(theta)
            style = "b-" if iteration > 0 else "r--"
            plt.plot(X_new, y_predict, style)
        gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)
        theta = theta - eta * gradients
        if theta_path is not None:
            theta_path.append(theta)
    plt.xlabel("$x_1$", fontsize=18)
    plt.axis([0, 2, 0, 15])
    plt.title(r"$\eta = {}$".format(eta), fontsize=16)

np.random.seed(42)
theta = np.random.randn(2,1)  # 随机初始化

plt.figure(figsize=(10,4))
plt.subplot(131); plot_gradient_descent(theta, eta=0.02)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.subplot(132); plot_gradient_descent(theta, eta=0.1, theta_path=theta_path_bgd)
plt.subplot(133); plot_gradient_descent(theta, eta=0.5)
plt.show()

这段代码展示了不同学习率对梯度下降过程的影响。

8. 随机梯度下降和小批量梯度下降

除了批量梯度下降，我们还可以实现随机梯度下降（SGD）和小批量梯度下降：

# 随机梯度下降
theta_path_sgd = []
m = len(X_b)
np.random.seed(42)
n_epochs = 50
t0, t1 = 5, 50  # 学习率调度超参数

def learning_schedule(t):
    return t0 / (t + t1)

theta = np.random.randn(2,1)  # 随机初始化
for epoch in range(n_epochs):
    for i in range(m):
        random_index = np.random.randint(m)
        xi = X_b[random_index:random_index+1]
        yi = y[random_index:random_index+1]
        gradients = 2 * xi.T.dot(xi.dot(theta) - yi)
        eta = learning_schedule(epoch * m + i)
        theta = theta - eta * gradients
        theta_path_sgd.append(theta)

# 小批量梯度下降
theta_path_mgd = []
n_iterations = 50
minibatch_size = 20
np.random.seed(42)
theta = np.random.randn(2,1)  # 随机初始化
t0, t1 = 200, 1000

def learning_schedule(t):
    return t0 / (t + t1)

t = 0
for epoch in range(n_iterations):
    shuffled_indices = np.random.permutation(m)
    X_b_shuffled = X_b[shuffled_indices]
    y_shuffled = y[shuffled_indices]
    for i in range(0, m, minibatch_size):
        t += 1
        xi = X_b_shuffled[i:i+minibatch_size]
        yi = y_shuffled[i:i+minibatch_size]
        gradients = 2/minibatch_size * xi.T.dot(xi.dot(theta) - yi)
        eta = learning_schedule(t)
        theta = theta - eta * gradients
        theta_path_mgd.append(theta)

9. 比较不同的梯度下降方法

最后，我们可以比较不同梯度下降方法的参数路径：

theta_path_bgd = np.array(theta_path_bgd)
theta_path_sgd = np.array(theta_path_sgd)
theta_path_mgd = np.array(theta_path_mgd)

plt.figure(figsize=(7,4))
plt.plot(theta_path_sgd[:, 0], theta_path_sgd[:, 1], "r-s", linewidth=1, label="Stochastic")
plt.plot(theta_path_mgd[:, 0], theta_path_mgd[:, 1], "g-+", linewidth=2, label="Mini-batch")
plt.plot(theta_path_bgd[:, 0], theta_path_bgd[:, 1], "b-o", linewidth=3, label="Batch")
plt.legend(loc="upper left", fontsize=16)
plt.xlabel(r"$\theta_0$", fontsize=20)
plt.ylabel(r"$\theta_1$   ", fontsize=20, rotation=0)
plt.axis([2.5, 4.5, 2.3, 3.9])
plt.show()

总结

在这篇博客中，我们学习了如何使用numpy手动实现线性回归，以及如何利用scikit-learn快速实现相同的功能。我们还深入探讨了不同的梯度下降方法，包括批量梯度下降、随机梯度下降和小批量梯度下降，并通过可视化比较了它们的性能。

通过这些实现和比较，我们不仅可以更深入地理解线性回归的原理，还能体会到使用成熟库的便利性，以及不同优化方法的特点。这些知识对于理解更复杂的机器学习算法和深度学习模型都是非常有帮助的。

希望这篇教程对你有所帮助！如果你有任何问题，欢迎在评论区留言。