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