ARL For Combined Individuals & Moving Range

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: SHWARL1                                             */
 /*   TITLE: ARL For Combined Individuals & Moving Range         */
 /* PRODUCT: QC                                                  */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: Shewhart Charts, Average Run Lengths,               */
 /*   PROCS: IML                                                 */
 /*    DATA:                                                     */
 /*                                                              */
 /*   NOTES:                                                     */
 /*    MISC:                                                     */
 /*                                                              */
 /*     REF: Crowder S. V. (1987), "A Program for the            */
 /*          Computation of ARL for Combined Individual          */
 /*          Measurement and Moving Range Charts"                */
 /*          Journal of Quality Technology, Vol. 19, pg. 103-106 */
 /*                                                              */
 /*          Crowder S. V. (1987), "Computation of ARL for       */
 /*          Combined Individual Measurement and Moving          */
 /*          Range Charts"                                       */
 /*          Journal of Quality Technology, Vol. 19, pg. 98-102  */
 /*                                                              */
 /****************************************************************/
 options ps=60 ls=80 nodate;

 proc iml;

   reset noname;

 /*  Standard Normal Density Function  */

   start p(x);
     y = 0.39894228#exp(-0.5#x#x);
     return(y);
   finish;

 /****************************************************************/
 /*   Input Parameters                                           */
 /*                                                              */
 /*   M = Control Limit Multiple for Individual Measurements     */
 /*   R = Control Limit Multiple for Moving Range Chart          */
 /*   D = Vector or Standarized Mean Shifts                      */
 /*   H = Subinterval Width for the Trapezoidal Quadrature       */
 /*       (Choose H so that M and R can be expressed as positive */
 /*       integer multiples of H)                                */
 /****************************************************************/

   m   = 2;
   r   = 3;
   d   = do(0,2.0,0.25)`;
   h   = 0.05;

   n   = 2*int(m/h);
   nn  = n+1;
   a   = j(nn,nn,0);
   b   = j(nn,1,-1);
   arl = j(nrow(d),1,0);
   ai  = j(nn,1,0);
   bi  = ai;
   n1  = ai;
   n2  = ai;

 /********************************************************/
 /*   Apply Trapezoidal Rule to replace equation (1) in  */
 /*   Crowder (1987) with a system of linear algebraic   */
 /*   equations                                          */
 /********************************************************/

   do i = 1 to nn;
     bi[i] = min( m,-m+h*(i-1)+r);
     ai[i] = max(-m,-m+h*(i-1)-r);
     n1[i] = int((ai[i]+m)/h)+1;
     n2[i] = int((bi[i]+m)/h)+1;
   end;

   do k = 1 to nrow(d);
     aux = d[k];
     a   = 0#a;
     do i = 1 to nn;
       do j = 1 to nn;
         if (j > n1[i] & j < n2[i]) then
           a[i,j] = h*p(ai[i]+h*(j-n1[i])-d[k]);
         else if (j = n1[j]) then
           a[i,j] = h/2*p(ai[i]-d[k]);
         else if (j = n2[i]) then
           a[i,j] = h/2*p(bi[i]-d[k]);
       end;
     end;

     a = a - i(nn);
     if (det(a) = 0) then
        print "Zero Determinant for linear system";
     else do;
         x = solve(a,b);

 /***********************************************************/
 /*   Approximate ARL in equation (2) of Crowder (1987) by  */
 /*   the solution of the linear system A*X=B               */
 /***********************************************************/

         arl[k] = (h/2)*(p(-m-d[k])*x[1]+p(m-d[k])*x[nn]);
         arg    = p((h#do(1,n-1,1))`-m-d[k]);
         x      = x[2:n];
         arl[k] = arl[k]+((arg`)*(h#x))+1;
       end;
   end;

 * Format Printed Output ;

  colnam = { "|Mu-Muo|/Sigma" "      ARL " };
  aux    = d||arl;
  print "     ARL for Combined Individuals and Moving Range Chart";
  print "         M =" m[format = 4.1] "and R =" r[format = 4.1];
  print colnam;
  print aux[format = 10.2];

 quit;

 title1;