mcgm.m

基于泰勒级数的灰度预测

function [Xmcgm,y,P1,ERROR]=mcgm(x)
nc=length(x);
k=0;
[y,e]=gm1(x,k) ;
%-------mcgm(1,1)----程序
  %%%----区间划分的边界点求取----- 分为4个区间
r=4;n=nc+k;
stats=zeros(1,r+1);
r1=0:r-1;
stats(end)=max(e);%避免舍入误差
stats(1:end-1)=min(e)+(max(e)-min(e))*r1/r;
L=stats(1:end-1);
U=stats(2:end);
for i=1:r+1
   Y(i,:)=y(1:n)-stats(i);
end

  %%% 寻找原始点所在的区域-----
ys=Y;
rate=zeros(r,nc);
for j=1:nc
   AAA={[ys(1,j),ys(2,j)] [ys(2,j),ys(3,j)] [ys(3,j),ys(4,j)] [ys(4,j),ys(5,j)]};%如果r要增加,这个地方需要修改
   for i=1:r      
      c=find((AAA{i}(1))>=x(j)&AAA{i}(2)<x(j));
      if c==1,rate(i,j)=1;end
   end
  %-----修正最下方数据点在区间端点时的情况--------
  c0=find((AAA{4}(2))>=x(j)&(AAA{4}(2))<=x(j));
  if c0==1,rate(i,j)=1; end                     
end
%----判断是否所有数据点都在区间内部-----------------
xc=sum(rate);
s=sum(sum(rate,2));
if s==nc
  sprintf('恭喜您! 所有原始点在划分区间内,总数据为%d个',n)
elseif s>n
  sprintf('数据点有重复计算,请检查')
else
 %s1=n-s;num=find(xc>=0&xc<=0);
    sprintf('注意! 有%d个原始点不在区间内,请进行处理,程序已停止',s1)
    %num
end
%---------转移概率的计算----- 转移概率为P1------
for i=1:nc
   c(i)=find(rate(:,i)>0);
end   
cc=5-c;
P0=zeros(4,4);
for i=1:18
    c1=cc(i);c2=cc(i+1);
    P0(c1,c2)= P0(c1,c2)+1;
end
sumP0=sum(P0,2);
for j=1:r
   if sumP0(j)==0
       P(j,:)=0;
   else
    P(j,:)=P0(j,:)/sumP0(j);
   end
end

ss=fliplr(-stats);
L1=ss(1:end-1);U1=ss(2:end);
vc=(U1+L1)/2;
P1=rot90(P,90);
YF=flipud(Y);
Xmcgm(1)=y(1);
for j=2:nc
 Xmcgm(j)=y(j)+vc*P1(cc(j),:)';
end
hold on
tt=1:n;
plot(tt,Y,1:nc,Xmcgm,'-.b') % 画出区间分布图, 需要加如注释语句
hold off
%%%%%-----误差分析---------------------
AEgm=abs(y(1:nc)-x');
AEmcgm=abs((Xmcgm(1:nc)-x'));
AMEgm=sum(AEgm)/nc;               %gm(1,1)的误差和
MSEgm=(y(1:nc)-x')*(y(1:nc)-x')'/n;%
MSEmcgm=(Xmcgm(1:nc)-x')*(Xmcgm(1:nc)-x')'/nc;
AMEmcgm=sum(AEmcgm)/nc;

ERROR=[MSEgm      AMEgm;MSEmcgm    AMEmcgm];