C Program Name : CUSUM C This program computes the probability function, C the cumulative distribution function, and C percentiles of the run length of a CUSUM chart C with parameters k and h. It also computes the ARL C and plots the probability function and cumulative C distribution function. C PARAMETER (NMAX1=3000,NPMAX=999,NPER1=10) DOUBLE PRECISION H,K,DELTA,U,ARL,PPP,PROB(NPMAX) &,PUSER(NPER1),CUM,P(25),W(25),PR(25) DIMENSION PROBN(NMAX1),CDF(NMAX1),RUN(NMAX1), &NQUANT(NPMAX),NUSER(NPER1) CHARACTER*1 GO CHARACTER*4 ITYPE CHARACTER*30 YLABEL CHARACTER*10 XLABEL CHARACTER*40 FFMMTT 100 NTIME=0 NRUN=0 NTYPE=1 ITYPE='S_0=' WRITE(*,400) READ(*,*) K IF(K.LT.0) THEN NTYPE=-1 ITYPE='T_0=' ENDIF K=ABS(K) WRITE(*,500) READ(*,*) H H=ABS(H) 200 INOTE=0 WRITE(*,600) READ(*,*) DELTA DELTA=ABS(DELTA) WRITE(*,700) READ(*,*) U U=ABS(U) NTIME=NTIME+1 IF(NTIME.EQ.1) THEN WRITE(*,800) NPER1 WRITE(*,900) READ(*,*) NPER IF(NPER.GT.0) THEN WRITE(*,1000) NPER WRITE(*,1100) READ(*,*) (PUSER(I),I=1,NPER) DO 300 I=1,NPER NUSER(I)=IDINT(1000.D0*PUSER(I)) 300 CONTINUE WRITE(*,1200) NMAX1 READ(*,*) NMAX ENDIF ENDIF C C Compute ARL C WRITE(*,*) IF(NPER.EQ.0) THEN WRITE(*,1300) ELSE WRITE(*,1400) ENDIF CALL SETUP(H,P,W) CALL CUSARL(K,DELTA,P,W,ARL,U) C C Compute probability function C IF(NPER.GT.0) THEN CUM=0.D0 DO 1600 I=1,NPMAX PROB(I)=DBLE(I)/1000.D0 NQUANT(I)=0 1600 CONTINUE KOUNT=1 DO 1800 N=1,NMAX IF(N.GE.NMAX) THEN WRITE(*,1500) NMAX NRUN=NMAX INOTE=1 GOTO 1900 ENDIF CALL CUSP(K,H,N,DELTA,P,W,U,PR,PPP) IF(PPP.LT.1.D-20) PPP=0.D0 CUM=CUM+PPP RUN(N)=REAL(N) PROBN(N)=PPP CDF(N)=CUM 1700 IF(KOUNT.LE.NPMAX.AND.CUM.GT.PROB(KOUNT)) & THEN NQUANT(KOUNT)=N KOUNT=KOUNT+1 GOTO 1700 ENDIF IF(CUM.GT.PUSER(NPER)) THEN NRUN=N GOTO 1900 ENDIF 1800 CONTINUE ENDIF C C Print menu C 1900 WRITE(*,2000) WRITE(*,2100) WRITE(*,2200) WRITE(*,2300) WRITE(*,2400) WRITE(*,2500) WRITE(*,2600) READ(*,*) IPICK IF(IPICK.EQ.1) THEN WRITE(*,2700) K*NTYPE WRITE(*,2800) H*NTYPE WRITE(*,2900) U*NTYPE WRITE(*,3000) DELTA*NTYPE WRITE(*,3100) ARL IF(NPER.GT.0) THEN IF(INOTE.EQ.0) THEN WRITE(*,3200) WRITE(*,3300) (PUSER(I),I=1,NPER) WRITE(*,3400) (NQUANT(NUSER(I)),I=1,NPER) ELSE DO 3700 I=1,NPER IF(NQUANT(NUSER(I)).EQ.0) THEN NOT0=I-1 GOTO 3800 ENDIF 3700 CONTINUE NOT0=NPER 3800 CONTINUE NOT0P1=NOT0+1 WRITE(*,3300) (PUSER(I),I=1,NOT0P1) IF(NOT0.EQ.0) THEN WRITE(*,3500) NMAX ELSE FFMMTT='(1X,'//CHAR(NOT0+ICHAR('0'))// & '(I5,1X),'//''' >'''//',I4)' WRITE(*,FMT=FFMMTT) (NQUANT(NUSER(I)), & I=1,NOT0),NMAX ENDIF ENDIF ENDIF WRITE(*,3600) READ(*,4300) GO C C Plot probability function C ELSEIF(IPICK.EQ.2) THEN IF(NRUN.LT.2) THEN WRITE(*,3900) ELSE ISCALE=0 XLABEL='RUN LENGTH' YLABEL='PROBABILITY' CALL PLOT(NRUN,RUN,PROBN,ISCALE,XLABEL, & YLABEL) IF(CUM.LT.0.5D0) THEN WRITE(*,4000) K*NTYPE,H*NTYPE,ITYPE, & U*NTYPE,ARL ELSE WRITE(*,4100) K*NTYPE,H*NTYPE,ITYPE, & U*NTYPE,ARL,NQUANT(500) ENDIF WRITE(*,4200) DELTA*NTYPE WRITE(*,3600) READ(*,4300) GO ENDIF C C Plot cumulative distribution function C ELSEIF(IPICK.EQ.3) THEN IF(NRUN.LT.2) THEN WRITE(*,3900) ELSE ISCALE=1 XLABEL='RUN LENGTH' YLABEL='CUMULATIVE PROBABILITY' CALL PLOT(NRUN,RUN,CDF,ISCALE,XLABEL,YLABEL) IF(CUM.LT.0.5D0) THEN WRITE(*,4000) K*NTYPE,H*NTYPE,ITYPE, & U*NTYPE,ARL ELSE WRITE(*,4100) K*NTYPE,H*NTYPE,ITYPE, & U*NTYPE,ARL,NQUANT(500) ENDIF WRITE(*,4200) DELTA*NTYPE WRITE(*,3600) READ(*,4300) GO ENDIF ELSEIF(IPICK.EQ.4) THEN GOTO 100 ELSEIF(IPICK.EQ.5) THEN GOTO 200 ELSEIF(IPICK.EQ.9) THEN GOTO 9999 ENDIF GOTO 1900 9999 CONTINUE 400 FORMAT(/1X,'ENTER THE CONTROL CHART', &' PARAMETER K ') 500 FORMAT(/1X,'ENTER THE CONTROL CHART', &' LIMIT H ') 600 FORMAT(/1X,'ENTER THE SHIFT IN PROCESS MEAN ') 700 FORMAT(/1X,'ENTER THE INITIAL CUSUM VALUE ') 800 FORMAT(/1X,'ENTER THE NUMBER OF PERCENTAGE', &' POINTS DESIRED (0-',I2,').') 900 FORMAT(1X,'(AN ENTRY OF 0 MEANS THAT ONLY', &' THE ARL IS COMPUTED.) ') 1000 FORMAT(/1X,'ENTER THE ',I2,' CUMULATIVE', &' PROBABILITIES IN INCREASING') 1100 FORMAT(1X,'ORDER (SEPARATED BY SPACES) ') 1200 FORMAT(/1X,'ENTER AN UPPER BOUND (1-',I4,') ON', &' THE VALUE OF THE MAXIMUM',/' PERCENTAGE', &' POINT DESIRED ') 1300 FORMAT(1X,'COMPUTING ARL') 1400 FORMAT(1X,'COMPUTING ARL AND PROBABILITY', &' FUNCTION, PLEASE WAIT') 1500 FORMAT(//1X,I5,'- MAXIMUM', &' PERCENTAGE POINT EXCEEDED') 2000 FORMAT(/1X,'MENU') 2100 FORMAT(1X,'1 PRINT PERCENTAGE POINTS AND ARL') 2200 FORMAT(1X,'2 PLOT PROBABILITY FUNCTION') 2300 FORMAT(1X,'3 PLOT CUMULATIVE DISTRIBUTION', &' FUNCTION') 2400 FORMAT(/1X,'4 NEW CUSUM CONTROL CHART') 2500 FORMAT(1X,'5 NEW SHIFT IN MEAN OR NEW INITIAL', &' CUSUM VALUE') 2600 FORMAT(1X,'9 STOP'/) 2700 FORMAT(1X,' K = ',F8.5,/) 2800 FORMAT(1X,' H = ',F8.5,/) 2900 FORMAT(1X,' U = ',F8.5,/) 3000 FORMAT(1X,'SHIFT IN MEAN = ',F8.5,/) 3100 FORMAT(1X,' ARL = ',F8.2,//) 3200 FORMAT(1X,'PERCENTAGE POINTS',/) 3300 FORMAT(1X,10(F5.3,1X)) 3400 FORMAT(1X,10(I5,1X)) 3500 FORMAT(1X,'>',I4) 3600 FORMAT(1X,//,1X,'PRESS RETURN KEY TO CONTINUE') 3900 FORMAT(1X,//,1X,'PLOT NOT CONSTRUCTED') 4000 FORMAT(1X,'K=',F6.3,3X,'H=',F6.3, &3X,A4,F5.2,3X,'ARL=',F7.2) 4100 FORMAT(1X,'K=',F6.3,3X,'H=',F6.3, &3X,A4,F5.2,3X,'ARL=', &F7.2,3X,'MEDIAN=',I4) 4200 FORMAT(1X,'SHIFT IN MEAN=',F6.3) 4300 FORMAT(A) STOP END C SUBROUTINE SETUP(H,PT,WT) C C Transforms the Gaussian quadrature points C from the interval (-1,1) to the interval C (0,H). C DOUBLE PRECISION H,P(25),W(25),PT(25),WT(25) DATA P/0.9951872199970213D0,0.9747285559713095D0 &,0.9382745520027327D0,0.8864155270044010D0, &0.8200019859739029D0,0.7401241915785543D0, &0.6480936519369755D0,0.5454214713888395D0, &0.4337935076260451D0,0.3150426796961634D0, &0.1911188674736163D0,0.0640568928626056D0, &0,0,0,0,0,0,0,0,0,0,0,0,0/ DATA W/0.0123412297999872D0,0.0285313886289337D0 &,0.0442774388174198D0,0.0592985849154368D0 &,0.0733464814110803D0,0.0861901615319533D0 &,0.0976186521041139D0,0.1074442701159656D0 &,0.1155056680537256D0,0.1216704729278034D0 &,0.1258374563468283D0,0.1279381953467521D0, &0,0,0,0,0,0,0,0,0,0,0,0,0/ NQUAD=24 ND2=NQUAD/2 NQP1=NQUAD+1 DO 200 I=1,ND2 P(NQP1-I)=-P(I) W(NQP1-I)=W(I) 200 CONTINUE DO 300 I=1,NQUAD WT(I)=H*W(I)/2.D0 PT(I)=H*P(I)/2.D0+H/2.D0 300 CONTINUE RETURN END C SUBROUTINE CUSP(K,H,N,DELTA,P,W,U,PR,PPP) C C Computes the probability function of the run length C DOUBLE PRECISION P(25),W(25),PR(25),PRX(25),K,H, &DELTA,U,ARG,F,PPP,DNML NQUAD=24 NQP1=NQUAD+1 IF(N.EQ.1) THEN DO 100 J=1,NQUAD ARG=H+K-P(J)-DELTA PR(J)=1.D0-DNML(ARG) 100 CONTINUE ARG=H+K-DELTA PR(NQP1)=1.D0-DNML(ARG) ENDIF IF(N.GE.2) THEN DO 200 J=1,NQP1 PRX(J)=PR(J) 200 CONTINUE DO 400 J=1,NQUAD ARG=K-P(J)-DELTA PR(J)=PRX(NQP1)*DNML(ARG) DO 300 I=1,NQUAD ARG=P(I)+K-P(J) PR(J)=PR(J)+W(I)*PRX(I)*F(ARG,DELTA,1.D0) 300 CONTINUE 400 CONTINUE ARG=K-DELTA PR(NQP1)=PRX(NQP1)*DNML(ARG) DO 500 I=1,NQUAD ARG=P(I)+K PR(NQP1)=PR(NQP1)+W(I)*PRX(I)* & F(ARG,DELTA,1.D0) 500 CONTINUE ENDIF IF(N.EQ.1) THEN ARG=H+K-U-DELTA PPP=1.D0-DNML(ARG) ENDIF IF(N.GE.2) THEN ARG=K-U-DELTA PPP=PRX(NQP1)*DNML(ARG) DO 600 I=1,NQUAD ARG=P(I)+K-U PPP=PPP+W(I)*PRX(I)*F(ARG,DELTA,1.D0) 600 CONTINUE ENDIF RETURN END C DOUBLE PRECISION FUNCTION F(X,XMU,SIGMA) C C Computes the normal density function C DOUBLE PRECISION X,XMU,SIGMA F=(3.989422804014327D-1/SIGMA)* &DEXP(-0.5D0*(X-XMU)*(X-XMU)/(SIGMA*SIGMA)) RETURN END C SUBROUTINE PLOT(N,XVAL,YVAL,ISCALE,XLABEL, &YLABEL) C C Source: Jensen, K. L., Crowder, S. V. and C Vardeman, S. B. (1988). "An Interactive C Probability Plotting Program". JQT 20, C pp. 196-210. C DIMENSION XVAL(N),YVAL(N),X(9) REAL LO INTEGER PRNFLG(75),COUNT CHARACTER*1 PRNCHR(2) CHARACTER*30 YLABEL CHARACTER*10 XLABEL DATA PRNCHR/' ','*'/ XMIN=1.E+10 XMAX=-1.E+10 DO 100 I=1,N IF(XVAL(I).GT.XMAX) XMAX=XVAL(I) IF(XVAL(I).LT.XMIN) XMIN=XVAL(I) 100 CONTINUE DX=(XMAX-XMIN)/64. YMIN=1.E+10 YMAX=-1.E+10 DO 120 I=1,N IF(YVAL(I).GT.YMAX) YMAX=YVAL(I) IF(YVAL(I).LT.YMIN) YMIN=YVAL(I) 120 CONTINUE IF(ISCALE.EQ.1) THEN YMAX=1. YMIN=0. ENDIF DY=(YMAX-YMIN)/16. WRITE(*,1000) YLABEL LO=YMAX COUNT=3 DO 300 I=1,17 DO 200 J=1,65 PRNFLG(J)=1 200 CONTINUE UP=LO LO=UP-DY DO 220 K=1,N IF((YVAL(K).GT.LO).AND.(YVAL(K).LE.UP)) THEN J=INT((XVAL(K)-XMIN)/DX)+1 PRNFLG(J)=2 ENDIF 220 CONTINUE IF(COUNT.EQ.3) THEN WRITE(*,1001) UP,(PRNCHR(PRNFLG(J)),J=1,65) COUNT=0 ELSE WRITE(*,1002) (PRNCHR(PRNFLG(J)),J=1,65) COUNT=COUNT+1 ENDIF 300 CONTINUE WRITE(*,1003) DX=(XMAX-XMIN)/8. DO 320 I=1,9 X(I)=XMIN+FLOAT(I-1)*DX 320 CONTINUE WRITE(*,1004) (INT(X(I)),I=1,9),XLABEL RETURN 1000 FORMAT(/,1X,A30,/,9X,'|') 1001 FORMAT(1X,F7.5,' +',1X,65A1) 1002 FORMAT(9X,'|',1X,65A1) 1003 FORMAT(10X,'-+-------+-------+-------+----', &'---+-------+-------+-------+-------+----') 1004 FORMAT(6X,9(I6,2X),/,62X,A10) END C SUBROUTINE CUSARL(K,DELTA,P,W,ARL,U) C C Computes the ARL of a CUSUM chart C DOUBLE PRECISION K,DELTA,ARL,U,F,A(25,25), &B(25),X(25),WK(25),P(25),W(25),DNORM,DNML,ARG INTEGER IPIVOT(25) NQP1=25 NQUAD=NQP1-1 C C Replace integral equation with a system of C linear equations AX=B C DO 200 I=1,NQP1 B(I)=-1.D0 DO 100 J=1,NQP1 IF(I.LE.NQUAD.AND.J.EQ.I) THEN ARG=P(J)+K-P(I) A(I,J)=W(J)*F(ARG,DELTA,1.D0)-1.D0 ELSEIF(I.LE.NQUAD.AND.J.LE.NQUAD) THEN ARG=P(J)+K-P(I) A(I,J)=W(J)*F(ARG,DELTA,1.D0) ELSEIF(I.LE.NQUAD.AND.J.EQ.NQP1) THEN DNORM=K-P(I)-DELTA A(I,J)=DNML(DNORM) ELSEIF(I.EQ.NQP1.AND.J.LT.I) THEN ARG=P(J)+K A(I,J) = W(J)*F(ARG,DELTA,1.D0) ELSEIF(I.EQ.NQP1.AND.J.EQ.I) THEN DNORM=K-DELTA A(I,J) = DNML(DNORM)-1.D0 ENDIF 100 CONTINUE 200 CONTINUE CALL FACTOR(A,NQP1,WK,IPIVOT,IFLAG) IF(IFLAG.EQ.0) THEN WRITE(*,300) STOP ENDIF CALL SUBST(A,IPIVOT,B,NQP1,X) C C obtain ARL(u) C ARL=1.D0 DNORM = K-U-DELTA ARL = ARL+DNML(DNORM)*X(NQP1) DO 400 J=1,NQUAD ARG=P(J)+K-U ARL=ARL+W(J)*F(ARG,DELTA,1.D0)*X(J) 400 CONTINUE RETURN 300 FORMAT(5X,'ZERO DETERMINANT FOR LINEAR SYSTEM') END C SUBROUTINE FACTOR(W1,N,D1,IPIVOT,IFLAG) C C Finds the solution of a system of linear equations C DOUBLE PRECISION D1(N),W1(N,N),AWIKOD,COLMAX, &RATIO,ROWMAX,TEMP,ZERO DIMENSION IPIVOT(N) IFLAG=1 ZERO=1.D-12 DO 200 I=1,N IPIVOT(I)=I ROWMAX=0.D0 DO 100 J=1,N ROWMAX=DMAX1(ROWMAX,DABS(W1(I,J))) 100 CONTINUE IF(DABS(ROWMAX).LT.ZERO) THEN IFLAG=0 ROWMAX=1.D0 ENDIF D1(I)=ROWMAX 200 CONTINUE IF(N.LE.1) RETURN N1=N-1 DO 700 K=1,N1 COLMAX=DABS(W1(K,K))/D1(K) ISTAR=K K1=K+1 DO 300 I=K1,N AWIKOD=DABS(W1(I,K))/D1(K) IF(AWIKOD.GT.COLMAX) THEN COLMAX=AWIKOD ISTAR=I ENDIF 300 CONTINUE IF(DABS(COLMAX).LT.ZERO) THEN IFLAG=0 ELSE IF(ISTAR.GT.K) THEN IFLAG=-IFLAG I=IPIVOT(ISTAR) IPIVOT(ISTAR)=IPIVOT(K) IPIVOT(K)=I TEMP=D1(ISTAR) D1(ISTAR)=D1(K) D1(K)=TEMP DO 400 J=1,N TEMP=W1(ISTAR,J) W1(ISTAR,J)=W1(K,J) W1(K,J)=TEMP 400 CONTINUE ENDIF K2=K+1 DO 600 I=K2,N W1(I,K)=W1(I,K)/W1(K,K) RATIO=W1(I,K) K3=K+1 DO 500 J=K3,N W1(I,J)=W1(I,J)-RATIO*W1(K,J) 500 CONTINUE 600 CONTINUE ENDIF 700 CONTINUE IF(DABS(W1(N,N)).LT.ZERO) IFLAG=0 RETURN END C SUBROUTINE SUBST(W1,IPIVOT,B,N,X2) C C Sources: Conte, S. D. and de Boor, C. (1980). C Elementary Numerical Analysis. McGraw-Hill, C New York, NY. C Crowder, S. V. (1987). "Average Run Lengths C of Exponentially Weighted Moving Average C Control Charts". JQT 19, pp. 161-164. C DOUBLE PRECISION B(N),W1(N,N),X2(N),SUM DIMENSION IPIVOT(N) IF(N.LE.1) THEN X2(1)=B(1)/W1(1,1) RETURN ENDIF IP=IPIVOT(1) X2(1)=B(IP) DO 200 I=2,N SUM=0.D0 I1=I-1 DO 100 J=1,I1 SUM=W1(I,J)*X2(J)+SUM 100 CONTINUE IP=IPIVOT(I) X2(I)=B(IP)-SUM 200 CONTINUE X2(N)=X2(N)/W1(N,N) I2=N-1 DO 400 ISTEP=1,I2 I=N-ISTEP SUM=0.D0 I3=I+1 DO 300 J=I3,N SUM=W1(I,J)*X2(J)+SUM 300 CONTINUE X2(I)=(X2(I)-SUM)/W1(I,I) 400 CONTINUE RETURN END C DOUBLE PRECISION FUNCTION DNML(X) C C Computes the cumulative distribution function C Pr[Y<=x] of a random variable Y having a C standard normal distribution. C C Source: Craig, R. J. (1984). "Normal Family C Distribution Functions: FORTRAN and BASIC C Programs". JQT 16, pp. 232-236. C DOUBLE PRECISION X,Y,S,RN,ZERO,ONE,ERF,SQRT2,PI DATA SQRT2,ONE/1.414213562373095D0,1.D0/ DATA PI,ZERO/3.141592653589793D0,0.D0/ IF(DABS(X).GT.8.3D0) GO TO 200 Y=X/SQRT2 IF(X.LT.ZERO) Y=-Y S=ZERO DO 100 N=1,37 RN=DBLE(N) S=S+DEXP(-RN*RN/25)/N*DSIN(2*N*Y/5) 100 CONTINUE S=S+Y/5 ERF=2*S/PI DNML=(ONE+ERF)/2 IF(X.LT.ZERO) DNML=(ONE-ERF)/2 RETURN 200 IF(X.LT.0.D0) DNML=ZERO IF(X.GT.0.D0) DNML=ONE RETURN END