%x,y應是原始數據，下面的數據僅是為了說明問題而自設的示例數據。
%f的形式應由具體的模型推導得到，要符合實際情況，不可牽強，否則擬合效果不好！
%x0的選擇也很重要，應和理論值比較接近，否則也可能擬合效果不好！
%p是根據數據和擬合的函數形式而得到的參數值。
x=0:0.1:1;
y=10./(1+exp(-0.01.*x))+1;
f=inline('p(1)./(1+exp(-p(2)*t))','p','t');
x0=[5,0.05];
p=nlinfit(x,y,f,x0);
plot(x,y,'ro') 
hold on
xx=0:0.01:1.2;
yy=f(p,xx);

plot(xx,yy) 






已知觀察數據如下表所示，按下屬方案求最小二乘擬合函數，并求出偏差平方和，比較擬合曲線的優劣。
  x:0 0.2 0.6 1.0 1.3 1.6 1.7 1.8 1.9 2.2 2.3 2.5 2.6
  y:0 -2.5 -4.0 -5.7 -3.5 -2.0 -1.0 2.0 3.5 4.0 7.0 7.5 9.9
  x:2.9 3.1 3.4 3.8 4.1 4.4 4.7 4.8 4.9 5.0 5.1 5.3
  y:10.9 11.9 13.5 13.0 11.9 9.0 6.5 4.0 1.5 0.0 -2.5 -5.0

 

 

%用離散正交多項式求三次擬合多項式
% x,y--表示原始數據的節點坐標
% w--表示權重系數
% N--表示要擬合的離散正交多項式的最高次數
% polyapproximate()--是自定義函數，可以求解多項式的系數
% 其返回值c為多項式系數，error為偏差平方和
x=[0 0.2 0.6 1.0 1.3 1.6 1.7 1.8 1.9 2.2 2.3 2.5 2.6 2.9 3.1 3.4 3.8 4.1 4.4 4.7 4.8 4.9 5.0 5.1 5.3];
nn=length(x);
for i=1:nn
    w(i)=1;
end
y=[0 -2.5 -4.0 -5.7 -3.5 -2.0 -1.0 2.0 3.5 4.0 7.0 7.5 9.9 10.9 11.9 13.5 13.0 11.9 9.0 6.5 4.0 1.5 0.0 -2.5 -5.0];
N=3;%此處可取3 or 4.
[c,error]=polyapproximate(x,y,w,N)
t=0:0.1:5.3;
u=polyval(c,t);
plot(t,u,x,y,'+')


 

%自定義函數polyapproximate（），用來做離散正交多項式擬合
% 此函數的作用是做不同次數的離散正交多項式的擬合
% X，Y 為原始數據的坐標值矩陣
% w 為權重系數
% N 為離散正交多項式的最高次數
function [C,E]=polyapproximate(X,Y,w,N)
M=length(X);
for i=1:N+1
    for j=1:i
        if j~=i
            P(i,j)=0;
        else
            P(i,j)=1;
        end
    end
end
S=0;
d(1)=0;
for i=1:M
    d(1)=d(1)+w(i);
    S=S+w(i)*X(i);
end
AF(1)=S/d(1);
P(2,1)=-AF(1);
for i=1:M
    PX(i,1)=1;
    PX(i,2)=X(i)-AF(1);
end
BA(1)=0;
for k=2:N+1
    S=0;
    dd=0;
    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-1)=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-1,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
d(N+1)=0;
for i=1:M
    d(N+1)=d(N+1)+w(i)*PX(i,N+1)*PX(i,N+1);
end
for i=1:N+1
    FM=0;
    for k=1:M
        FM=FM+w(k)*Y(k)*PX(k,i);
    end
    gp(i)=FM/d(i);
end
for i=1:N+1
    C(i)=0;
    for j=i:N+1
        C(i)=C(i)+gp(j)*P(j,i);
    end
end
C=flipud(C');
%C=C'
U=0;
for i=1:M
    U=U+w(i)*Y(i)*Y(i);
end
V=0;
for k=1:N+1
    V=V+gp(k)*gp(k)*d(k);
end
E=U-V;