predictbss01.m

该程序对一小区的电话流进行预测

         if( D_index-d_r<1 )
            '天预测准备数据不全，情况1'
             return
         end
         %若前d_r天该时刻的数据不为0，则取其数据
         if( origanatedata (D_index-d_r,point+1)~=0  )
             d_A(d_id)=origanatedata (D_index-d_r,point+1);
             d_id=d_id+1;
         end
         d_r=d_r+1; 
     end
     avg_d=sum(d_A(:,1))/d;  % 求得天取值的平均值
    j2=d+1;
    %下面判断所取的天的数据是否出现偶然大幅度波动的情况
    for i=1:d    %  B  for 4
     % 下面判断，若出现偶然涨幅很大的值，则重新取值
       if  or(d_A(i)>avg_d * 3,d_A(i)< avg_d *0.25)  % B  if 2
           for  k=i:d-1  %  B  for  5
               d_A(k)=d_A(k+1);
           end          %  E for  5
          if( D_index-j2<1 )
            '天预测准备数据不全，情况22'
             return
           end
           d_A(d)=origanatedata(D_index -j2,point+1);
           avg_d=sum(d_A(:,1))/d;  %再次求得平均值
           j2=j2+1;
           i=1;
       end      %  E  if 2
   end         %  E  for 4    
   % 计算各个天对应点的涨幅值
      for i=1:d
         dd_r(i)=( d_A(i,1) - f(point+1) )/f(point+1) ;  
      end          
   %计算预测天情况影响的涨幅值
     dh0= dd_r(1)*0.3551 + dd_r(2)*0.2732 + dd_r(3)*0.2101 + dd_r(4)*0.1616 ;   %得到本天该时刻的涨幅
  %2-------处理完近几天的数据涨幅----------
   
  
  %------下面开始处理近几周的数据涨幅--------
      %  取前3周该时刻的数据，放于数组w_A【1。。3】中
       w=3;
       w_A=zeros( w);
       ww_r=zeros(w);
       w_r=1;%天时间变量
       w_id=1; %天下表变量
      
       %下面取取前1~w=3周的数据
       while w_id<=w 
         if( D_index-w_r*7<1 )
            '周预测准备数据不全，情况1'
             return
         end
          %若前w_r周该时刻的数据不为0，则取其数据
         if( origanatedata (D_index-w_r*7,point+1)~=0  )
             w_A(w_id)=origanatedata (D_index-w_r*7,point+1);
             w_id=w_id+1;
          end
             w_r=w_r+1;
     end
     avg_w=sum(w_A(:,1))/w;  % 求得周取值的平均值
    j3=w+1;
    %下面判断所取的周的数据是否出现偶然大幅度波动的情况
    for i=1:w    %  B  for 4
     % 下面判断，若出现偶然涨幅很大的值，则重新取值
       if  or(w_A(i)>avg_w * 3,w_A(i)< avg_w *0.25)  % B  if 2
           for  k=i:w-1  %  B  for  5
               w_A(k)=w_A(k+1);
           end          %  E for  5
          if( D_index-j3*7<1 )
            '周预测准备数据不全，情况2'
             return
           end
           w_A(w)=origanatedata(D_index -j3*7,point+1);
           avg_w=sum(w_A(:,1))/w;  %再次求得平均值
           j3=j3+1;
           i=1;
       end      %  E  if 2
   end         %  E  for 4    
   % 计算各个周对应点的涨幅值
      for i=1:w
         ww_r(i)=( w_A(i,1) - f(point+1) )/f(point+1) ;  
      end          
   %计算预测周情况影响的涨幅值
     wh0= ww_r(1)*0.3551 + ww_r(2)*0.2732 + ww_r(3)*0.2101 + ww_r(4)*0.1616 ;   %得到本天该时刻的涨幅
  %2-------处理完近几周的数据涨幅----------
    
  %－－－－－－下面开始作预测计算－－－－－－－
        if (ballancefactor_h<0.01 )
            %得到该时刻的预测值
            lenarray(1,1)=lenarray(1,1)+1;
            R(point+1,1)=f(point+1)*(1 + hh0 ); 
        elseif ( ballancefactor_h>=0.10 )
            lenarray(2,1)=lenarray(2,1)+1;
            w0=0.40;     w1=0.60; 
            %得到该时刻的预测值
            R(point+1,1)=f(point+1)*(1 + dh0*w0 + hh0*w1 );
       else
           lenarray(3,1)=lenarray(3,1)+1;
           w0=0.25;      w1=0.35;       w2=0.40;
           %得到该时刻的预测值
           R(point+1,1)=f(point+1)*(1 + wh0*w0 + dh0*w1 + hh0*w2 ); 
        % 下面几行代码留着扩展备用    
        %else    
        %  lenarray(4,1)=lenarray(4,1)+1;
         %  w0=0.55;     w1=0.45;  
          %%得到该时刻的预测值
           %R(point+1,1)=f(point+1)*(1 + wh0*w0 + dh0*w1 );       
      end
