编写一个程序，可实现任意有限阶的多项式最小二乘曲线拟合功能。

编写一个程序用于实现任意有限阶的多项式最小二乘曲线拟合，通常需要几个关键步骤：

1. **数据准备**：首先，你需要输入一组(x, y)坐标对的数据，表示你想要拟合的点。

2. **选择多项式阶数**：确定你希望拟合的多项式的最高次数。比如，如果你选择了3次多项式，那么函数形式可能是y = ax^3 + bx^2 + cx + d。

3. **矩阵构造**：构建一个增广矩阵，它包含x的所有幂作为列，最后一列则是对应的y值。例如，对于3次多项式，矩阵会是 [[x1, x1^2, x1^3, 1], [x2, x2^2, x2^3, 1], ..., [xn, xn^2, xn^3, 1]]。

4. **求解线性方程组**：将增广矩阵转置并乘以其自身（也就是计算其自乘积），再加上一个常数项的矩阵I（单位矩阵），形成一个关于系数a、b、c、d的齐次线性方程组。然后，通过高斯消元法或其他线性代数方法求解这个方程组。

5. **计算残差**：用得到的系数构建多项式，并计算每个点的残差（预测值减去实际值）。这可以帮助评估拟合效果。

6. **结果输出**：输出多项式表达式以及它的系数，以及一个图表显示原始数据点和拟合曲线，直观地展示拟合效果。

下面是简单的Python代码示例（假设使用numpy库）：

import numpy as np

def poly_fit(data_points, degree):
    x, y = zip(*data_points)
    X = np.vander(x, degree + 1)  # Vandermonde matrix
    coefficients = np.linalg.lstsq(X, y)[0]  # Least squares solution
    return coefficients, np.polyval(coefficients, x)

# 使用示例
data_points = [(0, 0), (1, 1), (2, 4)]
degree = 2
coeffs, fitted_values = poly_fit(data_points, degree)
print(f"Polynomial coefficients: {coeffs}")