MATALB运用——最小二乘法拟合

实验一

编制以函数 { x k } k = 0 n \{x^k\}_{k=0}^n {xk}k=0n​为基的多项式最小二乘拟合程序，并用于对下表中的数据作3次多项式最小二乘拟合。

取权函数$w_i=1$，求拟合曲线${\phi }^n={\sum }_{k=0}^n{\alpha }_k^n{x}^k$中的参数 { α k } \{\alpha_k\} {αk​}、平方误差${\delta }^2$，并作离散数据 { x i , y i } \{x_i,y_i\} {xi​,yi​}的拟合函数$y={\phi }^{\ast }(x)$的图形。

1、思路

在MATLAB中，多项式最小二乘法的拟合可以采用库函数polyfit()
p = polyfit(x,y,n) :返回次数为 n 的多项式 p(x) 的系数，该阶数是 y 中数据的最佳拟合（在最小二乘方式中）。p 中的系数按降幂排列，p 的长度为 n+1。

2、程序

function test_1
format long
%利用库函数polyfit做3次多项式最小二乘拟合
%求参数{a_k}、平方误差、拟合函数的图形
x=[-1.0,-0.5,0.0,0.5,1.0,1.5,2.0];
y=[-4.447,-0.452,0.551,0.048,-0.447,0.549,4.552];
p=polyfit(x,y,3);%3次多项式最小二乘拟合
y1=polyval(p,x);
r=sum((y-y1).^2);%平方误差
clf;
plot(x,y,'*');
hold on;
x0=[-1:0.01:2];
plot(x0,polyval(p,x0),'r-');
axis([-1.5,2.5,-5,5]);
xlabel('x轴');ylabel('y轴');title('3次多项式最小二乘拟合');
disp(['平方误差：',sprintf('%g',r)]);
disp(['参数：',sprintf('%g\t',p)]);

3、运行结果

平方误差：2.17619e-05
参数：1.99911 -2.99767 -3.96825e-05

实验二

编制正交化多项式最小二乘法拟合程序，并用于求解上题中的3次多项式最小二乘法拟合问题，做拟合曲线的图形，计算平方误差，并于实验一进行比较。

1、思路

2、程序

function test_3_2
format long
%利用正交化求数据的二次多项式拟合
%输出：拟合系数参数、平方误差、拟合图像
x=[-1.0,-0.5,0.0,0.5,1.0,1.5,2.0];
y=[-4.447,-0.452,0.551,0.048,-0.447,0.549,4.552];
w=[1,1,1,1,1,1,1];
n=3;
[a,b,c,alph,r]=flp(x,y,w,n)
disp(['平方误差：',sprintf('%g',r)]);
disp(['参数：',sprintf('%g\t',alph)]);
clf;
plot(x,y,'*');
hold on;
x0=[-1:0.01:2];
plot(x0,polyval(alph,x0),'r-');
axis([-1.5,2.5,-5,5]);
xlabel('x轴');ylabel('y轴');title('正交化的3次多项式最小二乘拟合');


function [a,b,c,alph,r]=flp(x,y,w,n)
%求正交递推公式的参数a、b，拟合多项式的系数c，整理之后的降幂系数alpha
m=length(x);
s1=ones(1,m);
v1=sum(w);
d(1)=y*w';c(1)=d(1)/v1;
for k=1:n
    xs=x.*s1.^2*w';
    a(k)=xs/v1;
    if (k==1)
        s2=(x-a(k)).*s1;
    else
        b(k)=v2/v11;
        s2=(x-a(k)).*s1-b(k)*s0;
    end
    v2=s2.^2*w';
    d(k+1)=y.*s2*w';c(k+1)=d(k+1)/v2;
    v11=v1;v1=v2;s0=s1;s1=s2;
end

%求平方误差r
r=y.*y*w'-c*d';

%拟合多项式合并同类项后的降幂系数
syms x
p0=1;
T=c(1)*p0;
for k=2:n+1
    if (k==2)
        p2=x-a(k-1);
    else
        p2=(x-a(k-1))*p1-b(k-1)*p0;
        p0=p1;
    end
    T=T+c(k)*p2;
    p1=p2;
end
T=collect(T);%合并同类项
alph=sym2poly(T);

3、运行结果

可见，但多项式次数为3时，两种方法拟合的结果和误差基本一致。