用最小二乘法拟合二次曲线程序

假设有 $n$ 个点 $(x_1,y_1),(x_2,y_2),\cdots,(x_n,y_n)$，要用二次曲线 $y=a+bx+cx^2$ 来拟合这些点，求出系数 $a,b,c$。

我们可以将每个点代入二次曲线方程得到以下 $n$ 个方程：

$$
\begin{cases}
a+bx_1+cx_1^2=y_1 \\
a+bx_2+cx_2^2=y_2 \\
\cdots \\
a+bx_n+cx_n^2=y_n \\
\end{cases}
$$

将上述方程写成矩阵形式：

$$
\begin{bmatrix}
1 & x_1 & x_1^2 \\
1 & x_2 & x_2^2 \\
\vdots & \vdots & \vdots \\
1 & x_n & x_n^2 \\
\end{bmatrix}
\begin{bmatrix}
a \\
b \\
c \\
\end{bmatrix}
=
\begin{bmatrix}
y_1 \\
y_2 \\
\vdots \\
y_n \\
\end{bmatrix}
$$

设左边的矩阵为 $A$，右边的向量为 $\mathbf{y}$，未知向量为 $\mathbf{x}=(a,b,c)^T$，则以上方程可以写成 $A\mathbf{x}=\mathbf{y}$ 的形式。为了使拟合效果最好，我们需要找到一个最小二乘解 $\mathbf{x}_{ls}$，即使残差向量 $\mathbf{r}=\mathbf{y}-A\mathbf{x}_{ls}$ 的模长最小。

最小二乘解的求解公式为 $\mathbf{x}_{ls}=(A^TA)^{-1}A^T\mathbf{y}$。将系数代入二次曲线方程即可得到拟合结果。

下面是一个 Python 程序示例：

import numpy as np

def fit_quadratic_curve(x, y):
    n = len(x)
    A = np.zeros((n, 3))
    A[:, 0] = 1
    A[:, 1] = x
    A[:, 2] = x ** 2
    y = y.reshape((-1, 1))
    x_ls = np.linalg.inv(A.T @ A) @ A.T @ y
    a, b, c = x_ls.ravel()
    return a, b, c

其中，`x` 和 `y` 分别是长度为 $n$ 的一维数组，表示 $n$ 个点的横坐标和纵坐标。函数返回三个系数 $a,b,c$。