plesm.m

这是一个关于最小二乘法的程序

function [ g , E ] = Plesm( x,y,w,n )

%   正交多项式的最小二乘法曲线拟合
%   MATLAB FUNCTION Plesm.m
%
m = length(x);

sym y(m);
sym w(m);
sym n;

sym d(n+1);
sym PX(m,n+1);
sym P(n+1,n+1);

sym gp(n+1);
sym g(n+1);

for i = 1:n+1
    for j = 1:i
        if j ~= i
            P(i,j) = 0;
        else
            P(i,j) = 1;
        end
        
    end
end

sym s;

s = 0;
d(1) = 0;
for i = 1:m
    d(1) = d(1) + w(i);
    s = s + w(i).* x(i);
end

sym AF(n);
sym BA(n);
sym dd;

AF(1) = s./d(1);
P(2,1) = -AF(1);
for i = 1:m
    PX(i,1) = 1;
    PX(1,2) = x(i) - AF(1);
end

s = 0;
dd = 0;
BA(1) = 0;
for k = 2:n
    for i = 1:m
        s = s + w(i).*x(i).* PX(i,k).* PX(i,k);
        dd = dd + w(i).* PX(i,k).* PX(i,k);
    end
    d(k) = dd;
    AF(k) = s./d(k);
    BA(k) = d(k)./ d(k-1);
    P(k+1,1) = -AF(k - 1).* P(k,1) - BA(k - 1).* P(k-1,1);
    for i = 1:k-1
        j = k - i + 1;
        if j >= k
            t = 0;
        else
            t = P(k-j,j);    
        end
        P(k+1,j) = P(k,j-1) - AF(k-1).*P(k,j) - BA(k-1).* t;
    end
    
    for i = 1:m
        PX(i,k+1) = PX(i,k).* ( x(i) - AF(k-1) ) - BA(k-1).*PX(i,k-1);
    end
end

for k = 1:n+1
    temp = 0;
    for i = 1:m
        temp = temp + w(i).* PX(i,k).* PX(i,k);
    end
    d(k) = temp;
end

for i = 1:n+1
    temp = 0.0;
    for j = 1:m
        temp = temp + w(j).* y(j).* PX(j,i);
    end
    gp(i) = temp./ d(i);
end

for i = 1:n+1
    temp = 0.0;
    for j = i:n+1
        temp = temp + gp(j).* P(j,i);
    end
    g(i) = temp;
end

U = 0.0;
for i = 1:m
    U = U + w(i).* y(i).* y(i);
end

V = 0.0;
for i = 1:n+1
    V = V + gp(i).* gp(i).* d(i);
end

E = U - V;

%   OUT
g
E
%%%%%    为保证拟合的误差较小，取 n  的值稍微小点，如： n = 3
%  Drow out the picture
for i = 1:m
    plot(x(i),y(i),'b-*');
    xlabel('X Label');
    ylabel('Y Label');
    title('最小二乘法的曲线拟合');
    legend('蓝为已知点连线','红为拟合曲线',-1);
    grid on;
    hold on;
end
plot(x,y,'b--');

sym x0;
sym y0;
for i = 1:2.* m
    x0(i) = i.* 2./m;
    y0(i) = g(1);
    for k = 2:n+1
        y0(i) = y0(i) + g(k).* x0(i)^(k-1);
    end
end

plot(x0,y0,'r-*');


  
    

%  Testing

%  >> x = [0.24,0.65,0.95,1.24,1.73,2.01,2.23,2.52,2.77,2.99];
%  >> y = [0.23,-0.26,-1.10,-0.45,0.27,0.10,-0.29,0.24,0.56,1.00];
%  >> w = [1,1,0.8,0.9,1,1,1,1,0.9,0.9];
%  >> n = 5;
%  >> Plesm(x,y,w,n)
%  g =

%    2.4210   -9.1257   -2.0768    4.8929   -4.4954    1.0843


%  E =

%    2.3955