离散试验数据点的正交多项式最小二乘拟合
数值分析--32正交多项式与最小二乘拟合
j 0 n
{ j(x) = x j } 对应代数多项式 /* algebraic polynomial */
{ j(x) = cos jx }、{ j(x) = sin jx } { j(x), j(x) }对应三 角多项式 /* trigonometric polynomial */ { j(x) = e kj x , ki kj } 对应指数多项式 /* exponential
n
定理 Ba = c 存在唯一解 0(x), 1(x), … , n(x) 线性无关。
证明：若存在一组系数 {i } 使得 0 0 + 1 1 + ... + n n 0 则等式两边分别与0, 1, … , n作内积，得到：
0 ( 0 , 0 ) + 1 ( 1 , 0 ) + ... + n ( n , 0 ) 0 ( , ) + ( , ) + ... + ( , ) 0 1 1 1 n n 1 0 0 1 . . . ( , ) + ( , ) + ... + ( , ) 0 1 1 n n n n 0 0 n
1 2 49 3 y P( x) x + x 2 10 2
cond ( B ) 7623
|| B || 484,
|| B
1
||
63 4
§6 Orthogonal Polynomials & L-S Approximation
j 例：连续型拟合中，取 j ( x) x , ( x) 1, y( x) C[0, 1]
3.4 离散数据的曲线拟合——数值分析课件PPT
4
(P0 , P0 ) i P02 (xi ) 5 i0
4
(xP0 , P0 ) i xi P02 (xi ) 2.5 i0
a0
(xP0 , P0 ) (P0 , P0 )
0.5
P1(x) x a0 x 0.5
由此得 从而有
4
(P1, P1) i P12 (xi ) 0.625 i0
aj j (x)存在唯一；
j0
(b) p *(x)
n
aj j (x)的系数
a
j
n 可由法方程组
j0
j0
(0 ,0 ) (1 ,0 )
(n ,0 )
(0 ,1 ) (0 ,n ) a0 ( f ,0 )
(1 ,1 )
(1
,n
)
a1
( f
,1
)
(n ,1 )
(n ,n )an
i1
m
xi
i 1
m
xi2
i 1
m
xi3
i 1
m
xi
2
m
yi
i1 m
xi
3
a0
a1
i 1
m
xi yi
i1 m
xi 4
a2
i1
i1
m
xi
2
yi
i1
例 3.4.1 用多项式拟合表3-4中的离散数据。
表3-4
i
1
2
3
45
xi 0.00 0.25 0.50 0.75 1.00 yi 0.10 0.35 0.81 1.09 1.96
(
f
,
n
)
或Ga
d
最小二乘法多项式拟合实验报告.docx
连续系统仿真实验报告实验数据拟合建模姓名：专业：学号：时间：2013年5月1日实验单元二实验数据拟合建模一、实验目的1、 用C 语言实现最小二乘的多项式拟合和LU 分解法；2、 熟练掌握最小二乘拟合和LU 分解法的基本原理。
3、 体会用计算机编程解决计算问题的方法。
二、需求说明（一） 、需求阐述本次实验是要求根据己知的自变量和函数值，通过多项式拟合來分别计 算2、3、4阶拟合多项式，并根据拟合结果分別计算出待求点的函数值。
其中解 拟合系数方程组时采用LU 分解的方法计算拟合多项式的系数。
（二） 、实验公式m 次拟合函数公式为：（p （x ）=ao 七1対~・・・可点"计算系数4的方程组为:Sg a 0 +S] a 】 +...4-s ni a ni =t 0 < S]a ()+s?a]+...+s mF ]a m =t]k Sm a 0 +S mH a i +• •丹加^冃 其中 》= 士疋E ，i-0所以，在编程计算时，先计算出方程组①，再用LU 分解法计算求出耳的 值，即可得到拟合多项式。
LU 分解法的公式为：其中L 矩阵和U 矩阵的计算公式如下: 第一步，当k 二1,有：「1 0 0・・・01〔21 1 0-0^31 彳32 1 ••::::0 厶L ……1-i=0n-l最后求 u nn : U nn =a nn -^l m u m r=l三、设计说明（一） 、数据结构程序采用一维数组的形式来读取文件中给出的己知点处的值和要计算的未 知点处的H 变量值，最终的拟合计算结果也是采用一维数组的形式输出到文件中。
拟合多项式的系数a 和拟合系数方程组的参数t 都是采用一维数组來存储的，而 拟合系数方程组中的参数s 和L 、U 矩阵都是用二维数组來表示的。
由于要分别 计算2、3、4阶拟合结果，所以数组的规模取为5,矩阵的规模取为5*5.（二） 、算法设计及效率分析在进行LU 分解函数中，在计算L 矩阵和U 矩阵时，因为当k=2,3.-,n-l 时， 计算丈M 和土皿的循环条件不允许k=l 时进入,而正好k=l 时，计算1“和i 町不 x-1 r-1k-1 k ・l需要工1丿匕和工1以崎，因而对k=l 和k=2,3,-,n-l,就可以和在一起计算，这样就减少了 r=lr=l程序的长度。
实验题目用正交多项式做小二乘曲线拟合
实验题目：用正交多项式做小二乘曲线拟合实验题目： 用正交多项式做最小二乘的曲线拟合 学生组号： 6 完成日期： 2011/11/27 1 实验目的针对给定数据的煤自燃监测数据中煤温与NO 22,之间的非线性关系，用正交多项式做最小二乘曲线拟合。
2 实验步骤2.1 算法原理设给定n+1个数据点：（yx kk,），k=0,1,···,n ，则根据这些节点作一个m 次的最小二乘拟合多项式pm（x ）=a+x a x a a mm x +++ (2)21=x a jmj j ∑=0①其中，m ≤n,一般远小于n.。
若要构造一组次数不超过ｍ的在给定点上正交的多项式函数系｛)(x Qj（ｊ＝０，１，．．．，ｍ）｝，则可以首先利用｛)(x Qj（ｊ＝０，１，．．．，ｍ）｝作为基函数作最小二乘曲线的拟合，即pm（x ）＝)(...)()(11x x x Qq Q q Q q mm+++ ②根据②式，其中的系数qj（ｊ＝０，１，．．．，ｍ）为∑∑===nk kjnk kjkjx Q x Q y q2)()(，ｊ＝０，１，．．．，ｍ ③将④代入③后展开就成一般的多项式。
构造给定点上的正交多项式)(x Qj（ｊ＝０，１，．．．，ｍ）的递推公式如下：⎪⎪⎩⎪⎪⎨⎧-=--=-==-+1,...,2,1),()()()()()(1)(11010m j x x x x x x x QQ Q Q Q j jj j j βαα ④其中αj＝dx x jk j=0，ｊ＝０，１，．．．，ｍ－１ ⑤βj＝dd j j1-，ｊ＝１，２，．．．，ｍ－１ ⑥∑==nk k jjx Q d2)(，ｊ＝０，１，．．．，ｍ－１ ⑦则实际计算过程中，根据⑤式逐步求出个正交多项式)(x Qj，并用公式④计算出q j，并将每次计算展开后累加到拟合多项式①中。
2.2 算法步骤用三个向量Ｂ,Ｔ,Ｓ,存放多项式)(1x Qj -，)(x Q j，)(1x Qj +的系数。
数值分析(22)离散数据的最小二乘拟合
xi
m i0
(
m
n)上至多只有n个不同零点，则称0
(
x
),
1
(
x),
...,
n
(
x
)在点集
xi
m i0
上满足Haar条件。
可以证明，如果0
(
x
),
1
(
x),
...,
n
(
x
)在
xi
m i0
上满足
Haar条件，则法方程的系数矩阵非奇异。
数值分析
数值分析
法方程的另一种形式
记 A (0 , 1 , ..., n ), 由矩阵乘法可知 G ATWA，F ATWY
则
1
A
1
x0 x1
x02 x12
x0n x1n
1
xm
xm2
xmn
m
实际计算中要计算x0n，x1n，...，xmn , xik ,(k 0,1, ..., 2n)
i 1
当数据点数值很大时，可能会在计算中溢出，可采用
P240的方法处理。
数值分析
数值分析
三、法方程的求解
求解法方程AT AC ATY或GC F (取W I )
1.确定n+1个线性无关函数
j( x)
n 的具体形式
j=0
常用的代数多项式拟合是取
j( x)＝x j (j 0,1, , n) 最佳平方逼近函数的形式是
s( x)＝c0 c1 x c2 x2 cn xn
G
(1
最小二乘拟合多项式
则法方程系数矩阵为： T 常数项为：
T a y
T
n y i i 1 n x y i i y i 1 n xm y i i i 1
其他类型的拟合问题
最小二乘法并不只限于多项式，也可用于任何具体给出的 函数形式。特别重要的是有些非线性最小二乘拟合问题通 过适当的变换可以转化为线性最小二乘问题求解。
因此，我们需要一种新的逼近原函数的办法
解决方案：
1. 2. 不要求过所有数据点（可以消除误差影响）； 尽可能地刻画数据点的趋势，靠近这些数据点。
插值与拟合的关系： 问题：给定一组数据点，构造一个函数作为近似（或逼近）。 解决方案：
1. 若要求所求曲线通过给定的所有数据点，就是插值问题； 2. 若不要求曲线通过所有数据点，而是要求它反映数据点的整体变化趋势， 这就是数据拟合，又称曲线拟合，所求出的曲线称为拟合曲线。
a0 sn a1sn 1 an s2 n un
称为正规方程组（或法方程组）。
可以证明：当 x0 , x1, 是最小值问题的解。
法方程组 可写成以 下形式:
m m xk k 1 m n xk k 1
xn , 互异时，该方程组有唯一解，并
k
x
k 1 m k 1
m
2 x k
令
x
k 1
m
n 1 k
n x y i k 1 a0 i 1 m n xy n 1 x a k 1 i 1 i i k 1 an n m 2n xin yi xk k 1 i 1
曲线拟合
可做变换
Y ln y ,
X
1 x
,
A ln a ,
B b
Y A BX 就是个线性问题
将( xi , yi ) 化为( X i ,Yi ) 后易解 A 和B
a eA , b B , P(x) a eb/x
二、 一般的最小二乘法
2
2
m
2 i
m
m
[S*(
xi
)
yi
]2
min
S ( x )
则 勒 让 德 多 项 式 为:
P~n ( x)
n! (2n)!
dn dx n
[(x2
1)n ]
勒 让 德 多 项 式 有 以 下 几个 重 要 性 质 ：
性 质1： 正 交 性
1 1
Pn
(
x)Pm
(
x)d(
x)
0, 2
2n
1
,
m n; m n.
性 质2： 奇 偶 性 Pn ( x) (1)n Pn ( x) 性 质3： 递 推 关 系 (n 1)Pn1( x) (2n 1)xPn ( x) nPn1( x) 性 质4： 在 所 有 最 高 项 系 数 为1的n次 多 项 式 中 ， 勒 让 德 多项 式
条件.
可以证明,如果0(x), 1(x),… n(x)C[a,b]在{xi}0m上满
足哈尔(Haar)条件,则法方程(5.6)的系数矩阵G非奇异.
用最小二乘法得到的法方程组(3. 6)，其系数矩阵G是
病态的,但如果0(x), 1(x),… n(x)是关于点集 {xi}(i=0,
l, ..., m)带权(xi) (i=0,l,...,m) 正交的函数族，即
课堂练习
数据拟合的最小二乘方法的实现
数据拟合的最小二乘方法的实现班级学号姓名榴莲一、实验任务给定离散样本点，采用最小二乘方法拟合样本数据，涉及的线性方程组请用高斯列主元消去法求解。
求函数f(x)=cosπx,x∈[0,1]的一次和二次最佳平方逼近多项式。
二、编程环境Windows7，Codeblock.三、算法步骤1，采用目标函数对多项式系数求偏导，得到最优值条件，组成一个方程组；2，方程组的解法采用行列式变换（两次变换：普通行列式——三角行列式——对角行列式——求解），行列式的求解算法上优化过一次了，目前还没有更好的思路再优化运算方法，限幅和精度准备再修改修改。
四、程序流程图数据结构：五、程序#include<stdio.h>#include<math.h>double qwer(double a[],double b[],int n);void guess(double a[][10],double b[],double x[],int n); int main() {double a[10][10],b[10],asd[10];int i,j,k,n;int num = 6;double y[10] = {1.0,1.004,1.031,1.117,1.223,1.422};double x1[10] = {1.0,1.0,1.0,1.0,1.0,1.0};double x2[10] = {0,0.15,0.31,0.5,0.6,0.75};n = 2;a[0][0] = qwer(x1,x1,num);a[0][1] = qwer(x1,x2,num);a[1][0] = qwer(x2,x1,num);a[1][1] = qwer(x2,x2,num);b[0] = qwer(x1,y,num);b[1] = qwer(x2,y,num);printf("矩阵A:\n");for(i=0; i<n; i++) {for(j=0; j<n; j++) printf("%-11g",a[i][j]);printf("%-11g\n",b[i]);}guess(a,b,asd,n);printf("\n三角化矩阵A:\n");for(i=0; i<n; i++) {for(j=0; j<n; j++) printf("%-11g",a[i][j]);printf("%-11g\n",b[i]);}printf("\n系数值asd：");for(i=0; i<n; i++) printf("%8.4lf",asd[i]);printf("\n一次多项式为：f(x) = %.4lf+%.4lfx\n",asd[0],asd[1]);double x3[10] = {0*0,0.15*0.15,0.31*0.31,0.5*0.5,0.6*0.6,0.75*0.75};n = 3;a[0][0] = qwer(x1,x1,num);a[0][1] = qwer(x1,x2,num);a[0][2] = qwer(x1,x3,num);a[1][0] = qwer(x2,x1,num);a[1][1] = qwer(x2,x2,num);a[1][2] = qwer(x2,x3,num);a[2][0] = qwer(x3,x1,num);a[2][1] = qwer(x3,x2,num);a[2][2] = qwer(x3,x3,num);b[0] = qwer(x1,y,num);b[1] = qwer(x2,y,num);b[2] = qwer(x3,y,num);printf("\n\n矩阵A:\n");for(i=0; i<n; i++) {for(j=0; j<n; j++) printf("%-11g",a[i][j]);printf("%-11g\n",b[i]);}guess(a,b,asd,n);printf("\n三角化矩阵A:\n");for(i=0; i<n; i++) {for(j=0; j<n; j++) printf("%-11g",a[i][j]);printf("%-11g\n",b[i]);}printf("\n系数值asd：");for(i=0; i<n; i++) printf("%8.4lf",asd[i]);printf("\n二次多项式为：f(x) = %.4lf%.4lfx+%.4lfx^2\n",asd[0],asd[1],asd[2]);/* while(1) {double cc;scanf("%lf",&cc);printf("1=%lf\n",asd[0]+asd[1]*cc+asd[2]*cc*cc);printf("2=%lf\n",0.8115+1.3418*cc-0.9262*cc*cc);}*/}double qwer(double a[],double b[],int n) {int i,j;double s = 0;for(i=0; i<n; i++) {s += a[i]*b[i];}return s;}void guess(double a[][10],double b[],double x[],int n) { int k,i,j;//guess法化下三角矩阵for(k=0; k<n-1; k++) {//找主元double ma = a[k][k];int tab = k;for(i=k+1; i<n; i++) {if(fabs(ma)<fabs(a[i][k])) {ma = a[i][k];tab = i;}}double mid;mid = b[k];b[k] = b[tab];b[tab] = mid;for(i=k; i<n; i++) {mid = a[k][i];a[k][i] = a[tab][i];a[tab][i] = mid;}//三角化for(i=k+1; i<n; i++) {b[i]=b[i]-a[i][k]/a[k][k]* b[k];for(j=k+1; j<n; j++) {a[i][j]=a[i][j]-a[i][k]/a[k][k]*a[k][j];}a[i][k]=0;}}//回代求解x值for(k=n-1; k>=0; k--) {double s =0;for(j=k+1; j<n; j++) s += a[k][j]*x[j];x[k]=(b[k]-s)/a[k][k];}}六、实验结果及分析通过数据可以看到asd: 0,9295 0.5281。
matlab求解离散型最小二乘和连续性最小二乘求解
（一）实验目的：掌握并能够利用离散型最小二乘和连续性最小二乘求解问题。
（二）问题描述问题二1.以函数生成的6个数据点（0.25,23.1）,（1.0,1.68）, （1.5,1.0）,（2.0,0.84）,（2.4,0.826）,（5.0,1.2576）为例，运行程序得到不同阶对应的曲线拟合产生的多项式函数。
2. 例题：计算f(x)=exp(x)在[-1，1]上的二、三次最佳平方逼近多项式。
并画图进行比较。
（三）算法介绍问题二的（1）是离散最小二乘问题。
最小离散最小二乘就是根据一批有误的数据（）i=1,2，…，n确定参数，并通过均方误差来比较曲线拟合的优劣，在本题中通过画图来比较不同阶方程拟合效果的优劣。
选择两种方法实现离散最小二乘。
方法一，建立normal equation(法方程组)，求解k次多项式系数。
法方程组构造方法：…方法二：由于在matlab中存在ployfit函数，可以即为方便的用k次多项式拟合。
问题二的（2）是连续最小二乘问题。
在本题中我们选用首一的Legendre多项式实现二、三次最佳平方逼近。
对于二次最佳平方逼近：对于三次最佳平方逼近：其中（四）程序问题二（1）离散最小二乘，建立normalequation.m文件。
（方法一）function f=normalequation(n)syms txx=[0.25,1.0,1.5,2.0,2.4,5.0];yy=[23.1,1.68,1.0,0.84,0.826,1.2576];x=zeros(n,n);y=zeros(n ,1);for i=1:nfor k=1:6x(i,i)=xx(k).^(2*i-2)+x(i,i);y(i,1)=yy(k).*xx(k).^(i-1)+y(i,1);endendfor i=2:nfor j=1:i-1for k=1:6x(i,j)=xx(k).^(i+j-2)+x(i,j);x(j,i)=x(i,j);endendendz=x\y;for i=1:nr=z(i,1)*t^(i-1)+r;vpa(r,5)建立number2.m文件。
用正交多项式做最小二乘拟合
由于
eijx cos( jx) i sin( jx) ( j 0,1,, N 1, i 1),
所以函数族 {1, eix ,, ei(N 1)x}在区间 [0,2π] 上是正交的.
函数 eijx 在等距点集
xk
2π k (k 0,1,, N 1) 上的值
N
eijxk 组成的向量记作
利用（5.11）中 k表达式及以上结果，得
(Pk 1, Pk 1) (xPk , Pk 1) k (Pk 1, Pk 1)
(Pk , Pk ) (Pk , Pk ) 0.
7
最后，由 k和 k的表达式（5.11）有
(Pk 1, Pk ) (xPk , Pk ) k 1(Pk , Pk ) k (Pk , Pk 1)
（5.10）
(k 1,2,, n 1).
这里 Pk (x) 是首项系数为1的 k 次多项式， 根据 Pk (x) 的 正交性，得
3
m
k 1
( xi ) xi Pk2 ( xi )
i0
m
( xi )Pk2 ( xi )
( xPk ( x), Pk ( x)) (Pk ( x), Pk ( x))
l k, l k 0; l k 0; l k;
2m
cos lx j
j0
cos kxj
2m 1 2m2 1,
l k 0; l k 0;
2m
cos lx j sin kxj 0,
0 k, j m.
j0
i0
( xPk , Pk )
最小二乘法的基本原理和多项式拟合
最小二乘法的基本原理和多项式拟合 一 最小二乘法的基本原理从整体上考虑近似函数)(x p 同所给数据点),(i i y x (i=0,1,…,m)误差i i i y x p r -=)((i=0,1,…,m)的大小，常用的方法有以下三种：一是误差i i i y x p r -=)((i=0,1,…,m)绝对值的最大值im i r ≤≤0max ，即误差 向量T m r r r r ),,(10Λ=的∞—范数；二是误差绝对值的和∑=mi ir 0，即误差向量r 的1—范数；三是误差平方和∑=mi ir02的算术平方根，即误差向量r 的2—范数；前两种方法简单、自然，但不便于微分运算 ，后一种方法相当于考虑 2—范数的平方，因此在曲线拟合中常采用误差平方和∑=mi ir02来 度量误差i r (i=0，1，…，m)的整体大小。
数据拟合的具体作法是：对给定数据 ),(i i yx (i=0,1,…，m)，在取定的函数类Φ中,求Φ∈)(x p ,使误差i i i y x p r -=)((i=0,1,…,m)的平方和最小，即∑=mi ir 02=[]∑==-mi i i y x p 02min)(从几何意义上讲，就是寻求与给定点),(i i y x (i=0,1,…,m)的距离平方和为最小的曲线)(x p y =（图6-1）。
函数)(x p 称为拟合 函数或最小二乘解，求拟合函数)(x p 的方法称为曲线拟合的最小二乘法。
在曲线拟合中，函数类Φ可有不同的选取方法.6—1二 多项式拟合假设给定数据点),(i i y x (i=0,1,…,m)，Φ为所有次数不超过)(m n n ≤的多项式构成的函数类，现求一Φ∈=∑=nk k k n x a x p 0)(,使得[]min )(00202=⎪⎭⎫⎝⎛-=-=∑∑∑===mi mi n k i k i k i i n y x a y x p I (1)当拟合函数为多项式时，称为多项式拟合，满足式（1）的)(x p n 称为最小二乘拟合多项式。
一般最小二乘拟合
第三节 一般最小二乘拟合多项式拟合形式比较规范，方法也比较简单，但在实际应用中，针对所讨论问题的特点，拟合函数可能为其他类型，如指数函数、有理函数、三角函数等，这就是一般最小二乘拟合问题。
一 线性最小二乘拟合 设)(,),(),(10x x x n ϕϕϕ 为n+1个线性无关（与向量的线性无关定义类似）的连续函数，Φ为)(),(),(10x x x n ϕϕϕ 所张成的n+1维线性空间，即由其 所有线性组合∑==∈nk k kk n k R a x a 0),,1,0(),( ϕ构成的集 合，记作{})(,),(),(10x x x span n ϕϕϕ =Φ任取Φ∈)(x p ，则∑==nk k k x a x p 0)()(ϕ，它是关于n a a a ,,,10 的线性函数。
对已知数据点),,1,0)(,(m i y x i i =，在Φ中求一)(x p ，使得[]min )()(00202=⎥⎦⎤⎢⎣⎡-=-=∑∑∑===mi mi n k i i k k i i y x a y x p I ϕ (1)这就是一般线性最小二乘拟合问题。
同多项式拟合完全类似，上述问题归结为多元函数的极值问题。
由多元函数求极值的必要条件，可得nj x y x a a Im i i j nk i i k k j ,1,0,0)())((200==-=∂∂∑∑==ϕϕ即nj y x a x x m i i i j k nk m i i k i j ,,1,0,)()()(000 ==⎥⎦⎤⎢⎣⎡∑∑∑===ϕϕϕ (2)它是关于n a a a ,,,10 的线性方程组，即为一般线性最小二乘拟合的正规方程组 或 法方程组，系数矩阵为对称矩阵。
记T m T n T m k k k k y y y y a a a a n k x x x ),,,(,),,,(,,1,0,))(,),(),((101010 ====ϕϕϕϕ则nj y x y n k j x x mi i i j j mi i k i j k j ,,1,0,)(),(,,1,0,),()(),(0====∑∑==ϕϕϕϕϕϕ式（2）可用矩阵表示为⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡)()()(),(),(),(),(),(),(),(),(),(,,1,010101110101000y y y a a a n n n n n n n n ϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕ (4)式（3）也可表示为 y G Ga GT T= （5）如果G 的列向量组线性无关，即R(G)=n+1，则正规方程组（3）或（5）存在唯一解a=n a a a ,,,10 ，从而∑==nk k k x a x p 0)()(ϕ为满足式（1）的最小二乘拟合函数。
实验数据的最小二乘拟合算法与分析
龙源期刊网
实验数据的最小二乘拟合算法与分析
作者：龙辉平习胜丰侯新华
来源：《计算技术与自动化》2008年第03期
摘要：曲线拟合的优度，以及正交多项式最佳阶数的选择是程序研制过程中一个很关键的问题，用正交多项式最小二乘法进行曲线拟合是一种处理实验数据的重要手段。
通过严格的数学推导，选择正交多项式，采用最小二乘逼近作为拟合实验数据的方法，利用F检验，确定实验数据拟合优度及最佳阶数，用编制的程序来拟合实验数据，得到了较好的结果。
关键词：正交多项式；最小二乘法；拟合；F检验；程序。
数值分析(21)离散数据的最小二乘拟合
逼近未知函数f ( x),使
m
i ( f ( xi ) s( xi ))2 min
(1)
i0
此问题称为最小二乘曲线拟合，又称为离散数据的
最佳平方逼近。
使拟合误差的平方和最小——最小二乘原理
两种拟合问题
1. 线性最小二乘曲线拟合 如：取s( x, c) cn xn cn1 xn1 c1 x c0 s( x, c)是关于系数c (c0 , c1 , , cn1 , cn )T的线性函数。 这是多项式拟合。 若取s( x, c) cne xn c2e x2 c1e x c0，这也是 关于系数的线性拟合。
若s*( x)存在,则称其为f ( x)在[a, b]上的最佳平方逼近 函数（最小二乘拟合曲线）。
在C[a, b]中，定义带权i (i 0,1, ..., m)的内积
m
f ( x),( x) i f ( xi )( xi ) Y TW （Y , ) i0
其中 W diag(0 ,1, ...,m )
线性最小二乘问题：求矛盾方程组 AC Y 的最小二乘解。
m
AC Y 2 min 2
(yi s( xi ))2 min
i0
数值分析
连续函数最佳平方逼近问题的一般提法
在内积空间C[a, b]中，设f ( x) C[a, b], 但f ( x) ,
n
在中寻找一个函数s*( x) c j j ( x)
2
a
离散数据的最佳平方逼近问题的一般提法
在内积空间C[a, b]中，设f ( x) C[a, b], 但f ( x) ,
n
在中寻找一个函数s*( x) c j j ( x) j0
使得 f ( x) s*( x) 2 min 2
基于正交多项式拟合的气象数据处理方法
基于正交多项式拟合的气象数据处理方法＊孙宝京【摘 要】针对高空气象探测数据变化规律复杂、突变情况不可预测、数据量大等特点，提出了基于离散数据正交多项式最小二乘的分段拟合方法.首先，从计算稳定性角度，说明了采用离散数据的正交多项式最小二乘法的原因；其次，从曲线拟合的保形性角度，详细阐述了二次分段拟合原则；最后，以具有典型特征的气温探测数据为例，证明了采用基于离散数据正交多项式最小二乘的分段拟合方法能够获得明显优于传统内插方法的拟合精度，提高高空气象探测数据处理的精度和自动化处理程度.【期刊名称】吉首大学学报（自然科学版）【年(卷),期】2013(034)005【总页数】4【关键词】离散数据；正交多项式；最小二乘；分段拟合；气象探测当前气象探测站整理高空气象探测数据时，主要采用的方法是将整个探测空间分层，计算规定高度的中间数值，作为计算层数据，其他数据则采用依托计算层数据内插的方法计算.上述计算方法使用高度上的数据均由计算层数据插值得到，存在一定的截断误差和舍入误差；另外在较少的计算层数据中，必然会丢弃绝大多数探测数据，造成探测资源的极大浪费.随着气象探测及作业装备自动化水平和数据处理能力的不断提高，对高空气象探测数据处理精度和时效性提出了更高的要求，这也使得从根本上改变气象参数的计算方法成为可能［1］.由于气象要素随时间、空间在不断地变化，而且探测仪器本身也受到许多复杂因素的影响，因此各种高空气象探测数据普遍具有数据量大、变化规律复杂、突变情况较多且不可预测等特点［2］.以弹道气象诸元中最重要的高空实时气温／高度曲线为例，当探测高度达到2万m时，采用电子探空仪探测的数据超过2 000组，虽然高空实时气温总体呈下降趋势，但是在局部经常出现不可预测的无规则变化.针对高空气象探测数据特点，笔者基于离散数据正交多项式最小二乘拟合思想，提出根据数据变化率的变化情况，分段拟合各种高空气象探测数据曲线的方法.借助编程软件的自动处理，经过多次试验证明，此方法可直接得出气象通报所需要的气（虚）温－高度、气压－高度、风速－高度和风向－高度函数方程（高度为自变量），并且对高空气象探测数据处理的精确性、实效性和计算简捷性都有了很大的提高.因为高空实时气温／高度曲线在弹道气象诸元中最为重要，所以笔者将以此作为说明实例.1 离散数据分段最小二乘拟合方法用插值方法所找出的近似曲线，通过所有己知点，但是由观测或实验所获得的数据，不可避免地含有误差，这就使插值所找出的近似曲线保存所有观测误差，导致所得结果可能偏离了实际情况，出现曲线拟合的过度现象.而根据实际标准尽可能地通过节点的函数值近旁，以部分地抵消原始数据所包含的观测误差，则更具有实用价值，最小二乘拟合就是其中的一个重要方法［3］.1.1 离散数据的正交多项式最小二乘拟合在线性空间由Φ＝span（φ0，φ1，...，φn）求解最小二乘曲线拟合问题，经常会出现法方程组是病态方程组的情况，特别是当样本组数n偏大时更是如此，同时在求解法方程组时系数矩阵或右端项微小的扰动都可能导致解函数有很大的误差.为了确保计算过程的稳定性和精确性，笔者利用基于高空气象探测数据的正交多项式作为拟合函数系，使法方程组的系数矩阵变为对角矩阵，从而解得正交多项式的系数［4－5］.根据点集｛x1，x2，...，xm｝，构造出相应的正交多项式系｛p0（x），p1（x），...，pn（x）｝.由正交多项式三项递推关系，最高次项系数为1的正交多项式系｛pk（x）｝（k＝0，1，2，...，n）有如下递推关系：其中pk（x）为最高项系数为1的k次多项式.由｛pk（x）｝正交性可知，1.2 二次分段拟合原则因一次完整综合气象探测的数据量大，若对整组数据进行一次性拟合，则拟合精度和计算效果都比较差.同时，为克服传统分段曲线拟合方法中对数据点分段时经验成分较多的缺点，笔者根据数据曲线斜率的变化情况对曲线进行分段拟合，并且加入了使相邻曲线连续，即曲线边界点必须连续的约束条件［6－7］.考虑到高度是单调递增，为了适应高空气象数据曲线在局部随高度变化的无规则性，采取二次分段原则.一次分段：在气温和高度曲线中，设Δti＝（ti＋1－ti）／（hi＋1－hi），其中ti和hi分别为探测数据中第i个探测点的气温和高度，ti＋1和hi＋1分别为探测数据中第i＋1个探测点的气温和高度.若Δti－1为负，Δti，Δti＋1，Δti＋2，...，Δti＋a（a≥1）为正，Δti＋a＋1，...，Δti＋a＋b（b≥1）为负，Δti＋a＋b＋1为正，则｛ti，...，ti＋a＋b｝为第i个分段，第i＋1个分段从点ti＋a＋b开始.二次分段：在气温和高度曲线中，第i个分段已经由一次分段确定，设δi＋c＝Δti＋c－Δti＋c－1（1≤c≤a＋b－1）.若δi＋1，δi＋2，...，δi＋d（d≥1）均大于标准值δ0（为负数），δi＋d＋1≤δ0，且δi＋d＋2，...，δi＋d＋e（e≥1）小于或等于（－δ0），δi＋d＋e＋1＞（－δ0），则ti＋d＋e＋1为分段点，即｛ti，...，ti＋d＋e＋1｝为一个分段，下一个分段以ti＋d＋e＋1为起点，从δi＋d＋e＋2开始判断；若δi＋1≤δ0，且δi＋2，...，δi＋e（e≥1）小于或等于（－δ0），δi＋e＋1＞（－δ0），则ti＋e＋1为分段点，即｛ti，...，ti＋e＋1｝为一个分段，下一个分段以ti＋e＋1为起点，从δi＋e＋2开始判断［8－9］.为使在单个分段内拟合曲线在误差允许范围内最大程度地保持数据曲线的形状，并且减少数据处理程序的计算量，实际计算过程对原始探测数据中心化计算后，采用如下做法进行拟合：若a＋b≥5，则用四阶正交多项式为拟合函数系，即当一个分段包含4个以上的点时用四阶正交多项式拟合，若a＋b＝4，则用三阶正交多项式为拟合函数系，即当一个分段包含4个点时用三阶正交多项式拟合，若a＋b＝3，则用二阶正交多项式为拟合函数系，即当一个分段包含3个点时用二阶正交多项式拟合，同时，使参与第i个分段拟合的第i个分段中的最后一个数据点，作为第i＋1个分段的第1个数据点参与第i＋1个分段的拟合，这样第i个分段的拟合函数和第i＋1个分段的拟合函数在分界点上的值相等，从而实现所有分段拟合函数在整体上保持连续［10］，即2 计算实例以编号为087233的GZZ8型电子探空仪探测数据为例.图1是该组探测数据曲线中一部分，温度t∈［－61.31，－62.27］，高度H∈［15 125.89，1 5 494.56］.根据文中的分段规则，数据处理系统自动将此段数据曲线分为5段进行拟合.第1个分段为前13个点，h∈［15 125.89，15 230.23］，拟合函数为平方误差为0.08.第2个分段为后续7个点，h∈［15 230.23，15 289.47］，拟合函数为平方误差为0.03.第3个分段为后续10个点，h∈［15 289.47，15 397.85］，拟合函数为平方误差为0.02.第4个分段为后续5个点，h∈［15 397.85，15 439.32］，拟合函数为平方误差为0.01.第5个分段为最后7个点，h∈［15 439.32，15 494.56］，拟合函数为平方误差为0.03.根据传统高空气象探测数据处理方法，14 000～16 000m高度范围内的相关气象数据从标准高度层14 000m和16 000m处的实际探测数据进行线性内插得到.在此例中，14 000m和16 000m处的气温分别为－56.16℃和－62.37℃，从而内插得到高度为15 142.00，15 276.4，15 378.67，15 417.06，1 5 479.1m处的气温.各高度对应的实际探测值、内插结果和拟合结果见表1.由表1可见，采用基于离散数据正交多项式最小二乘的分段拟合法，可以获得比传统内插法更好的拟合精度.3 结语针对高空气象探测数据量大、变化规律复杂、突变情况较多且不可预测等特点，提出了基于离散数据正交多项式最小二乘的分段拟合法.经过实践证明，该方法能够根据算法自动分段，拟合出符合精度要求的曲线，在处理程序自动计算过程中确保计算的稳定性和高效性，避免传统高空气象探测数据处理方法中不断累计的舍入误差对结果的影响，提高高空气象探测数据处理结果的精度和自动化处理程度.参考文献：［1］ 孙宝京.炮兵防空兵气象信息系统研究［D］.沈阳：东北大学，2004.［2］ 曲延禄.外弹道气象学概论［M］.北京：气象出版社，1987：59－176.［3］ FRANK R.Scattered Data Interpolation：Tests of Some Methods［J］.Mathematical Computation，1982，38：181－200.［4］ MICCHELLI C.Interpolation of Scattered Data：Distance Matrix and Conditionally Positive Definite Functions［J］.Constructive Approximation，1986，2：11－22.［5］ 蔺小林，蒋耀林.现代数值分析［M］.北京：国防工业出版社，2004：212－255.［6］ WU Zong-min，SCHABACK R.Shape Preserving Properties and Convergence of Univariate Multi-Quadric Quasti-Interpolation［J］.ACTA Mathematice Application Sinica，1994，10：441－446.［7］ GREGORY J A.Shape Preserving Spline Interpolation［J］.Computer-Aided Design，1986，18（1）：53－57.［8］ KAUFMANN E，KLASS R.Smoothing Surfaces Using Reflection Lines for Families of Splines［J］.Computer-Aided Design，1988，20（6）：312－316.［9］ HOSCHEK J.Smoothing Curves and Surfaces［J］.Computer-Aided Geometric Design，1985，2（1－3）：97－105.［10］ JORG PRTERS.Local Smooth Surface Interpolation：A Classification［J］.Computer-Aided Geometric Design，1 990（7）：191－195.（责任编辑 向阳洁）。
最小二乘法的本原理和多项式拟合
第一节最小二乘法的基本原理和多项式拟合 最小二乘法的基本原理从整体上考虑近似函数P(x )同所给数据点(X i ，y i )(i=o,l,…,m)误差riP (x i ) y i (i=o,i,…,m)的大小，常用的方法有以下三种：一是误差rip(xi) y i (i=o,i,…,m)绝对值的最大值maXr ^，即误差 向量mr (ro，r1，咕厂的范数；二是误差绝对值的和i o *，即误差向量r 的1 —m2 ri范数；三是误差平方和i o 的算术平方根，即误差向量r 的2—范数；前两种 方法简单、自然，但不便于微分运算，后一种方法相当于考虑2—范数的平方, m 2 ri因此在曲线拟合中常采用误差平方和io 来 度量误差ri (i=O , 1,…，m)的整 体大小。
数据拟合的具体作法是：对给定数据 (Xi,yi ) (i=o,i,…，m)，在取定的函数类 中,求P (x ),使误差r iP (X i ) y i (i=o,i,…,m)的平方和最小，即 mm 2 /、 2r iP (xJ y i mini oi o 拟合多项式。
特别地，当n=1时，称为线性拟合或直线拟合。
P n (x)成的函数类，现求一mIP n (X i )i oka k Xk o \ m2y ii onk a k Xiy i当拟合函数为多项式时，称为多项式拟合，满足式(min(1)1) 的 Pn(x)称为最小二乘从几何意义上讲，就是寻求与给定点(xi,%)(i=O,1,…,m)的距离平方和为最 小的曲线 y P (x )(图6-1 )。
函数P (x )称为拟合 函数或最小二乘解，求拟 合函数P (x )的方法称为曲线拟合的最小二乘法。
在曲线拟合中，函数类可有不同的选取方法•多项式拟合假设给定数据点(x i ,y i )(i=O,1,…,m)，为所有次数不超过n(n m)的多项式构n显然 I ( a k X i k y i )2i 0 k 0 为a o , a i , a n 的多元函数，因此上述问题即为求 由多元函数求极值的必要条件，得 m n—2 ( a k X k yj 彳 0,a j i 0 k 0即 n mmj k j(X i )a kX i y i ,k 0 i 0i 01I(a0,a1, an)的极值问题。