第二讲 数据拟合方法
在实验中,实验和戡测常常会产生大量的数据。为了解释这些数据或者根据这些数据做
出预测、判断,给决策者提供重要的依据。需要对测量数据进行拟合,寻找一个反映数据
变化规律的函数。数据拟合方法与数据插值方法不同,它所处理的数据量大而且不能保证
每一个数据没有误差,所以要求一个函数严格通过每一个数据点是不合理的。数据拟合方
法求拟合函数,插值方法求插值函数。这两类函数最大的不同之处是,对拟合函数不要求
它通过所给的数据点,而插值函数则必须通过每一个数据点。例如,在某化学反应中,测
得生成物的质量浓度 y (10 –3 g/cm3)与时间 t ( min )的关系如表所示
t
y
1
2
3
4
6
8
10
12
14
16
4.00
6.41
8.01
8.79
9.53
9.86
10.33
10.42
10.53
10.61
显然,连续函数关系 y(t)是
客观存在的。但是通过表中
的数据不可能确切地得到这
种关系。何况,由于仪器和
环境的影响,测量数据难免
有误差。因此只能寻求一个
近拟表达式
y =
(t)
寻求合理的近拟表达式,以
反映数据变化的规律,这种
方法就是数据拟合方法。数
据拟合需要解决两个问题:
第一,选择什么类型的函数
)(t 作为拟合函数(数学模
型);第二, 对于选定的拟合
函数,如何确定拟合函数中的参数。
11
10
9
8
7
6
5
4
0
2
4
6
8
10
12
14
16
数学模型应建立在合理假设的基础上,假设的合理性首先体现在选择某种类型的拟合
。拟合函数的选择比较灵活,可以选择线
函数使之符合数据变化的趋势(总体的变化规律)
性函数、多项式函数、指数函数、三角函数或其它函数,这应根据数据分布的趋势作出选
择。为了问题叙述的方便,将例
x3
y3
一.线性拟合(线性模型)
假设拟合函数是线性函数,即拟合函数的图形是一条平面上的直线。而表中的数据点
1 的数据表写成一般的形式
x1
y1
x4
y4
x2
y2
t
y
x5
y5
x6
y6
x7
y7
x8
y8
x9
y9
x10
y10
未能精确地落在一条直线上的原因是实验数据的误差。则下一步是确定函数
y= a + b x
中系数 a 和 b 各等于多少?从几何背景来考虑,就是要以
平面直线使得表中数据所对应的
全部落在这条直线上,如果第
线的方程,即
a 和 b 作为待定系数,确定一条
10 个点尽可能地靠近这条直线。一般来讲,数据点将不会
k 个点的数据恰好落在这条直线上,则这个点的坐标满足直
如果这个点不在直线上,则它的坐标不满足直线方程,有一个绝对值为
a
bx
k
y
k
的差
异(残差) 。于是全部点处的总误差是
a + b xk = y k
10
1k
a
bx
k
y
k
这是关于 a 和 b 的一个二元函数,合理的做法是选取
是在实际求解问题时为了操作上的方便,常常是求
a 和 b ,使得这个函数取极小值。但
a 和 b 使得函数
baF
(
,
)
10
k
1
(
a
bx
k
y
k
2)
达到极小。为了求该函数的极小值点,令
F
a
0
,
得
10
1k
(2
a
bx
k
y
k
)
0
,
F
b
10
k
1
0
,
(2
a
bx
k
y
k
)
x
k
0
这是关于未知数 a 和 b 的线性方程组。它们被称为法方程,又可以写成
10
k
1
bx
k
10
k
1
y
k
10
a
10
1
求解这个二元线性方程组便得待定系数
直线是数据的线性拟合的结果。
k
ax
k
10
k
1
2
bx
k
10
k
1
yx
k
k
a 和 b,从而得线性拟合函数 y = a + b x 。下图中
12
11
10
9
8
7
6
5
4
0
2
4
6
8
10
12
14
16
二.二次函数拟合(二次多项式模型)
假设拟合函数不是线性函数,而是一个二次多项式函数。即拟合函数的图形是一条平
面上的抛物线,而表中的数据点未能精确地落在这条抛物线上的原因是实验数据的误差。
则下一步是确定函数
y = a0 + a 1 x + a2 x 2
中系数 a0、a1 和 a2 各等于多少?从几何背景来考虑,就是要以
确定二次曲线使得表中数据所对应的
不会全部落在这条曲线上,如果第
次曲线的方程,即
10 个点尽可能地靠近这条曲线。一般来讲,数据点将
k 个点的数据恰好落在曲线上,则这个点的坐标满足二
a0、a1 和 a2 为待定系数,
如果这个点不在曲线上,则它的坐标不满足曲线方程,有一个误差(残差)
的总误差用残差平方和表示
。于是全部点处
a0 + a 1 xk + a 2 xk
2 = y k
aaaF
(
,
,
1
0
)
2
10
k
1
[(
a
0
xa
1
k
xa
2
k
2
)
2
y
k
]
这是关于 a0、a1 和 a2 的一个三元函数,合理的做法是选取
取极小值。为了求该函数的极小值点,令
a0、a1 和 a2 ,使得这个函数
得
F
0a
0
,
F
1a
0
,
F
2a
0
10
k
1
10
k
1
10
k
1
[(2
a
0
xa
1
k
xa
2
k
2
)
y
k
]
0
[(2
a
0
xa
1
k
xa
2
k
2
)
xy
k
]
k
0
[(2
a
0
xa
1
k
xa
2
k
2
)
2
xy
k
]
k
0
这是关于待定系数 a0、a1 和 a2 的线性方程组,写成等价的形式为
10
a
0
10
k
1
ax
1
k
10
k
1
2
ax
k
2
10
k
1
10
k
1
ax
k
0
2
ax
k
0
10
k
1
10
2
ax
k
1
3
ax
k
1
10
k
1
10
y
k
10
10
k
1
2
yx
k
k
3
ax
k
k
1
10
2
2
x
k
y
k
4
ax
k
k
1
k
1
k
1
这就是法方程,求解这一方程组可得二次拟合函数中的三个待定系数。下图反映了例题所
给数据的二次曲线拟合的结果
11
10
9
8
7
6
5
4
0
2
4
6
8
10
12
14
16
三. 数据的 n 次多项式拟合
x
f(x)
x2
x2
已知函数在个离散点处的函数值,假设拟合函数是
下面的函数
x1
y1
……
……
xm
ym
n 次多项式,则需要用所给数据来确定
y = a 0 + a 1 x + a 2 x 2 + …… + an x n
这里要做一个假设,即多项式的阶数
类似前面的推导,可得数据的
n 应小于题目所给数据的数目 m(例题中 m = 10)。
n 次多项式拟合中拟合函数的系数应满足的正规方程组如下
m
1
k
m
k
1
x
k
2
x
k
m
x
k
m
k
1
m
n
x
k
m
n
1
x
k
k
1
k
1
m
1
k
m
k
1
m
k
1
n
x
k
n
1
x
k
2
n
x
k
a
0
a
1
a
n
m
k
1
m
k
1
m
k
1
y
k
yx
k
k
n
x
k
y
k
从这一方程组可以看出,线性拟合方法和二次拟合方法是多项式拟合的特殊情况。从算法
上看,数据最小二乘拟合的多项式方法是解一个超定方程组
n
2
a
0
xa
11
xa
12
xa
n
1
y
1
a
0
xa
21
2
xa
2
2
n
xa
n
2
y
2
( m > n)
a
0
xa
1
m
2
xa
2
m
n
xa
mn
y
m
的最小二乘解。而多项式拟合所引出的正规方程组恰好是用超定方程组的系数矩阵的转置
矩阵去左乘超定方程组左、右两端所得。正规方程组的系数矩阵是一个病态矩阵,这类方
程组被称为病态方程组。当系数矩阵或者是右端向量有微小的误差时,可能引起方程组准
确解有很大的误差。为了避免求解这样的线性方程组,在做多项式拟合时可以将多项式中
的各次幂函数做正交化变换,使得所推出的正规方程的系数矩阵是对角矩阵。
四.点集 {x1,x2,……, xm} 上的正交多项式系
多项式 q0(x), q1(x),q2(x),……, qn(x)在点集 { x1,x2,……, xm} 上的正交
m
(
q
k
,
q
j
)
xqxq
i
(
(
)
k
j
i
)
正交多项式系可以认为是幂函数系:
正交多项式系构造的方法如下:
i
1
1,x,x 2
,……, x n
通过正交变换而得到的一组函数。
q0(x)=1,q0(x) = x – a1 ,(a1 =
m
i
1
nx
i /
),
qk(x) = (x - ak) qk -1(x) - bk qk-2(x) , ( k = 2,3,……, n)
其中,
a
k
(
xq
,
q
/()
q
,
q
)
k
1
k
1
k
1
k
1
m
1
i
m
b
k
(
q
,
q
k
1
/()
q
k
1
k
2
,
q
k
)
2
q
k
1
i
1
五.用正交多项式系组成拟合函数的多项式拟合
2
qx
i
k
1
(
x
i
/)
m
i
1
2
q
k
1
(
x
i
)
2
(
x
i
/)
m
i
1
2
q
k
2
(
x
i
)
)(
xqa
0
0
(
x
考虑拟合函数:
)
x
f(x)
中的数据代入,得超定方程
0
0
xqa
1
xqa
(
(
0
0
2
)
)
xqa
(
11
1
(
xqa
11
2
)
)
x1
y1
)(
xqa
11
x2
x2
)
)
2
2
xqa
(
2
1
(
xqa
2
2
)(
xqa
n
n
……
……
,将数据表
xm
ym
n
(
xqa
1
n
(
xqa
n
n
)
2
)
y
1
y
2
(m > n)
xqa
(
0
0
)
m
xqa
11
(
m
)
xqa
(
2
2
)
m
xqa
m
(
n
n
)
y
m
其系数矩阵为
0
(
xq
1
xq
(
0
2
)
)
)
(
xq
1
1
xq
1
(
2
)
)
(
xq
1
2
xq
2
(
2
xq
(
0
)
m
xq
(
1
)
m
xq
(
2
xq
(
n
1
xq
(
n
)
2
)
xq
(
n
)
m
)
)
m
由于多项式 q0(x),q1(x),q2(x),……, qn(x)在点集 {x1,x2,……, xm} 上的正交,所以
超定方程组的系数矩阵中不同列的列向量是相互正交的向量组。于是用这一矩阵的转置矩
阵去左乘超定方程组左、右两端得正规方程组
)
)
(
(
)
aqq
0
0
aqq
)
1
,
,
1
1
0
(
,
yq
yq
)
0
,
1
(
0
a
a
1
,
(
yq
0
yq
(
,
1
/()
/()
,
qq
0
,
qq
1
1
0
)
=>
其中,
(
qq
,
k
(
m
i
1
)
k
aqq
)
n
,
n
(
n
yq
,
n
)
a
n
(
yq
,
n
/()
qq
,
n
)
n
2
q
k
(
x
i
)
,
(
yq
k
,
)
m
i
1
yxq
(
)
k
i
。因为正规方程组中每一个方程都是
i
一元一次方程可以直接写出原超方程组的最小二乘解,所以拟合函数为
(
,
yq
n
qq
(
,
n
(
,
yq
1
qq
(
,
1
1
(
,
yq
0
qq
(
,
xq
1
xq
)
)
)(
x
)
)
)
(
)
(
0
0
0
n
)
)
xq
(
n
)
这一结果与用次多项式拟合所得结果在理论是完全一样的,只是形式上不同、算法实现上
避免了解病态方程组。
六.指数函数的数据拟合
问题 1:世界人中预测问题
下表给出了本世纪六十年代世界人口的统计数据
(单位:亿 )
年
人口
1960
29.72
1961
30.61
1962
31.51
1963
32.13
1964
32.34
1965
32.85
1966
33.56
1967
34.20
1968
34.83
有人根据表中数据,预测公元
难以置信,但现在已成为事实。试建立数学模型并根据表中数据推算出
数量。
2000 年世界人口会超过 60 亿。这一结论在六十年代末令人
2000 年世界人口的
根据马尔萨斯人口理论,人口数量按指数递增的规律发展。记人口数为
数函数 N
了计算方便, 对表达式两边取对数, 得 ln N a bt ,令
ea bt 。现需要根据六十年代的人口数据确定函数表达式中两个常数
y
ln 。于是
N
(1)计算出表中人口数据的对数值 yk = ln Nk
( k = 1, 2,……, 9)
N(t),则有指
a、 b。为
bt
a
ty )(
。
(2) 根据表中数据写出关于两个未知数
个数的方程组)
a 、b 的 9 个方程的超定方程组 (方程数多于未知数
a + b t k = y k
( k = 1, 2,……, 9)
其中, t1 =1960,t2 =1961,t3 =1962,……, t9 =1968;
y1= ln29.72,y2 = ln 30.61,……, y9 = ln34.83。
(3) 利用 MA TLAB 解线性方程组 Ax=c 的命令 A c 计算出 a 、b 的值,并写出人口增长函
数。利用人口增长函数计算出
2000 年世界人口数据: N(2000)
七.多元线性函数的数据拟合
问题 2 人的耗氧能力的数据拟合。
人的耗氧能力 y (ml/min ·kg)与下列变量有关
x1 年龄
x2 体重
x3 1.5 英里跑步所用时间
x4 静止时心速
x5 跑步时最大心速
某健身中心对 31 个自愿者进行测试,得到 31 组数据(每一组数据有 6 个数)
yk
x1k
x2k
x3k
x4k
x5k
(k = 1,… 31)
令耗氧能力为因变量,其它的指标为自变量,建立线性模型
y = a 0 + a1 x1+ a2 x2+ a3 x3+ a4 x4+ a5 x5
为了确定 6 个系数,利用已记录的数据得超定方程组
a0 + a1 x1k+ a2 x2k+ a3 x3k+ a4 x4k+ a5 x5k = y k
( k = 1, 2,……, 31)
这一方程组包含 6 个未知数 a0 ,a1, a2 , a3 ,a4 ,a5 ,但却有 31 个方程。写出超
定方程组的系数矩阵和右端向量如下
A
1
1
1
x
11
x
12
x
21
x
22
x
x
31
32
x
x
41
42
x
51
x
52
,
y
x
31,1
x
31,2
x
31,3
x
31,4
x
31,5
y
1
y
2
y
31
由最小二乘法可得正规方程组
其中, X = [a0 a1 a2 a3 a4 a5] T
T
AXA
T
yA