2023年8月3日发(作者:)

西北师大数学系

基于最小二乘法的多项式拟合

一 最小二乘法的基本原理设已知某物理过程

yf(x)的一组观测数据

(xi,f(xi)),

i1,2,,m. (1)

要求在某特定函数类(x)寻求一个函数(x)作为yf(x)的近似函数,使得二者在xi上的残差

i(xi)f(xi),

i1,2,,m (2)

按某种度量标准为最小,这就是拟合问题.

[0,1,,m]Tii要求残差按某种度量标准为最小,即要求由残差构成的残差向量的某种范数为最小,要求1,或m即

1i(xi)f(xi)i0i0m

maximax(xi)f(xi)ii

为最小,这本来都是很自然的,可是计算不太方便.通常要求:

或者

2(){[(xi)f(xi)]}2i2i0i0m12m12

22[(xi)f(xi)]22ii0i0mm (3)

为最小.这种要求误差平方和最小的拟合称为曲线拟合的最小二乘法.就是说,最小二乘法提供了一种数学方法,利用这种方法可以对实验数据实现在最小平方误差意义下的最好拟合.在曲线拟合中,函数类可有不同的选取方法.下面就常用的多项式拟合做介绍。

二 多项式拟合

(x,y)假设给定数据点ii(i=0,1,…,m),为所有次数不超过n(nm)的多项式构成的函数类,现求一pn(x)akxkk0n,使得

m2n2Ipn(xi)yiakxikyimini0i0k0 (4)

p(x)称为最小二乘拟合多项式。特别地,当n=1当拟合函数为多项式时,称为多项式拟合,满足式(4)的n时,称为线性拟合或直线拟合。显然

mnmI(akxikyi)2i0k0 西北师大数学系

n的多元函数,因此上述问题即为求为01必要条件,得

a,a,aII(a0,a1,an)的极值 问题。由多元函数求极值的mnI2(akxikyi)xij0,aji0k0j0,1,,n (5)

(xk0i0nmjki)akxijyi,i0mj0,1,,n (6)

(6)是关于a0,a1,an的线性方程组,用矩阵表示为

m1mxii0mxini0xxi0i0mmi2in1ixi0mmxyiai0i00mmn1axi1xiyii0i0ammn2nnxiyixii0 (7)

i0nim式(6)或式(7)称为正规方程组或法方程组。

可以证明,方程组(7)的系数矩阵是一个对称正定矩阵,故存在唯一解。从式(7)中解出n),从而可得多项式

