### Python 中基于最小二乘法的曲线拟合实现 #### 使用 NumPy 和 Matplotlib 进行最小二乘法曲线拟合 在 Python 中，可以通过 `NumPy` 的 `polyfit` 函数轻松实现基于最小二乘法的多项式曲线拟合。该函数会返回一组系数，表示最佳拟合多项式的参数[^1]。 以下是具体实现过程： ```python import numpy as np import matplotlib.pyplot as plt # 假设我们有一些数据点 (x, y) x = np.array([0, 1, 2, 3, 4, 5]) y = np.array([0, 0.8, 0.9, 0.1, -0.8, -1]) # 定义多项式的阶数 n n = 2 # 可以调整为任意正整数 # 调用 polyfit 执行最小二乘法拟合 coefficients = np.polyfit(x, y, n) # 创建一个多项式对象以便后续绘图 polynomial = np.poly1d(coefficients) # 计算拟合后的 y 值 x_fit = np.linspace(min(x), max(x), 100) y_fit = polynomial(x_fit) # 绘制原始数据点和拟合曲线 plt.scatter(x, y, color="red", label="Data Points") # 数据点 plt.plot(x_fit, y_fit, label=f"Polynomial Fit (degree {n})") # 拟合曲线 plt.legend() plt.xlabel("X-axis") plt.ylabel("Y-axis") plt.title("Least Squares Polynomial Fitting with NumPy") plt.show() print(f"Fitted Polynomial Coefficients: {coefficients}") ``` 此代码片段展示了如何使用 `np.polyfit` 来执行最小二乘法拟合，并绘制出对应的拟合曲线。 --- #### 自定义实现最小二乘法 如果希望手动实现最小二乘法而不依赖内置库，则可以按照以下方式完成： 假设我们需要拟合一元二次多项式 \( f(x) = a_0 + a_1 x + a_2 x^2 \)，则目标是最小化残差平方和 \( S(a_0, a_1, a_2) = \sum_{i=1}^{N}(y_i - f(x_i))^2 \)[^2]。 下面是具体的自定义实现： ```python def least_squares_custom(x_data, y_data, degree): """ 自定义最小二乘法实现。 参数: x_data : 输入变量数组 y_data : 输出变量数组 degree : 多项式次数 返回: coefficients : 最佳拟合多项式的系数向量 """ X_matrix = np.vander(x_data, N=degree+1) # 构造范德蒙矩阵 A = np.dot(X_matrix.T, X_matrix) # 正规方程中的 A 部分 B = np.dot(X_matrix.T, y_data) # 正规方程中的 b 部分 coefficients = np.linalg.solve(A, B) # 解正规方程得到系数 return coefficients[::-1] # 测试自定义实现 custom_coefficients = least_squares_custom(x, y, n) custom_polynomial = np.poly1d(custom_coefficients) # 绘制结果 y_custom_fit = custom_polynomial(x_fit) plt.scatter(x, y, color="blue", label="Custom Data Points") plt.plot(x_fit, y_custom_fit, label=f"Custom Polynomial Fit (degree {n})", linestyle="--") plt.legend() plt.show() print(f"Custom Fitted Polynomial Coefficients: {custom_coefficients}") ``` 在此部分中，通过构建范德蒙矩阵并解决正规方程组的方式实现了最小二乘法[^2]。 --- #### RMS 误差评估 为了判断拟合效果的好坏，通常采用均方根误差（Root Mean Square Error, RMSE）。RMSE 表达式如下所示： \[ RMSE = \sqrt{\frac{1}{N}\sum_{i=1}^{N}(y_i - \hat{y}_i)^2} \] 其中 \( y_i \) 是真实值，\( \hat{y}_i \) 是预测值[^2]。 以下是计算 RMSE 的代码示例： ```python def calculate_rmse(y_true, y_pred): """计算均方根误差""" mse = np.mean((y_true - y_pred)**2) rmse = np.sqrt(mse) return rmse # 对比两种方法的结果 y_polyfit_pred = polynomial(x) rmse_polyfit = calculate_rmse(y, y_polyfit_pred) y_custom_pred = custom_polynomial(x) rmse_custom = calculate_rmse(y, y_custom_pred) print(f"RMSE of PolyFit Method: {rmse_polyfit:.4f}") print(f"RMSE of Custom Method: {rmse_custom:.4f}") ``` 以上代码可用于比较不同拟合方法的效果。 ---