最小二乘法进行曲线拟合

工作需求，这里记录一下数值插值和数值分析方面的算法，希望和大家一起进步。

曲线拟合的最小二乘定义

求一条曲线,使数据点均在离此曲线的上方或下方不远处,所求的曲线称为拟合曲线,
它既能反映数据的总体分布,又不至于出现局部较大的波动,更能反映被逼近函数的特性,
使求得的逼近函数与已知函数从总体上来说其偏差按某种方法度量达到最小,
这就是最小二乘法.

与函数插值不同，曲线拟合不要求曲线通过所有已知点，而是要求得到的近似函数能反映数据的基本关系，
在某种意义上,曲线拟合更有实用价值.

已知离散观测数据点f(x):(x0,y0)、(x1,y1)、(x2,y2)⋯(xn,yn),求得曲线拟合函数p(x)，记p(x)在xi处的残差为：ϵi=p(xi)−f(xi)(i=0,1,2⋯,n)记向量e=[ϵ0,ϵ1⋯,ϵn]若e的2−范数：Σi=0nϵi2=Σi=0n[p(xi)−f(xi)]2为最小则p(x)为最小二乘法拟合得来的曲线.已知离散观测数据点 f(x):(x_0,y_0)、(x_1,y_1)、(x_2,y_2)\cdots(x_n,y_n),求得曲线拟合函数p(x)，记p(x)在x_i处的残差为：\\ \epsilon_i = p(x_i) - f(x_i) (i=0,1,2\cdots,n)\\ 记向量e=[\epsilon_0,\epsilon_1\cdots,\epsilon_n]\\ 若e的2-范数：\Sigma^n_{i=0}\epsilon^2_i = \Sigma^n_{i=0}[p(x_i)-f(x_i)]^2为最小\\ 则p(x)为最小二乘法拟合得来的曲线. 已知离散观测数据点f(x):(x0,y0)、(x1,y1)、(x2,y2)⋯(xn,yn),求得曲线拟合函数p(x)，记p(x)在xi处的残差为：ϵi=p(xi)−f(xi)(i=0,1,2⋯,n)记向量e=[ϵ0,ϵ1⋯,ϵn]若e的2−范数：Σi=0nϵi2=Σi=0n[p(xi)−f(xi)]2为最小则p(x)为最小二乘法拟合得来的曲线.

最小二乘法–多项式拟合

条件：

对于给定的一组数据(xi,yi),i=1,2,...,m,寻求次数不超过n(n<<m)的多项式:y=a0+a1x+a2x2+...+anxn来拟合给定的数据，使偏差的平方和Q=Σi=1m(yi−Σj=0najxij)2为最小.对于给定的一组数据(x_i,y_i),i=1,2,...,m,寻求次数不超过n (n<<m)的多项式:\\ y = a_0 + a_1x + a_2x^2 + ... + a_nx^n 来拟合给定的数据，使偏差的平方和\\ Q = \Sigma^m_{i=1}(y_i-\Sigma^n_{j=0}a_jx^j_i)^2 为最小. 对于给定的一组数据(xi,yi),i=1,2,...,m,寻求次数不超过n(n<<m)的多项式:y=a0+a1x+a2x2+...+anxn来拟合给定的数据，使偏差的平方和Q=Σi=1m(yi−Σj=0najxij)2为最小.

Q可以看做是关于系数a的多元函数，所以拟合多项式的构造问题可归结为多元函数的极值问题：

∂Q∂ak=0,k=0,1,2,...,n即：Σi=1m(yi−Σj=0najxij)xik=0,k=0,1,2,...,n\frac {\partial Q } {\partial a_k }=0,\quad k=0,1,2,...,n\\ 即：\\ \Sigma^m_{i=1}(y_i-\Sigma^n_{j=0}a_jx^j_i)x^k_i=0 ,\quad k=0,1,2,...,n\\ ∂ak∂Q=0,k=0,1,2,...,n即：Σi=1m(yi−Σj=0najxij)xik=0,k=0,1,2,...,n

上式是关于系数a的线性方程组，有唯一解.

算法示例

多项式拟合–直线拟合

取数据点：

(1,1),(2.5,2),(3,3.3),(3.5,4),(5.6,5),(9,8),(11,9)设拟合曲线p(x)=a0+a1x(1,1),(2.5,2),(3,3.3),(3.5,4),(5.6,5),(9,8),(11,9)\\ 设拟合曲线 p(x) = a_0+a_1x (1,1),(2.5,2),(3,3.3),(3.5,4),(5.6,5),(9,8),(11,9)设拟合曲线p(x)=a0+a1x

求e -2范数关于a0,a1的偏导数，得到下列方程组：

{2Σi=1m(a0+a1xi−yi)=02Σi=1m(a0+a1xi−yi)xi=0得到：{a0m+a1Σi=1mxi=Σi=1myia0Σi=1mxi+a1Σi=1mxi2=Σi=1mxiyi\begin{cases} 2\Sigma^m_{i=1}(a_0+a_1x_i-y_i)=0\\ \\ 2\Sigma^m_{i=1}(a_0+a_1x_i-y_i)x_i=0\\ \end{cases} \\ 得到：\\ \begin{cases} a_0m+a_1\Sigma^m_{i=1}x_i=\Sigma^m_{i=1}y_i\\ \\ a_0\Sigma^m_{i=1}x_i+a_1\Sigma^m_{i=1}x^2_i=\Sigma^m_{i=1}x_iy_i\\ \end{cases} ⎩⎪⎨⎪⎧2Σi=1m(a0+a1xi−yi)=02Σi=1m(a0+a1xi−yi)xi=0得到：⎩⎪⎨⎪⎧a0m+a1Σi=1mxi=Σi=1myia0Σi=1mxi+a1Σi=1mxi2=Σi=1mxiyi

总共有7个数据点，取m=7,使用高斯消元法解上述方程组，解得：

a0=0.5466852879821776a1=0.7998090725877739a_0 = 0.5466852879821776\\ a_1 = 0.7998090725877739 a0=0.5466852879821776a1=0.7998090725877739

python图像为：

多项式拟合–抛物线拟合

设拟合曲线p(x)=ax2+bx+c设拟合曲线 p(x) = ax^2+bx+c 设拟合曲线p(x)=ax2+bx+c

求e -2范数关于a,b,c的偏导数，三个方程，利用高斯消元法，即可解得抛物线方程。
步骤和直线拟合类似，不再叙述.