Cusum ARLs Using Markov Chain
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: CUSARLM */
/* TITLE: Cusum ARLs Using Markov Chain */
/* PRODUCT: QC */
/* SYSTEM: ALL */
/* KEYS: SQC */
/* PROCS: IML */
/* DATA: */
/* */
/* REF: Hawkins, D. M. (1992). "Evaluation of Average Run */
/* Lengths of Cumulative Sum Charts for an Arbitrary */
/* Data Distribution." Communications in Statistics-- */
/* Simulation. 21(4) 1001-1020. */
/* */
/* NOTES: This example provides a SAS/IML program to compute */
/* the ARL of cusum schemes using a Markov chain */
/* method. The program is a translation of a program */
/* due to Hawkins (1992). */
/* */
/* To support distributions other than the normal, */
/* replace the two calls to probnorm and make */
/* modifications to denrat as described by Hawkins. */
/* */
/****************************************************************/
proc iml;
start arlcomp;
/* input parameters */
di = 5.0; /* decision interval h (H=) */
ref = 0.5; /* reference value k (K=) */
denrat = 0.0; /* common denominator */
arg2 = 0.0; /* distribution parameter (DELTA=) */
winsor = 10.0; /* winsorization value (W=) */
eps = 0.001;
/* initialize output parameters */
regarl = 0.0;
firarl = 0.0;
/* constants */
esterr = 0.0;
maxhlf = 7;
maxpow = 2##maxhlf; /* 2**maxhlf */
toler = 1.0e-12; /* subtractive cancellation */
/* working matrices */
tran = j( maxpow + 2, maxpow + 2, 0 );
answer = j( maxpow + 2, 1, 0 );
copy = j( maxpow + 2, maxpow + 2, 0 );
phitab = j( 4#maxpow + 5, 1, 0 );
richar1 = j( maxhlf, maxhlf, . );
richar2 = j( maxhlf, maxhlf, . );
/* check parameters for feasibility */
if ( di <= 0.0 )
then do;
print , "ERROR: Nonpositive decision interval";
stop;
end;
if ( winsor <= ref )
then do;
print , "ERROR: Winsorization value less than reference value";
stop;
end;
if ( denrat < 0 )
then do;
print , "ERROR: Common denominator is less than zero.";
stop;
end;
/* are the data continuous? */
if ( denrat = 0 ) then realx = 1; else realx = 0;
/* other initialization */
mbig = 1;
idenrn = 0;
/* subscript adjustments */
phadj = 2 # maxpow + 3;
if ( realx = 0 )
then do;
if ( ref <= 0.0 )
then knumer = int( ref*denrat - 0.5 );
else
knumer = int( ref#denrat + 0.5 );
idenrn = int( denrat );
mbig = int( di#denrat - 0.5 );
iwintg = int( winsor#denrat + 0.5 );
mb1 = mbig + 1;
if ( mbig >= maxpow )
then do;
print , "ERROR: di#denrat exceeds maxpow";
stop;
end;
deviat = abs( knumer - ref#denrat ) +
abs( mbig - di#denrat + 1 ) +
abs( iwintg - winsor#denrat ) ;
if ( deviat > 0.01 )
then do;
print , "ERROR: di, ref or winsor not a multiple of 1/denrat";
stop;
end;
end;
answer[ 1 ] = 0;
/* set up coefficient matrix of linear equations */
do mloop = 1 to maxhlf;
if ( realx = 1 )
then do;
mbig = 2 # mbig;
mb1 = mbig + 1;
delta = di / mbig;
mtwo = 2 # mbig;
lastnx = mtwo + 2;
fact = 4.0;
end;
else lastnx = mbig;
do i= -lastnx - idenrn to lastnx;
iadj = i + phadj;
if ( realx = 1 )
then do;
fact = 5.0 - fact;
arg = 0.5 # delta # i + ref;
if ( arg > winsor ) then
phitab[ iadj ] = fact / 6.0;
else
phitab[ iadj ] = probnorm(arg - arg2) # fact / 6.0;
end;
else do;
iarg = i + knumer;
if ( mod( iarg, idenrn ) = 0 )
then do;
arg = iarg / idenrn;
phitab[ iadj ] = 1.0;
if ( arg < winsor ) then
phitab[ iadj ] = probnorm(arg - arg2);
end;
else do;
phitab[ iadj ] = phitab[ iadj - 1 ];
if ( iarg >= iwintg ) then
phitab[ iadj ] = 1.0;
end;
end;
end;
do j = 1 to mbig;
if ( realx = 1 )
then do;
inx = 2 # j;
tran[ j+1, 1 ] = phitab[ -inx + phadj ] +
phitab[ -inx + 1 + phadj ] +
phitab[ -inx + 2 + phadj ] ;
tran[ 1, j+1 ] = 6.0 # ( phitab[ inx + phadj ] -
phitab[ inx - 2 + phadj ] );
tran[ mb1 + 1, j + 1 ] = 0.0;
do i = 1 to mbig;
inx = 2 # ( i - j );
tran[ j+1, i+1 ] = ( phitab[ inx + 2 + phadj ] -
phitab[ inx + phadj ] ) +
( phitab[ inx + 1 + phadj ] -
phitab[ inx - 1 + phadj ] ) +
( phitab[ inx + phadj ] -
phitab[ inx - 2 + phadj ] ) ;
end;
end;
else do;
tran[ j+1, 1 ] = phitab[ -j + phadj ];
tran[ 1, j+1 ] = phitab[ j + phadj ] - phitab[ j-1 + phadj ];
do i = 1 to mbig;
tran[ j+1, i+1 ] = phitab[ i-j + phadj ] -
phitab[ i-j - 1 + phadj ] ;
end;
end;
end;
tran[ 1, 1 ] = phitab[ phadj ];
if ( realx = 1 ) then tran[ 1, 1 ] = 6.0 # phitab[ phadj ];
tran[ mb1 + 1, mb1 + 1 ] = 1.0;
tran[ mb1 + 1, 1 ] = 0.0;
do j = 0 to mb1;
copy[ j + 1, mb1 + 1 ] = 1.0;
do i = 0 to mbig;
copy[ j + 1, i + 1 ] = -tran[ j + 1, i + 1 ];
if ( i = j ) then
copy[ j + 1, i + 1 ] = copy[ j + 1, i + 1 ] + 1.0;
end;
end;
/* solve the equations */
result = sweep( copy, 1 : mbig + 1 );
answer = result[ , mb1 + 1 ];
/* exit loop if data are integer */
if ( realx = 0 )
then do;
regarl = answer[ 1 ];
firarl = answer[ ( mbig + 1 ) / 2 + 1 ];
print , "Algorithm terminated for integer case";
print , regarl;
print , firarl;
stop;
end;
/* for a head start a*di other than di/2, set the subscript */
/* above to the integer part of a*(mbig+1) */
richar1[ mloop, 1 ] = answer[ 1 ];
richar2[ mloop, 1 ] = ( answer[ mbig / 2 + 1 ] +
answer[ mbig / 2 + 2 ] ) / 2.0 ;
/* for a head start a*di other than di/2, use linear inter- */
/* polation of the table entries bracketing a*mbig */
/* perform Richardson extrapolation */
pow = 1.0;
do jl = 2 to mloop;
pow = pow # 4;
richar1[ mloop, jl ] = ( pow # richar1[ mloop, jl-1 ] -
richar1[ mloop-1, jl-1 ] ) / ( pow - 1 );
richar2[ mloop, jl ] = ( pow # richar2[ mloop, jl-1 ] -
richar2[ mloop-1, jl-1 ] ) / ( pow - 1 );
end;
if ( mloop ^= 1 )
then do;
oacce1 = richar1[ mloop - 1, mloop - 1 ];
oacce2 = richar2[ mloop - 1, mloop - 1 ];
acce1 = richar1[ mloop, mloop ];
acce2 = richar2[ mloop, mloop ];
esterr = max( abs( acce1 - oacce1 ),
abs( acce2 - oacce2 ) ) / abs( acce1 );
/* test for convergence */
if ( esterr < eps )
then do;
/* diagnostic */
print , richar1, richar2;
regarl = acce1;
firarl = acce2;
print , "Algorithm converged";
print , regarl;
print , firarl;
stop;
end;
end;
end;
/* diagnostic */
print , richar1, richart2;
print , "ERROR: Convergence to requested precision not obtained";
finish arlcomp;
run arlcomp;
quit;
run;