matlab nihe.txt

用matlab实现拟合的几个例子源码

%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;