nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

JMP equivalent to enclosed SAS script for evaluating Run Length distribution?

jerry_cooper

Hi, all,

I have a customer asking if there is an existing JMP script or process for evaluating the Run Length Distribution of a Shewhart Control Chart with Supplementary Runs Rules? The enclosed SAS script does this, anyone know of a JMP equivalent?

Thanks!

Jerry

SAS script for RL distribution.zip

Accepted answers

ian_jmp

FYI, I translated this to JSL yesterday, and it was an interesting exercise. I didn't make a UI since I'm assuming that (given the specialist nature of the topic) anyone who needs it would be happy to make small changes to the code. It's had no testing apart from trivial cases.

It's another good illustration of the 'flaw of averages', and the dangers of selecting 'good' charts using just the ARL.

Cheers,

Ian

Shewhart Run Lenth Distribution.jsl

All comments

Jeff_Perkinson

Here's the SAS Pogram that Jerry references:


 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: SHWARL2                                             */
 /*   TITLE: ARL With Supplementary Run Rules (SHWARL2)          */
 /* 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;

ian_jmp

It's another good illustration of the 'flaw of averages', and the dangers of selecting 'good' charts using just the ARL.

Cheers,

Ian

Shewhart Run Lenth Distribution.jsl

jerry_cooper

Ian,

Thank you very much for doing this! Unfortunately, I can't seem to access the script... my browser locks up when trying to open or save it to my computer. Also, how would you suggest making this available to the user (and others) that requested it? Is it something you'd be willing to put on the File Exchange?

Thanks!

Jerry

jerry_cooper

Update: for the record, this appears to have been a transient problem for me... I am now able to download the file.