end     %   E  for  1
  %－－－－－－预测计算处理完毕－－－－－－－
 end

 
 %得到该天的预测值
predict_result=R
lenarray;

%－－－－预测误差计算统计处理开始－－－－－
   avg_predict=0;
   max_predict=zeros(1,6);
   max_predict(1,6)=0;
   min_predict=zeros(1,6);
   min_predict(1,6)=100;
%下面以该天的数据对模型做检验，画出走势图型
 %该天的实际值数据
 innormal_max=0;
 innormal_0=0;
 innormal=zeros(1,6);
 zerospoint=zeros(1,5);
 realdata=zeros(24,1);
   for i1=1:24
       realdata(i1)=origanatedata(base,i1);
       if ( realdata(i1)~=0 )
          pp0=abs ((  realdata(i1)-f(i1) ) / f(i1) *100);
          pp=abs (( predict_result(i1) - realdata(i1) ) / realdata(i1) *100);
          if ( pp0<=50 )
            avg_predict = avg_predict + pp;
            if pp>max_predict(1,6);
                max_predict(1,1)=base;
                max_predict(1,2)=i1;
                max_predict(1,3)=f(i1);
                max_predict(1,4)=realdata(i1);
                max_predict(1,5)=predict_result(i1);
                max_predict(1,6)=pp;
            end
            if pp < min_predict(1,6)
               min_predict(1,1)=base;
               min_predict(1,2)=i1;
               min_predict(1,3)=f(i1);
               min_predict(1,4)=realdata(i1);
               min_predict(1,5)=predict_result(i1);
               min_predict(1,6)=pp;               
            end
          else
            innormal(1,1)=base;
            innormal(1,2)=i1;
            innormal(1,3)=f(i1);
            innormal(1,4)=realdata(i1);
            innormal(1,5)=predict_result(i1);
            innormal(1,6)=pp;
            tubiandian=innormal;
            innormal_max=innormal_max+1;
         end
      else 
            zerospoint(1,1)=base;
            zerospoint(1,2)=i1;
            zerospoint(1,3)=f(i1);
            zerospoint(1,4)=0;
            zerospoint(1,5)=predict_result(i1);
            zeropoint_0=zerospoint;
            innormal_0=innormal_0+1;
      end
   end
  
%'标志序号1 突变点个数：',innormal_max;
%'标志序号2 零点个数：',innormal_0;
%'标志序号3 正常情况下平均预测误差为(%)： ',avg_predict/24;  
%'标志序号4 正常情况下最大预测误差为(%)： ', max_predict  ;
%'标志序号5 正常情况下最小预测误差为(%)： ',min_predict;
%－－－－预测误差计算统计处理结束－－－－－


% 做实际走势图
plot(x,realdata,'r')
xlabel('0 -- 23  ')
title('预测2002年8月19日，BSSO1_171通话走势图 （模型检验）')
  
 % 做预测走势图
hold on
plot(x,predict_result,'g')

 % 做平均走势图
hold on
plot(x,f,'b')
legend('红色为实际走势图','绿色为预测走势图','蓝色为平均走势图',3);     
hold off