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;