c语言中线性与非线性,最小二乘法 线性与非线性拟合

最小二乘法实现数据拟合

最小二乘法原理

函数插值是差值函数p(x)与被插函数f(x)在节点处函数值相同，即p( )=f( ) (i=0,1,2,3……,n)，而曲线拟合函数 不要求严格地通过所有数据点( )，也就是说拟合函数 在 处的偏差

=

不都严格地等于零。但是，为了使近似曲线能尽量反应所给数据点的变化趋势，要求| |按某种度量标准最小。

即

=

为最小。这种要求误差平方和最小的拟合称为曲线拟合的最小二乘法。

(一)线性最小二乘拟合

根据线性最小二乘拟合理论，我们得知关于系数矩阵A的解法为A＝R\Y。

例题 假设测出了一组，由下面的表格给出，且已知函数原型为

y(x)=c1+c2*e^(-3*x)+c3*cos(-2*x)*exp(-4*x)+c4*x^2

x

0

0.2

0.4

0.7

0.9

0.92

0.99

1.2

1.4

1.48

1.5

y

2.88

2.2576

1.9683

1.9258

2.0862

2.109

2.1979

2.5409

2.9627

3.155

3.2052

试用已知数据求出待定系数的值。

在Matlab中输入以下程序

x=[0,0.2,0.4,0.7,0.9,0.92,0.99,1.2,1.4,1.48,1.5]';

y=[2.88;2.2576;1.9683;1.9258;2.0862;2.109;

2.1979;2.5409;2.9627;3.155;3.2052];

A=[ones(size(x)) exp(-3*x),cos(-2*x).*exp(-4*x)

x.^2];

c=A\y;

c'

运行结果为

ans =

1.2200 2.3397 -0.6797 0.8700

下面画出由拟合得到的曲线及已知的数据散点图

x1=[0:0.01:1.5]';

A1=[ones(size(x1))

exp(-3*x1),cos(-2*x1).*exp(-4*x1) x1.^2];

y1=A1*c;

plot(x1,y1,x,y,'o')

事实上，上面给出的数据就是由已知曲线

y(x)= 0.8700-0.6797*e^(-3*x)+

2.3397*cos(-2*x)*exp(-4*x)+ 1.2200*x^2

产生的，由上图可见拟合效果较好。

多项式最小二乘拟合

在Matlab的线性最小二乘拟合中，用得较多的是多项式拟合，其命令是

A=polyfit(x,y,m)

其中 表示函数中的自变量矩阵， 表示因变量矩阵，是输出的系数矩阵，即多项式的系数。

多项式在自变量x处的函数值y可用以下命令计算：

y=polyval(A,x)

例题 对下面一组数据作二次多项式拟合，即要求出二次多项式 中的 ，使最小。

x

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

y

－0.447

1.978

3.28

6.16

7.08

7.34

7.66

9.56

9.48

9.30

11.2

在Matlab中输入以下命令

x=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)

运行结果为

a =

-9.8108 20.1293 -0.0317

f=vpa(poly2sym(a),5)

%vpa(polyval2sym()，n)只适用于关于多项式函数的拟合。因为此函数对于自变量统一规定为“x”，将由polyfit()所得出的系数按自变量幂次升降放在相应的位置。

运行结果为

f =

-9.8108*x^2+20.129*x-.31671e-1

下面画出由拟合得到的曲线及已知的数据散点图

y1=polyval(a,x);

plot(x,y,'o',x,y1)

(二)非线性最小二乘拟合

(1)lsqcurvefit( )

lsqcurvefit( )是非线性最小二乘拟合函数，其本质上是求解

最优化问题。其使用格式为

x=lsqcurvefit(fun,x0,xdata,ydata)

其中，fun是要拟合的非线性函数，x0是初始参数，xdata，ydata是拟合点的数据，该函数最终返回系数矩阵。

例题 假设已知

并已知该函数满足原型为 ，其中 为待定系数。

在Matlab中输入以下命令

x=0:.1:10;

y=0.12*exp(-0.213*x)+0.54*exp(-0.17*x).*sin(1.23*x);

f=inline('a(1)*exp(-a(2)*x)+a(3)*exp(-a(4)*x).*sin(a(5)*x)','a','x');

%建立函数原型，则可以根据他来进行下面的求取系数的计算

[a,res]=lsqcurvefit(f,[1,1,1,1,1],x,y);

a',res

运行结果为

ans =

0.1197

0.2125

0.5404

0.1702

1.2300

res =

7.1637e-007

所求得的系数与原式中的系数相近。

如果向进一步提高精度，则需修改最优化的选项，函数的调用格式也将随之改变。

在Matlab中输入以下命令

ff=optimset;ff.TolFun=1e-20;ff.TolX=1e-15;%修改精度，即是修改其限制条件

[a,res]=lsqcurvefit(f,[1,1,1,1,1],x,y,[],[],ff);%两个空矩阵表示系数向量的上下限

a',res

运行结果为

ans =

0.1200