ak(k=0,1,…,pn(x)akxkk0n (8)

可以证明,式(8)中的为最小二乘拟合多项式pn(x)满足式(4),即pn(x)为所求的拟合多项式。我们把i0pn(x)的平方误差,记作

222pmn(xi)yi2称r由式(5)可得

mpn(xi)yii0nmm

r22yak(xikyi)2ii0k0i0 (9)

多项式拟合的一般方法可归纳为以下几步:

(1) 由已知数据画出函数粗略的图形——散点图,确定拟合多项式的次数n;

(2) 列表计算xi0mji(j0,1,,2n)和xi0mjiyi(j0,1,,2n);

(3) 写出正规方程组,求出a0,a1,an;nk0

(4) 写出拟合多项式在实际应用中,npn(x)akxk。

m或nm;当nm时所得的拟合多项式就是拉格朗日或牛顿插值多项式。R()如表6-1,求电阻R与温度 T的近似函数关系。T例1 测得铜导线在温度i(℃)时的电阻ii 0

19.1

1

25.0

2

30.1

3

36.0

4

40.0

5

45.1

6

Ti(℃)

50.0 西北师大数学系

Ri()

76.30 77.80 79.25 80.80 82.35 83.90 85.10

解 画出散点图(图6-2),可见测得的数据接近一条直线,故取n=1,拟合函数为

Ra0a1T列表如下

i

0

1

2

3

4

5

6

Ti

19.1

25.0

30.1

36.0

40.0

45.1

50.0

245.3

Ri

76.30

77.80

79.25

80.80

82.35

83.90

85.10

565.5

Ti2

364.81

625.00

906.01

1296.00

1600.00

2034.01

2500.00

9325.83

TiRi

1457.330

1945.000

2385.425

2908.800

3294.000

3783.890

4255.000

20029.445

正规方程组为

解方程组得

245.3a0565.57245.39325.83a20029.4451

a070.572,故得R与T的拟合直线为

利用上述关系式,可以预测不同温度时铜导线的电阻值。例如,由R=0得T=-242.5,即预测温度时,铜导线无电阻。

a10.921

T=-242.5℃R70.5720.921T

6-2

例2 已知实验数据如下表

i 0

1

10

1

3

5

2

4

4

3

5

2

4

6

1

5

7

1

6

8

2

7

9

3

8

10

4

xi

yi

试用最小二乘法求它的二次拟合多项式。解 设拟合曲线方程为

ya0a1xa2x2

列表如下 西北师大数学系

I

0

1

2

3

4

5

6

7

8

xi

1

3

4

5

6

7

8

9

10

53

yi

10

5

4

2

1

1

2

3

4

32

xi2

1

9

16

25

36

49

64

81

100

381

xi3

1

27

64

125

216

343

512

729

1000

3017

xi4

1

81

256

625

1296

2401

4096

6561

10000

25317

xiyi

10

15

16

10

6

7

16

27

40

147

xi2yi

10

45

64

50

36

49

128

243

400

1025

得正规方程组

解得

52381a0329523813017a1471381301725317a21025

a013.4597,故拟合多项式为

a13.6053a20.2676

y13.45973.60530.2676x2

三 多项式拟合的MATLAB实现

用Polyfit函数P=polyfit(x,y,n)对数据进行拟合,返回n次多项式的系数,并用降序排列的向量表示,长度为n+1.

p(x)p1xnp2xn1pnxpn1

[p,s]=polyfit(x,y,n)返回多项式系数向量p和矩阵s。s与polyval函数一起用时,可以得到预测值的误差估计。

例如在MATLAB界面中输入一下命令

>> x=[0 0.0385 0.0963 0.1925 0.2888 0.385];

>> y=[0.042 0.104 0.186 0.338 0.479 0.612];

>> [p,s,mu]=polyfit(x,y,5)

输出结果为:

p =

Columns 1 through 5

0.0193 -0.0110 -0.0430 0.0073 0.2449

Column 6

0.2961

说明拟合的多项式为:0.0193x0.0110x0.043x0.0073xs =

R: [6x6 double]

df: 0

normr: 2.3684e-016

mu =

0.1669

0.1499

自由度为 0

标准偏差为 2.3684e-016

54320.2449x0.2961

2023年8月3日发(作者:)

西北师大数学系

基于最小二乘法的多项式拟合

一 最小二乘法的基本原理设已知某物理过程

yf(x)的一组观测数据

(xi,f(xi)),

i1,2,,m. (1)

要求在某特定函数类(x)寻求一个函数(x)作为yf(x)的近似函数,使得二者在xi上的残差

i(xi)f(xi),

i1,2,,m (2)

按某种度量标准为最小,这就是拟合问题.

[0,1,,m]Tii要求残差按某种度量标准为最小,即要求由残差构成的残差向量的某种范数为最小,要求1,或m即

1i(xi)f(xi)i0i0m

maximax(xi)f(xi)ii

为最小,这本来都是很自然的,可是计算不太方便.通常要求:

或者

2(){[(xi)f(xi)]}2i2i0i0m12m12

22[(xi)f(xi)]22ii0i0mm (3)

为最小.这种要求误差平方和最小的拟合称为曲线拟合的最小二乘法.就是说,最小二乘法提供了一种数学方法,利用这种方法可以对实验数据实现在最小平方误差意义下的最好拟合.在曲线拟合中,函数类可有不同的选取方法.下面就常用的多项式拟合做介绍。

二 多项式拟合

(x,y)假设给定数据点ii(i=0,1,…,m),为所有次数不超过n(nm)的多项式构成的函数类,现求一pn(x)akxkk0n,使得

m2n2Ipn(xi)yiakxikyimini0i0k0 (4)

p(x)称为最小二乘拟合多项式。特别地,当n=1当拟合函数为多项式时,称为多项式拟合,满足式(4)的n时,称为线性拟合或直线拟合。显然

mnmI(akxikyi)2i0k0 西北师大数学系

n的多元函数,因此上述问题即为求为01必要条件,得

a,a,aII(a0,a1,an)的极值 问题。由多元函数求极值的mnI2(akxikyi)xij0,aji0k0j0,1,,n (5)

(xk0i0nmjki)akxijyi,i0mj0,1,,n (6)

(6)是关于a0,a1,an的线性方程组,用矩阵表示为

m1mxii0mxini0xxi0i0mmi2in1ixi0mmxyiai0i00mmn1axi1xiyii0i0ammn2nnxiyixii0 (7)

i0nim式(6)或式(7)称为正规方程组或法方程组。

可以证明,方程组(7)的系数矩阵是一个对称正定矩阵,故存在唯一解。从式(7)中解出n),从而可得多项式

