ARL With Supplementary Run Rules

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: SHWARL2                                             */
 /*   TITLE: ARL With Supplementary Run Rules                    */
 /* PRODUCT: QC                                                  */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: Shewhart Charts, Average Run Lengths,               */
 /*   PROCS: IML GPLOT                                           */
 /*    DATA:                                                     */
 /*                                                              */
 /*   NOTES:                                                     */
 /*    MISC:                                                     */
 /*                                                              */
 /*     REF: SAS/QC Software:  Usage and Reference, Version 6,   */
 /*          First Edition, Volume 1 and Volume 2                */
 /*                                                              */
 /*          Champ C. W. and Woodall W. H. (1990), "A Program to */
 /*          Evaluate the Run Length Distribution of a Shewhart  */
 /*          Control Chart with Supplementary Run Rules"         */
 /*          Journal of Quality Technology, Vol. 22, pg. 68-73.  */
 /*                                                              */
 /*          Champ C. W. and Woodall W. H. (1987), "Exact        */
 /*          Results for Shewhart Control Chart With             */
 /*          Supplementary Runs Rules". Technometrics, Vol. 29,  */
 /*          pg. 393-399.                                        */
 /*                                                              */
 /*  Examples in the documentation were created using the        */
 /*    statement:  ods listing style=statistical;                */
 /*                                                              */
 /****************************************************************/
 options  ps=60 ls=80 nodate;


 title1 'ARL with Supplementary Runs Rules';

 * Input the Rules Defining the Chart ;

 data rules;
   input k m a b;
   cards;
   1  1   3   9
   2  3   2   3
   2  3  -3  -2
   1  1  -9  -3
   ;

 proc iml;

  use rules;  read all var {k m a b};

 /********************************************************/
 /*  Input the Standarized Shift,                        */
 /*  the Upper Limit for the Calculation of ARL and Std, */
 /*  and the Upper Limit for the CDF Graph               */
 /********************************************************/

  shf    = 0.6;
  uppcal = 0.99999;
  uppgra = 0.99;

 * Find the Regions of the Chart ;

  r  = -9//a//b//9;
  r  = unique(r);
  nt = nrow(k);
  nr = ncol(r)-1;

 /***************************************************/
 /*   Find the Length of the State Vector           */
 /*   and the pointers to the End of the Subvector  */
 /*   Associated with Each Runs Rule                */
 /***************************************************/

  nv = sum(m)-nt;
  if(sum(m-k) = 0) then nv = nv+1;
  di = j(nt,1,0);
  di[1] = m[1]-1;
  do i = 2 to nt;
    if(m[i] > 1) then di[i] = di[i-1]+m[i]-1;
  end;

 /**************************************************/
 /*   Determine the Value of the Present           */
 /*   Indicator Variable by Runs Rule and Region   */
 /*   Combination                                  */
 /**************************************************/

  a1 = j(nt,1,1);
  a1 = a1*r;
  x  = j(nt,nr,0);
  do j = 1 to nr;
    x[,j] = (a <= a1[,j])#(b >= a1[,j+1]);
  end;

 * Determine the State to State Transitions by Regions ;

  ps   = j(nv,1,0);
  s    = j(nt,1,0);
  a1   = j(1,nr,1);
  nx   = j(nv,1,0);
  qqns = 2**nv-1;
  ns   = 1;
  h    = 1;
  qq   = 0;
  q    = a1;
  do while ( h <= ns );
    qh = qq[h];
    do L = 1 to nv;
      ps[L] = qh-2*int(qh/2);
      qh = int(qh/2);
    end;
    do i = 1 to nt;
      s[i] = 0;
      if (m[i] > 1) then do
        L = di[i]-m[i]+2 to di[i];
        s[i] = s[i]+ps[L];
      end;
    end;
    do j = 1 to nr;
      sg = 0;
      do i = 1 to nt;
        if (sg = 0) then do;
          if (s[i]+x[i,j] >= k[i]) then sg = 1;
          else do;
            if (m[i] > 1) then nx[di[i]-m[i]+2] = x[i,j];
            if (m[i] > 2) then do
              L = di[i]-m[i]+3 to di[i];
              nx[L] = ps[L-1];
            end;
          end;
          if (x[i,j] = 0 & m[i] > 1) then do;
            tmp = s[i]-ps[di[i]]+1;
            L = di[i];
            ck = 0;
            do while ( ck = 0 & L >= di[i]-m[i]+2 );
              if (nx[L] = 1) then do;
                ck = 1;
                if (tmp < k[i]) then do;
                  nx[L] = 0;
                  tmp = tmp-1;
                  ck = 0;
                end;
              end;
              L = L-1;
              tmp = tmp+1;
            end;
          end;
        end;
      end;
      if (sg = 0) then do;
        qh = nx[1];
        do L = 2 to nv;
          qh = qh+nx[L]*(2**(L-1));
        end;
        ck = 0;
        do L = 1 to ns;
          if (ck = 0 & qh = qq[L]) then do;
            q[h,j] = qq[L];
            ck = 1;
          end;
        end;
        if (ck = 0) then do;
          ns = ns+1;
          qq = qq//qh;
          q[h,j] = qh;
        end;
      end;
      else
        q[h,j] = qqns;
    end;
    h = h+1;
    q = q//(qqns*a1);
  end;
  ns = ns+1;
  qq = qq//qqns;

  free / q qq ns nr r shf qqns uppcal uppgra;

 /******************************************/
 /*   Sort the States in Ascending Order   */
 /*     of Their Base Two Representations  */
 /******************************************/

  a1 = q;
  a1[rank(qq),] = q;
  q  = a1;
  a1 = qq;
  a1[rank(qq),] = qq;
  qq = a1;

 * Remove any Duplicate States ;

  ck = 1;
  do while (ck = 1);
    ck = 0;
    i  = 1;
    do while (i <= ns-1);
      j = i+1;
      do while (j <= ns);
        if (q[i,] = q[j,]) then do;
            ck = 1;
            q  = q[remove(1:ns,j),];
            a  = (q = (0*q+qq[j]));
            q  = (q#(1-a))+(qq[i]#a);
            qq = qq[remove(1:ns,j),];
            j  = j-1;
            ns = ns-1;
          end;
        j = j+1;
      end;
      i = i+1;
    end;
  end;

 * Number the Next-State Transitions ;

  a1 = j(1,ns,1);
  a2 = (qq*a1)`;
  a3 = j(ns,nr,0);
  do j = 1 to nr;
    a3[,j] = ((q[,j]*a1) = a2)*((1:ns)`);
  end;
  a1 = (q < (0*q+qqns));
  q  = (ns#(1-a1))+(a3#a1);

  free qq qqns ck a2 a3;

 /**********************************************************/
 /*   Compute the Cumulative Distribution,                 */
 /*   Function, the Average Run Length (ARL) and the       */
 /*   Standard Deviation (STD) of the Run Length for a     */
 /*   Given Standarized Shift (SHF) in the Mean.           */
 /*                                                        */
 /*   The Iterative Process for the Calculation of the     */
 /*   CDF, ARL and STD stops when the CDF is > UPPCAL      */
 /*                                                        */
 /*   The Number of Recorded Points for the CDF Graph      */
 /*   Stops When CDF is > UPPGRA                           */
 /**********************************************************/

  pr = probnorm((r-shf)`);
  pr = (pr[2:nr+1,]//0)-pr;
  pr = pr[1:nr,];

  pt = j(ns,ns,0);
  do i = 1 to ns;
    do j = 1 to nr;
      pt[i,q[i,j]] = pt[i,q[i,j]]+pr[j];
    end;
  end;

  a1  = pt[1:ns-1,ns];
  pt  = pt[1:ns-1,1:ns-1];
  cdf = a1[1];
  acu = a1[1];
  n   = 1;
  do while ( acu <= uppcal );
    a1  = pt*a1;
    acu = acu+a1[1];
    if ( acu <= uppgra ) then n=n+1;
    cdf = cdf//a1[1];
  end;

  rl  = do(1,nrow(cdf),1);
  arl = rl*cdf;
  std = sqrt(((rl##2)*cdf)-(arl*arl));
  cdf = cusum(cdf);
  rl  = rl[,1:n]`;
  cdf = cdf[1:n,];

  create par var{arl std shf}; append; close par;
  create cdf var{rl cdf};      append; close cdf;

  reset center noname pagesize = 800 linesize = 80;
  c1  = {"ARL" "STD"};
  arl = arl//std;
  print arl[rowname = c1 format = 10.6];

 quit;

 * Print Every 10th Run Length and associated CDF value ;

 title1 'Run Length CDF';
 proc print data=cdf noobs label;
  label rl='Run Length' cdf='CDF of Run Length';
  where mod(rl,10)=0;
  run;

 /******************************************************/
 /*   Annotate Dataset that will Superimpose a Table   */
 /*  Containing the Mean and Standard Deviation of the */
 /*  Run Length                                        */
 /******************************************************/

  data anno;
    set par;
    length text function style color $ 8;
    shf  = put(shf,7.3);
    arl  = put(arl,7.3);
    std  = put(std,7.3);

    xsys = '1';
    ysys = '1';
    hsys = '1';

    function = 'move';     x = 1;        y = 99;    output;

    style    = 'solid';
    color    = 'blue';
    when     = 'A';
    xsys     = 'A';
    ysys     = 'A';
    hsys     = '6';

    function = 'poly';     x = 0;        y =  0;    output;
    color    = 'white';
    function = 'polycont'; x = 0;        y = -7;    output;
    function = 'polycont'; x = 17;       y = -7;    output;
    function = 'polycont'; x = 17;       y =  0;    output;

    xsys = '1';
    ysys = '1';
    hsys = '1';

    function = 'move';     x = 1;        y = 99;    output;
    function = 'cntl2txt';                          output;

    color    = 'white';
    line     = 1;
    style    = 'none';
    size     = 1;

    position = 'c';

    xsys = 'A';
    ysys = 'A';
    hsys = '6';

    x = 1;
    function = 'label';  y = -2; text = 'SHIFT';    output;
    function = 'cntl2txt';                          output;
    function = 'label';  y = -4; text = 'ARL';      output;
    function = 'cntl2txt';                          output;
    function = 'label';  y = -6; text = 'STD';      output;

    x = 7;
    function = 'cntl2txt';                          output;
    function = 'label';  y = -2; text = '=';        output;
    function = 'cntl2txt';                          output;
    function = 'label';  y = -4; text = '=';        output;
    function = 'cntl2txt';                          output;
    function = 'label';  y = -6; text = '=';        output;

    position = 'a';

    x = 16;
    function = 'cntl2txt';                          output;
    function = 'label';  y = -2; text = shf;        output;
    function = 'cntl2txt';                          output;
    function = 'label';  y = -4; text = arl;        output;
    function = 'cntl2txt';                          output;
    function = 'label';  y = -6; text = std;        output;
  run;

 /*************************************************************/
 /*   Display Graph of the Run Length Cumulative Distribution */
 /*   Function and the Annotate Dataset                       */
 /*************************************************************/

  title1  'Run Length Cumulative Distribution Function';
  symbol1 color = white width = 2 interpol = join value = none
          line  = 1 ;
  axis1   width = 1 label = ('Run Length');
  axis2   width = 1 order = (0.0 to 1.0 by 0.1)
          label = (justify = c a=-90 r=90 'Probability');

  proc gplot data = cdf;
    plot cdf*rl /
    frame
    annotate = anno
    vminor   = 0
    hminor   = 0
    haxis    = axis1
    vaxis    = axis2
    vref     = 0.2 to 0.8 by 0.2
    lvref    = 3
    cvref    = white
    cframe   = gray
    nolegend;
  run;
  quit;

goptions reset=all;