0.2130

0.5400

0.1700

1.2300

res =

9.5035e-021

下面绘图

x1=0:0.01:10;

y1=f(a,x1);

plot(x1,y1,x,y,'o')

(2)lsqnonlin( )

lsqnonlin( )函数是另一种求最小二乘拟合的函数，其本质上是求解最优化问题

最优化解。它的应用格式为

x=lsqnonlin(fun,x0)

其中，fun的定义与lsqnonlin( )函数中fun的定义有差别，

x0仍为初始参数向量，将输出的系数结果放在变量 中。

说明lsqnonlin()函数的使用方法。

首先编写目标函数 (curve_fun.m)

function y=curve_fun(p)%非线性最小二乘拟合函数

x=[0.02 0.02 0.06 0.06 0.11 0.11 0.22 0.22 0.56

0.56 1.10 1.10];

y=[76 47 97 107 123 139 159 152 191 201 207

200];

y=p(1)*x./(p(2)+x)-y;

再用lsqnonlin()函数求解，输入

[p,resnorm,residual]=lsqnonlin(@curve_fun,[200,0.1])

运行结果为

p =

212.6836 0.0641

resnorm =

1.1954e+003

residual =

Columns 1 through 11

-25.4339 3.5661 5.8111 -4.1889 11.3617 -4.6383 5.6847 12.6847 -0.1671 -10.1671 -6.0313

Column 12

0.9687

上面的两种方法都可以作非线性最小二乘曲线拟合

(3)非线性函数的线性化

在进行非线性拟合时，以往由于计算机和相关软件水平有限，常常先把非线性的曲线拟合线性化，然后再进行拟合。下面比较一下先线性化再进行拟合和直接进行非线性拟合的差异。

例题 已知数据

t

0.25

0.5

1

1.5

2

3

4

6

8

c

19.21

18.15

15.36

14.10

12.89

9.32

7.45

5.24

3.01

满足曲线 通过数据拟合求出参数和 。

方法一：先将非线性函数转化为线性函数

编写Matlab程序如下

d=300;

t=[0.25 0.5 1 1.5 2 3 4 6 8];

c=[19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24

3.01];

y=log(c);

a=polyfit(t,y,1)

运行结果为

a =

-0.2347 2.9943

k=-a(1)

k =

0.2347

v=d/exp(a(2))

v =

15.0219

由此也可以求出相关系数。

方法二：应用非线性拟合直接求解系数

建立m文件：

function f=curvefun3(x,tdata)

d=300

f=(x(1)\d)*exp(-x(2)*tdata)%x(1)=v;x(2)=k

运行程序

tdata=[0.25 0.5 1 1.5 2 3 4 6 8];

cdata=[19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24

3.01];

x0=[10 0.5];

x=lsqcurvefit('curvefun3',x0,tdata,cdata)

运行结果为

x =

14.8212 0.2420

下面绘图

f=curvefun3(x,tdata);

plot(tdata,cdata,'o',tdata,f)

我们发现两种求法求出的系数很接近。

(三)线性拟合和非线性拟合区别与联系

在许多实际问题中，变量之间内在的关系并不想前面说的那样简单。呈线性关系，但有些非线性拟合曲线可以通过适当的变量替换转化为线性曲线，从而用线性拟合进行处理。对于一个实际的曲线拟合问题，一般先根据观测值在直角坐标平面上描出散点图，看一看散点的分布同哪类曲线图形接近，让后选用相接近的曲线拟合方程，再通过适当的变量替换转化为线性拟合问题，按线性拟合解出后再还原为原变量所表示的曲线拟合方程。

表1.1

线性拟合方程

变量变换

变换后线性拟合方程

Y=

,

Y=a

Y=a

Y=

,

Y=

Y=

例题

测出一组实际数据 见下表 是对其进行函数拟合。

X

1.1052

1.2214

1.3499

1.4918

1.6478

3.6693

Y

0.6795

0.6006

0.5309

0.4693

0.4148

0.1546

X

1.8221

2.0138

2.2255

2.4596

2.7183

Y

0.3666

0.3241

0.2864

0.2532

0.2238

>>x=[1.1052,1.2214,1.3499,1.4918,1.6478,1.8221,2.0138,2.2255,2.4596,2.7183,3.6693];

>>y=[0.6795,0.6006,0.5309,0.4693,0.4148,0.3666,0.3241,0.2864,0.2532,0.2238,0.1546];

>>plot(x,y,x,y,'*')

见下图

由上图可以看出是非线性曲线y=a  由表1.1第一行可知

Y=

,

分别对x，y进行对数变换

>>x1=log(x);y1=log(y);plot(x1,y1)

用线性函数拟合的方法可以求出现行参数，使得ln y=aln x+b,即y= ,求解系数a,b及 。

>>

A=[x1',ones(size(x1'))];c=[A\y1']

c =

-1.2338

-0.2631

>> exp(c(2))

ans =

0.7686

拟合函数为y(x)=0.768 703 388 199 24