最早接触到最小二乘算法是在大四上学期的系统建模与辨识的课上，此课也是大四为数不多的认真做作业的课程之一。最小二乘算法原理简单，实用性好，编程很容易实现。

    最小二乘算法有递推最小二乘、辅助变量最小二乘、增广矩阵法等，多阶段最小二乘算法等，大四上学期和研一上学期的系统辨识中对递推、增广的算法曾编程过，这里我们讨论一种很简单的多项式最小二乘的实现。其实这个题目也是袁牛要做的数值分析中的一个。

    题目：给定数表：x=-0.75 -0.5 -0.25 0 0.25 0.5 0.75，y=0.33 0.88 1.44 2.00 2.56 3.13 3.71，试分别用一次、二次、三次多项式根据最小二乘原则拟合这些数据，并比较优劣。

    原理：已知X1 X2...Xm和Y1 Y2...Ym，求一个简单易算的近似函数y(x)去代替数据反映的关系f(x)，对于近似函数，当然可用插值方法进行处理，但是当数据量太大的时候，差值方法有着明显的不足。（1）大量的实验数据难以保证每个数据的精确性，当其中数据存在一定误差的时候，由于插值条件的要求，其误差将完全被插值函数进一步继承。（2）即使所有数据都精确，但是多项式次数过高或者过多的分段处理，都会带来相应的问题。其实，我们没有必要要求近似函数和原数据在采样点值相同，这时需要引入拟合的思想，使得在采样点的误差组合最小即可。

    确定多项式y(x)=C0+C1*X+C2*X^2+...Cn*X^n，使得在采样点处(y(Xi)-Yi)^2的总和最小。对于具体的公式，随便找本书都有的，在这里需要说明的是，所需要的点数即等式的个数要多于需要估计的参数的个数。我一向觉得最重要的是把公式Estimate=inv(H'H)H'y中的H找出来，待估计量Estimate是一个列向量，Yi=[1 Xi Xi^2...Xi^n]*[C0 C1 C2...Cn]'，等式右半部份的前面即是H，从而我们得到如下的程序：

clc;
close all;

%原始数据
Y=[0.33 0.88 1.44 2.00 2.56 3.13 3.71];%给定Y数值
X=[-0.75 -0.5 -0.25 0 0.25 0.5 0.75];  %给定X数值
L=length(Y);                           %给定的数据个数

for i=1:1:L
     row1(i)=X(i);
     row2(i)=X(i)^2;
     row3(i)=X(i)^3;
end

X1=[ones(L,1),row1'];             %构造矩阵L×2
X2=[ones(L,1),row1',row2'];       %构造矩阵L×3
X3=[ones(L,1),row1',row2',row3']; %构造矩阵L×4

%一阶估计
Estimate1=inv(X1'*X1)*X1'*Y';   
for i=1:1:L
    EZ1(i)=Estimate1(2,1)*X(i)+Estimate1(1,1); 
end

%二阶估计
Estimate2=inv(X2'*X2)*X2'*Y';
for i=1:1:L
    EZ2(i)=Estimate2(3,1)*X(i)^2+Estimate2(2,1)*X(i)+Estimate2(1,1);
end

%三阶估计
Estimate3=inv(X3'*X3)*X3'*Y';
for i=1:1:L
    EZ3(i)=Estimate3(4,1)*X(i)^3+Estimate3(3,1)*X(i)^2+Estimate3(2,1)*X(i)+Estimate3(1,1);
end

figure(1);
plot(X,Y,'+',X,EZ1,'-',X,EZ2,'.',X,EZ3,'^');
xlabel('时间');
ylabel('数据曲线');
title('最小二乘法估计曲线');
legend('原系统','一阶估计','二阶估计','三阶估计','Location','NorthWest');
figure(2);
plot(X,Y,'-',X,EZ1,'^');
xlabel('时间');
ylabel('数据曲线');
title('一阶最小二乘法估计曲线');
legend('原系统','一阶估计','Location','NorthWest');
figure(3);
plot(X,Y,'-',X,EZ2,'+');
xlabel('时间');
ylabel('数据曲线');
title('二阶最小二乘法估计曲线');
legend('原系统','二阶估计','Location','NorthWest');
figure(4);
plot(X,Y,'-',X,EZ3,'+');
xlabel('时间');
ylabel('数据曲线');
title('三阶最小二乘法估计曲线');
legend('原系统','三阶估计','Location','NorthWest');

%平方误差信息
error1=0;
error2=0;
error3=0;
for i=1:1:L
    error1=error1+(Y(i)-EZ1(i))^2;
    error2=error2+(Y(i)-EZ2(i))^2;
    error3=error3+(Y(i)-EZ3(i))^2;
end
e1=sqrt(error1)
e2=sqrt(error2)
e3=sqrt(error3)

%输出系数信息
%digits(n) 控制总位数 vpa(s,n)控制有效位数
Estimate1
Estimate2
Estimate3

    在主窗口得到的输出信息如下：
e1 =

   0.022677868380554

e2 =

   0.006172133998484

e3 =

   0.004629100498863

Estimate1 =

   2.007142857142857
   2.251428571428571

Estimate2 =

   1.997619047619047
   2.251428571428571
   0.038095238095238

Estimate3 =

   1.997619047619047
   2.243650793650794
   0.038095238095238
   0.017777777777781

    分别为平方误差和、系数，可见随着阶数的增加，拟合误差减小。

http://s2/bmiddle/4b013fb1n8b00a5971391&690

图1 数据拟合曲线图

    当然了，matlab中也有自带的最小二乘拟合函数polyfit(X,Y,n)，拟合所给数据，三次拟合的时候，主窗口如下输出：0.017777777777778 0.038095238095239 2.243650793650793 1.997619047619048 系数是从高到低输出，1.997619047619047 2.243650793650794  0.038095238095238 0.017777777777781 系数是从低到高输出，两者几乎完全相同。

email:dingqian12345@126.com

2010年07月10日上午于njust automation 101 房间

参考文献

[1] 李言俊，张科.系统辨识理论及应用[M].国防工业出版社，2006.

[2] 李建良，蒋勇，汪光先.计算机数值方法[M].东南大学出版社，2000.

CopyRight：版权所有若需转载或使用请联系作者

Email:dingqian12345@126.com