R r1 ( xn ) rm ( xn ) nm
1
1
m
1
我们一般选取多项式做拟合：f(x)=a1xm+ …+amx+am+1，即 取{1,x,x2,·,xm-1}。此时矩阵R含有一个m阶子式是范德蒙行 · · 列式，从而有惟一解。
用MATLAB做多项式最小二乘拟合
数学建模的基本方法
----数据处理和拟合模型
拟合模型原理
如果想确定具有因果关系的变量之间的确 定函数关系，可通过先测量一组数据，再 通过数学方法得到具体的函数表达式模型， 利用该模型可进行解释所研究的问题，可 进行预测。这个过程就称为拟合过程，得 到的模型称为拟合模型。
拟 合 模 型 实 例 1
可化简为一次线性的非线性最小二乘法 • 10. y=a+b1f1(x)+b2 f2(x)+…+bn fn(x) • 令 ui= fi(x), 则有 • y=a+b1u1+…+bnun. • 20. y=a ebx . • 令 z=ln y, 则有 • z = ln a + b x = a* + b x . • 3 0. y = a x b . • 令 z = ln y, u = ln x, 则有 • z = ln y = lnａ＋b ln x = a*+ b u
k 1
m
(2)
问题归结为，求 a1,a2, …am 使 J(a1,a2, …am) 最小。
线性最小二乘法的求解－超矩阵解法
m n 1 r ( xi )[ ak rk ( xi ) yi ] 0 J k 1 i 1 0 (3) a j n m ( j 1, m) rm ( xi )[ ak rk ( xi ) yi ] 0 i 1 k 1
上机作业：下面是美国黄松的数据：其中x表示树身 中部测得的直径（单位：英寸）；y是体积的度量。
x 17 19 20 22 23 25 28
y
x y
19
31 140
25
32 153
32
33 187
51
36 192
57
37 205
71
39 250
113
42 260
请按下面给出的函数类型用最小二乘法则进行参数估计，并 判断优劣（使用计算机）：
1）输入以下命令： x=0.1:0.1:1.1; y=[-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2]; R=[(x.^2)’ x’ ones(11,1)];％第三列全是1 A=R\y’ 2）计算结果： Ａ = -9.8108 20.1293 -0.0317
+ i (x+ i) i,y
+
+
+
+
y=f(x)
x i =|yi－f(xi)|为点（xi,yi) 与曲线 y=f(x) 的距离
建立拟合模型需解决的问题 •如何选定拟合的函数类型？例如用多项式型函数 还是指数型函数去拟合？ •在某个拟合的准则意义下如何在一类函数（带有 参数的函数族）中选择最佳的函数？即：从该类 函数中选出最佳的函数（确定函数中的具体参 数），使之在此准则的意义下最精确地代表了数 据。 •如何从一些已经拟合好的类型中选择最合适的？ 例如判断最佳的指数型函数是否比最佳的多项式 型函数更合适？
Min ∑ | yi－f(xi) |2
绝对偏差的平方和是优化的指标； 优化标准是指标越小越好。 问题转化求使得该指标最小的参数值，即求优 化指标是参数函数的最小值点的问题。
拟合模型最常用的准则——线性最小二乘法的基本思路
先选定一组函数 r1(x), r2(x), …rm(x), m<n, 令 注意：此优 f(x)=a1r1(x)+a2r2(x)+ …+amrm(x) 化函数是以 其中 a1,a2, …am 为待定系数。 参数为自变 量的函数。 确定a ,a , …a 的准则（最小二乘准则）：
1 2 m
（1）
函数类型是这些 使n个已知的数据点（xi,yi) 与曲线 y=f(x) 的距离i 的平方和 最小 。 r(x)的线性组合 n n 记 J (a , a , a ) 2 [ f ( x ) y ]2
1 2 m

i 1 n i 1
i

