fit.m

不含Matlab 内部函数

function Fit(X,Y,maxM)
%进行一倒maxM阶多项式拟合fit(X,Y) 中参数代表x,y矩阵 maxM是最大阶数
%%输入数据检验begin
if nargin~=3
    error('函数缺少参数，参数分别为x,y矩阵和最大阶数，参数数目为三！')
end
[mx,nx]=size(X);
[my,ny]=size(Y);
if nx~=1|ny~=1|my~=mx|my<3
    error('A或B参数输入数据不合法！，程序退出！注意：X,Y需是列向量！');
    return;
end
if maxM<2
    error('最大阶数太小！')
end
%%输入数据检验end
N=mx;%数据数目
%M=1-maxM变化拟合阶数
for M=1:maxM
  a=FitProc(X,Y,M,N);%%%得到拟合系数
%%%%%%%%%%%%%%求残差平方和---begin%%%%%%%%%%%%%
QQ=0;%残差平方和
for i=1:N
    QQ=QQ+(Y(i)-FitFun(M,a,X(i)))^2;
end
qq=QQ/(N-M-1);
if(M==1)
     qq_return=qq;
     continue;
else
    if abs(qq)<abs(qq_return)
        qq_return=qq;
        M_return=M;
        QQ_return=QQ;
    end
end
%%%%%%%%%%%%%%求残差平方和---end%%%%%%%%%%%%%
end
%%%%%%%%%%%结果输出-----begin%%%%%%%%%%%%%
clc;%清屏
disp('%%%%%%%%%%%% 多项式结果输出 %%%%%%%%%%%%%%')
a=FitProc(X,Y,M_return,N);
disp('拟合方差：');
disp(qq_return);
disp('拟合阶数：');
disp(M_return);
disp('拟合系数：');
disp(a);
%%%%%%%%%输出方程式---begin
str=sprintf('y=%f',a(1));
for i=2:M_return+1
     str=strcat(str,sprintf('+(%f)x^%d',a(i),(i-1)));
end
disp('拟合方程：');
disp(str)
%%%%%%%%%输出方程式---end
%%%%%%%作图%%%%%%%%%%%
t=floor(min(X)):0.1:ceil((max(X)));
[temp1,temp2]=size(t);
z=zeros(1,temp2);
for i=1:temp2
    z(i)=FitFun( M_return,a,t(i));
end
plot(X',Y','*',t,z);
str=sprintf('%d阶多项式拟合图像',M_return);
title(str);
xlabel('X');
Ylabel('Y');
%%%%%%%作图%%%%%%%%%%%
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
%%%%%%%%%%%结果输出-------end%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%主函数结束%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%定义子函数%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%  FitProc  %%%%%%%%%%%%%%%%%%%%%%
function a=FitProc(X,Y,M,N)
%多项式拟和子函数
C=zeros(M+1,M+1);%系数矩阵
B=zeros(M+1,1);%常数矢量
A=zeros(M+1,1);%a0 a1 ....am
%%%%%取得系数矩阵C--begin
SUMXX=zeros(2*M);
%1次到2m次x的平方和SUMXX(m)=x^m的平方和
for i=1:2*M
    SUMXX(i)=SumXX(X,i);
end
C(1,1)=N;
for i=2:M+1
    C(1,i)=SUMXX(i-1);
end
for i=2:M+1
    k=i-1;
    for j=1:M+1
        C(i,j)=SUMXX(k);
        k=k+1;
    end
end
%%%%%取得系数矩阵C--end
%%%%%取得常数矢量矩阵B--begin
for i=1:M+1
    sum=0;
    for j=1:N
        sum=sum+Y(j)*X(j)^(i-1);
    end
    B(i)=sum;
end
%%%%%取得常数矢量矩阵B--end
%%%%%%a0 a1 ....am矩阵A--begin
A=inv(C)*B;
%%%%%%a0 a1 ....am矩阵A--begin
a=A;
return;
%%%%%%%%%%%%%%%%%  SumXX  %%%%%%%%%%%%%%%%%%%%%%
function Sum_XX=SumXX(X,m)
%求X矩阵中x的m次方的平方和
[mm,nn]=size(X);
sum=0;
for i=1:mm
    sum=sum+power(X(i),m);
end
Sum_XX=sum;
return;
%%%%%%%%%%%%%%%%%  FitFun  %%%%%%%%%%%%%%%%%%%%%%%
%%%%%根据拟合常数a0,a1.....求拟合的值
function y=FitFun(M,a,x)
sum=0;
for i=1:M+1
    sum=sum+a(i)*x^(i-1);
end
y=sum;
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%定义子函数%%%%%%%%%%%%%%%%%%%%%%%%%