数据点的最小二乘多项式拟合(改进算法)

一、背景：

文章《数据点的最小二乘多项式拟合》描述了详细的多项式拟合算法，在2012-8-19之前，统计计算器也应用了该算法来实现二次函数和三次函数的拟合。为什么没有更高阶的多项式拟合功能呢？因为发现从4次开始，该算法的JavaScript计算结果非常不理想。原因是多项式阶数越高，涉及到的矩阵维度更高。在高维度的矩阵计算时，涉及到的非常多的浮点计算误差，就变得越来越不可忽略，最终导致计算结果与实际值相差甚远。于是将多项式拟合功能暂时搁置了。

在网友Goldman发信来希望在统计计算器里增加四次函数拟合功能后，我开始想有没有其他的拟合算法。最后在王宇曦同学的提醒下，了解到可以通过使用多元统计回归的途径来实现多项式拟合。下面将《数据点的最小二乘多项式拟合》文所用的方法称之为多项式拟合的普通算法，而将本文将要描述的称为多元回归算法。仔细研究后发现，同样阶数的多项式拟合，使用多元回归算法，要比普通算法节省很多运算。原因是，它用的矩阵维数更少。举例来说，4次多项式拟合，普通算法涉及到5维矩阵，而多元回归算法只涉及到4维矩阵。

2012-8-19，我将统计计算器的多项式拟合模块的算法由原来的普通算法更新为多元回归算法，Goldman的4次多项式拟合功能，于是很自然的实现了。忙完了最近的一些事，今天有空终于可以将多元回归算法记录下来，以备参考。

当然，这个算法完全不能解决JavaScript的浮点计算误差问题，只不过它使用更少的矩阵维度，得以让这个问题在低次的多项式拟合过程中，不至于被暴露出来。经过实现，统计计算器现在可以支持6阶及以下的多项式拟合。

要实现更高阶的多项式拟合，估计非要写个JavaScript的数学运算库，能够没有任何误差不可。想了想，应该需要能够直接进行分数、根号等等运算（或者说，支持符号运算）。或者直接支持符号运算。像MatLab、Mathematica等软件，是有完整的符号运算系统的，不知道有没有相应的JavaScript版本。有了解的朋友请分享给我，不甚感激。

6阶以上的多项式拟合，实用性估计也不怎么高，于是，这个功能就又暂时放下了。先把这个多元回归的多项式拟合算法记录下来吧。

二、问题：

任何算法总从问题始，这里的问题同《数据点的最小二乘多项式拟合》文里所说的。

给定一个数据点的集合：

{ (x_i, y_i) | i = 1, 2, …, n}

现在使用m次多项式函数 y = b_mx^m + b_m-1x^m-1 + … + b₂x² + b₁x + b₀ 对其进行拟合，使得拟合的误差平方和最小。

求能达到此目标的参数 [b_m , b_m-1 , …, b₂ , b₁ , b₀]^T 的值。

三、求解过程：

以4阶多项式为例。

(1) 将数据整理成下表，就像是多元统计回归中那样：

表 1：观测数据

序号 \ 变量	y	x	x2	x3	x4
1	y1	x1	x12	x13	x14
2	y2	x2	x22	x23	x24
…	…	…	…	…	…
n	yn	xn	xn2	xn3	xn4

(2) 偏微分求极小值

回归方程为，现欲求出使得

的最小值，分别对b₀, b₁, b₂, b₃, b₄求偏微分并令它们等于0：

整理得：

计算得

记： ， Error creating node with http://latex.codecogs.com/gif.latex?SS_y=\sum_{j=1}^n\left(y_i-\overline{y}\right)^2. url passed to createRemoteFileNode is either missing or not a proper web uri: http://latex.codecogs.com/gif.latex?SS_y=\sum_{j=1}^n\left(y_i-\overline{y}\right)^2，

Error creating node with http://latex.codecogs.com/gif.latex?SP_{ik}=\sum_{j=1}^n\left(x_j^i-\overline{x^i}\right)\left(x_j^k-\overline{x^k}\right)=SP_{ki}. url passed to createRemoteFileNode is either missing or not a proper web uri: http://latex.codecogs.com/gif.latex?SP_{ik}=\sum_{j=1}^n\left(x_j^i-\overline{x^i}\right)\left(x_j^k-\overline{x^k}\right)=SP_{ki}，

Error creating node with http://latex.codecogs.com/gif.latex?SP_{i0}=\sum_{j=1}^n\left(x_j^i-\overline{x^i}\right)\left(y_j-\overline{y}\right)=SP_{0i}. url passed to createRemoteFileNode is either missing or not a proper web uri: http://latex.codecogs.com/gif.latex?SP_{i0}=\sum_{j=1}^n\left(x_j^i-\overline{x^i}\right)\left(y_j-\overline{y}\right)=SP_{0i}。

将b₀代入上述方程组，可以消去一个方程，再将方程组整理成矩阵形式：

将上式简写为：

Ab = C

求得b = A^-1C。

这正是目前统计计算器多项式拟合所进行的计算。