i 1
i
i
[ ak rk ( xi ) yi ]2
7
6.5
6
5.5
5
4.5 二次拟合 三次拟合 六次拟合 4
3.5 0
1
2
3
4
5
6
7
8
9
10
为了准确判断拟合效果，需计算“节点处的总 误差”： （续前面程序） wc1=sqrt(sum((polyval(a1,x0)-y0).^2)) wc2=sqrt(sum((polyval(a2,x0)-y0).^2)) wc3=sqrt(sum((polyval(a3,x0)-y0).^2)) wc6=sqrt(sum((polyval(a6,x0)-y0).^2)) 执行得： wc1 =0.4188 wc2 =0.0565 wc3 =0.0078 wc6 =0.000705
y=ax+b; y=ax2+b y=ax2+bx+c; y=ax2+bx, y=axb;
已知一室模型快速静脉注射下的血药浓度数据(t=0注射300mg) t (h) 0.25 0.5 1 1.5 2 3 4 6 8
c (g/ml) 19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01
求血药浓度随时间的变化规律c(t).
10
2
c (t ) c0 e kt
clear hold on x0=1:9;y0=[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50]; for i=1:9 plot(x0(i),y0(i),'+') end a1=polyfit(x0,y0,1),a2=polyfit(x0,y0,2) ， a3=polyfit(x0,y0,3),a6=polyfit(x0,y0,6) x=0:0.1:10; y1=polyval(a1,x);y2=polyval(a2,x); y3=polyval(a3,x);y6=polyval(a6,x);plot(x,y1,x,y2,'g',x,y3,'b',x,y6,'r') hold off gtext('二次拟合'),gtext('三次拟合'),gtext('六次拟合')
-0.447 1.978
9.30 11.2
即要求 出二次多项式:
f ( x) a1x 2 a2 x a3
中 的 A (a1 , a2 , a3 ) 使得:
[ f ( xi ) yi ]2
i 1
11
最小
解法1．用解超定方程的方法
此时 x12 R 2 x11 x1 x11 1 1
f ( x) 9.8108x 2 20.1293x 0.0317
12
解法2：用多项式拟合的命令
10 8 6
1）输入以下命令： x=0:0.1:1; y=[-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2]; A=polyfit(x,y,2) z=polyval(A,x); plot(x,y,‘k+’,x,z,‘r’) %作出数据点和拟合曲线的图 形“k+”指数据点(x,y)是＋黑色；“r”指线是红色
温度t(0C) 20.5 32.7 51.0 73.0 95.7 已知热敏电阻数据： 电阻R() 765 求600C时的电阻R。
1100 1000 900 800 700 20
826
873
942 1032
确定函数类型为 R=at+b a,b为待定系数， 如何估计？
40 60 80 100
拟 合 模 型 实 例 2
温度t(0C) 20.5 32.7 51.0 73.0 95.7 电阻R() 765 826 873 942 1032
拟合R=a1t+a2
用命令 polyfit(x,y,m) 得到 a1=3.3940, a2=702.4918
例1： 对下面一组数据作二次多项式拟合
xi yi 0.1 0.2 0.3 3.28 0.4 6.16 0.5 7.08 0.6 7.34 0.7 7.66 0.8 9.56 0.9 9.48 1.0 1.1
4 2 0 -2 0 0.2 0.4 0.6 0.8 1
2）计算结果： Ａ = -9.8108
20.1293
-0.0317
f ( x) 9.8108x 2 20.1293x 0.0317
例2：对函数C=C(t)测量得下面一组数据： t ： 1 2 3 4 5 6 7 8 9 C：4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 试分别用1次、2次、3次、6次多项式作拟合，并画图显示 拟合效果。
T T
(4) 超矩阵解法
称为正规方程组或法方程组.当RTR可逆时，（4）有唯一解：
a ( RT R) 1 RT y
(5)
线性最小二乘拟合 f(x)=a1r1(x)+ …+amrm(x)中 函数{r1(x), …rm(x)}的选取
怎样选择{r1(x), …rm(x)}，以保证系数{a1,…am}有唯一解？ 提示 {a1,…am}有唯一解 RTR可逆 Rank(RTR)=m Rank(R)=m R列满秩 列向量组{r1(x), …rm(x)} r ( x ) r ( x ) 线性无关