ak(k=0,1,…,pn(x)akxkk0n (8)

可以证明,式(8)中的为最小二乘拟合多项式pn(x)满足式(4),即pn(x)为所求的拟合多项式。我们把i0pn(x)的平方误差,记作

222pmn(xi)yi2称r由式(5)可得

mpn(xi)yii0nmm

r22yak(xikyi)2ii0k0i0 (9)

多项式拟合的一般方法可归纳为以下几步:

(1) 由已知数据画出函数粗略的图形——散点图,确定拟合多项式的次数n;

(2) 列表计算xi0mji(j0,1,,2n)和xi0mjiyi(j0,1,,2n);

(3) 写出正规方程组,求出a0,a1,an;nk0

(4) 写出拟合多项式在实际应用中,npn(x)akxk。

m或nm;当nm时所得的拟合多项式就是拉格朗日或牛顿插值多项式。R()如表6-1,求电阻R与温度 T的近似函数关系。T例1 测得铜导线在温度i(℃)时的电阻ii 0

19.1

1

25.0

2

30.1

3

36.0

4

40.0

5

45.1

6

Ti(℃)

50.0 西北师大数学系

Ri()

76.30 77.80 79.25 80.80 82.35 83.90 85.10

解 画出散点图(图6-2),可见测得的数据接近一条直线,故取n=1,拟合函数为

Ra0a1T列表如下

i

0

1

2

3

4

5

6

Ti

19.1

25.0

30.1

36.0

40.0

45.1

50.0

245.3

Ri

76.30

77.80

79.25

80.80

82.35

83.90

85.10

565.5

Ti2

364.81

625.00

906.01

1296.00

1600.00

2034.01

2500.00

9325.83

TiRi

1457.330

1945.000

2385.425

2908.800

3294.000

3783.890

4255.000

20029.445

正规方程组为

解方程组得

245.3a0565.57245.39325.83a20029.4451

a070.572,故得R与T的拟合直线为

利用上述关系式,可以预测不同温度时铜导线的电阻值。例如,由R=0得T=-242.5,即预测温度时,铜导线无电阻。

a10.921

T=-242.5℃R70.5720.921T

6-2

例2 已知实验数据如下表

i 0

1

10

1

3

5

2

4

4

3

5

2

4

6

1

5

7

1

6

8

2

7

9

3

8

10

4

xi

yi

试用最小二乘法求它的二次拟合多项式。解 设拟合曲线方程为

ya0a1xa2x2

列表如下 西北师大数学系

I

0

1

2

3

4

5

6

7

8

xi

1

3

4

5

6

7

8

9

10

53

yi

10

5

4

2

1

1

2

3

4

32

xi2

1

9

16

25

36

49

64

81

100

381

xi3

1

27

64

125

216

343

512

729

1000

3017

xi4

1

81

256

625

1296

2401

4096

6561

10000

25317

xiyi

10

15

16

10

6

7

16

27

40

147

xi2yi

10

45

64

50

36

49

128

243

400

1025

得正规方程组

解得

52381a0329523813017a1471381301725317a21025

a013.4597,故拟合多项式为

a13.6053a20.2676

y13.45973.60530.2676x2

三 多项式拟合的MATLAB实现

用Polyfit函数P=polyfit(x,y,n)对数据进行拟合,返回n次多项式的系数,并用降序排列的向量表示,长度为n+1.

p(x)p1xnp2xn1pnxpn1

[p,s]=polyfit(x,y,n)返回多项式系数向量p和矩阵s。s与polyval函数一起用时,可以得到预测值的误差估计。

例如在MATLAB界面中输入一下命令

>> x=[0 0.0385 0.0963 0.1925 0.2888 0.385];

>> y=[0.042 0.104 0.186 0.338 0.479 0.612];

>> [p,s,mu]=polyfit(x,y,5)

输出结果为:

p =

Columns 1 through 5

0.0193 -0.0110 -0.0430 0.0073 0.2449

Column 6

0.2961

说明拟合的多项式为:0.0193x0.0110x0.043x0.0073xs =

R: [6x6 double]

df: 0

normr: 2.3684e-016

mu =

0.1669

0.1499

自由度为 0

标准偏差为 2.3684e-016

54320.2449x0.2961