/* Repetitive Moving Median */
new; cls;
let data[50,1]=
      60
      61
      60
      62
      60
      63
      64
      65
      69
      70
      69
      66
      64
      64
      62
      64
      65
      55
      53
      57
      60
      58
      63
      63
      63
      68
      69
      69
      68
      69
      70
      67
      66
      67
      68
      69
      69
      77
      80
      80
      78
      76
      76
      72
      71
      70
      73
      71
      78
      80
;
data=rev(data);               /* original data in reverse order */
t=seqa(1,1,rows(data));

library pgraph;
graphset;
_pltype=6;
_plegctl=1;
_plegstr="actual\000Single MMed(k)\000Repetitive MMed(k)";
title("Repetitive Moving Median");
k=7;
xy(t,data~cmmed(data,k)~repmmed(data,k));




proc repmmed(x,k);
   local m,i,j;
   if k%2==0;
      errorlog "ERROR: Parameter k must be odd number.";
      retp(".");
   endif;
   m=zeros(rows(x),cols(x)); 
   j=1;
   do until x==m;
      if j>=2;
         x=m;
      endif;
      i=k;
      do while i<=rows(x);
         m[i-(k-1)/2,.]=median(x[i-k+1:i,.])';
         i=i+1;
      endo;
      m[1:(k-1)/2,.]=x[1:(k-1)/2,.];
      m[rows(x)-(k-1)/2+1:rows(x),.]=x[rows(x)-(k-1)/2+1:rows(x),.];
      print/rz "iteration" j;
      j=j+1;
   endo;
   m[1:(k-1)/2,.]=miss(zeros((k-1)/2,1),0);
   m[rows(x)-(k-1)/2+1:rows(x),.]=miss(zeros((k-1)/2,1),0);
   retp(m);
endp;




/* the procedure below is not directly used by repmmed */
proc cmmed(x,k);
   local m,i;
   if k%2==0;
      errorlog "ERROR: Parameter k must be odd number.";
      retp(".");
   endif;
   m=zeros(rows(x),cols(x));
   i=(k-1)/2+1;
   do while i<=rows(x)-(k-1)/2;
      m[i,.]=median(x[i-(k-1)/2:i+(k-1)/2,.])';
      i=i+1;
   endo;
   m[1:(k-1)/2,.]=miss(m[1:(k-1)/2,.],0);
   m[rows(x)-(k-1)/2+1:rows(x),.]=miss(m[rows(x)-(k-1)/2+1:rows(x),.],0);
   retp(m);
endp;