C C Program Name : SAMPLE C THIS PROGRAM IS TAKEN FROM AN ARTICLE ENTITLED "MINIMUM COST C SAMPLING PLANS USING BAYESIAN METHODS," BY THOMAS J. LORENZEN,GENERAL C MOTORS RESEARCH PUBLICATION GMR-4300, MARCH 2,1983. C IMPLICIT REAL*8(A-H,O-Z) DIMENSION ISTO(2,500),STO(2,500) REAL*8 N,LGAMMA CHARACTER*8 ACC,REJ,ALW,ZER,PRIN1,PRIN2,PRIN3 WRITE(*,99) 99 FORMAT(////,1X,'MINIMUM COST SAMPLING PLANS',////) C C INPUTS TO PROGRAM: C C A0=COST ASSOCIATED WITH ACCEPTING A GOOD ITEM C A1=COST ASSOCIATED WITH ACCEPTING A DEFECTIVE ITEM C R0=COST ASSOCIATED WITH REJECTING A GOOD ITEM C R1=COST ASSOCIATED WITH REJECTING A DEFECTIVE ITEM C S0=COST ASSOCIATED WITH SAMPLING A GOOD ITEM C S1=COST ASSOCIATED WITH SAMPLING A DEFECTIVE ITEM C QM=MOST LIKELY PERCENT DEFECTIVE SHIPPED BY SUPPLIER C QO=LOWEST PERCENT DEFECTIVE SHIPPED BY SUPPLIER C QP=HIGHEST PERCENT DEFECTIVE SHIPPED BY SUPPLIER C N=LARGEST LOT SIZE OF INTEREST C WRITE(*,100) 100 FORMAT(/,1X,'ENTER COST OF ACCEPTING A GOOD ITEM:',/) READ(*,*) A0 WRITE(*,101) 101 FORMAT(/,1X,'ENTER COST OF ACCEPTING A DEFECTIVE:',/) READ(*,*) A1 WRITE(*,102) 102 FORMAT(/,1X,'ENTER COST OF REJECTING A GOOD ITEM:',/) READ(*,*) R0 WRITE(*,103) 103 FORMAT(/,1X,'ENTER COST OF REJECTING A DEFECTIVE ITEM:',/) READ(*,*) R1 WRITE(*,104) 104 FORMAT(/,1X,'ENTER COST OF SAMPLING A GOOD ITEM:',/) READ(*,*) S0 WRITE(*,105) 105 FORMAT(/,1X,'ENTER COST OF SAMPLING A DEFECTIVE ITEM:',/) READ(*,*) S1 WRITE(*,111) 111 FORMAT(/,1X,'WE ASSUME THAT THE SUPPLIER QUALITY CAN BE DESCRIBED 1 ',/,1X,'BY A BETA(A,B) DISTRIBUTION.',//,1X, 1 'IF YOU WISH TO SPECIFY THE PARAMETERS (A,B) OF THE BETA',/, 1 1X,'DISTRIBUTION DIRECTLY, ENTER 1 AND RETURN.',//,1X, 1 'IF YOU WISH TO SPECIFY:',//, 1 5X,'1. MOST LIKELY % DEFECTIVE SHIPPED BY SUPPLIER',/, 1 5X,'2. LOWEST % DEFECTIVE SHIPPED BY SUPPLIER',/, 1 5X,'3. HIGHEST % DEFECTIVE SHIPPED BY SUPPLIER',//, 1 1X,'THE PARAMETERS (A,B) WILL BE APPROXIMATED FOR YOU. IF SO,',/, 1 1X,'ENTER 2 AND RETURN.',/) READ(*,*) PATH IF (PATH .EQ. 1.0D0) THEN WRITE(*,112) 112 FORMAT(/,1X,'ENTER A:',/) READ(*,*) A WRITE(*,113) 113 FORMAT(/,1X,'ENTER B:',/) READ(*,*) B GO TO 110 END IF WRITE(*,106) 106 FORMAT(/,1X,'ENTER MOST LIKELY PERCENT DEFECTIVE SHIPPED BY',/, 1 1X,'SUPPLIER (10 PERCENT = .10):',/) READ(*,*) QM WRITE(*,107) 107 FORMAT(/,1X,'ENTER LOWEST PERCENT DEFECTIVE SHIPPED BY SUPPLIER:' 1 ,/) READ(*,*) QO WRITE(*,108) 108 FORMAT(/,1X,'ENTER HIGHEST PERCENT DEFECTIVE SHIPPED BY SUPPLIER: 1 ',/) READ(*,*) QP 110 WRITE(*,109) 109 FORMAT(/,1X,'ENTER LARGEST LOT SIZE OF INTEREST:',/) READ(*,*) N C C THE PRIOR ESTIMATE OF SUPPLIER QUALITY (GIVEN BY QM,QO,QP) CAN C BE UPDATED USING RESULTS FROM ALL PREVIOUS SAMPLING. THE UPDATE IS C ACCOMPLISHED BY SPECIFYING UPDAN AND UPDAA, WHERE : C C UPDAN=NUMBER OF SAMPLES TO BE USED IN THE UPDATE C UPDAA=NUMBER OF FAILURES THAT HAVE BEEN OBSERVED C UPDAN=0.0D0 UPDAA=0.0D0 ACC=' ACCEPT ' REJ=' REJECT ' ALW=' ALWAYS' ZER=' 0 ' PRIN1=' N ' PRIN2=' N' PRIN3='A ' IFLAG=0 IF (PATH .EQ. 2.0D0) THEN U=1.0D0-(QP+4.0D0*QM+QO)/6.0D0 SIGSQ=((QP-QO)/6.0D0)**2 AL=U*U*(1.0D0-U)/SIGSQ-U+UPDAN-UPDAA BE=U*(1.0D0-U)**2/SIGSQ+U-1.0D0+UPDAA ELSE AL=A BE=B END IF S=(AL*S0+BE*S1)/(AL+BE) M=0 STO0=S AM=0.0D0 TEMP1=(AL*A0+BE*A1)/(AL+BE) IF (TEMP1 .GT. S) GO TO 10 PRIN1=ZER PRIN2=ALW PRIN3=ACC STO0=TEMP1 10 TEMP=(AL*R0+BE*R1)/(AL+BE) IF (STO0 .LE. TEMP) GO TO 20 PRIN1=ZER PRIN2=ALW PRIN3=REJ STO0=TEMP 20 STO(1,2)=N D0=STO0 1 M=M+1 AM=M SCOST=AM*S DCOST=N*STO0 IF (IFLAG .GT. 0) DCOST=ISTO(1,IFLAG+1)*S+(N-ISTO(1,IFLAG+1))* $ STO(2,IFLAG+1) IF (SCOST .GT. DCOST) GO TO 5 IF (AM .GT. N) GO TO 5 INT=((AM+AL)*(R0-A0)+BE*(R1-A1))/(R0+A1-R1-A0) IF (INT .LT. 0) INT=0 IF (INT .GT. M) INT=M LOOP=INT+1 D1=TEMP CONST=LGAMMA(AM+1.0D0)+LGAMMA(AL+BE)-LGAMMA(AL)-LGAMMA(BE) $ -LGAMMA(AM+AL+BE) DO 200 I=1,LOOP AJ=I-1 PROD=((A0-R0)*(AM+AL-AJ)+(A1-R1)*(BE+AJ))/(AM+AL+BE) IF (PROD .EQ. 0) GO TO 200 PLOG=DLOG(DABS(PROD)) ARG=CONST+LGAMMA(AM+AL-AJ)+LGAMMA(BE+AJ)-LGAMMA(AJ+1.0D0)- $ LGAMMA(AM+1.0D0-AJ) IF ((ARG+PLOG) .LT. -180.20D0) GO TO 200 C C CHECK TO PREVENT MACHINE UNDERFLOW. ERRORS ARE NOT CUMULATIVE C SO SET THAT TERM TO 0.0 WITHOUT SACRIFICING ACCURACY OF THE FINAL C RESULT. C D1=D1+PROD*DEXP(ARG) 200 CONTINUE IF (D1 .GE. D0) GO TO 1 2 IF (IFLAG .EQ. 0) GO TO 4 AN=(S*(AM-AM0)+AM0*D0-AM*D1)/(D0-D1) IF (AN .LT. AM) AN=AM IF (AN .GE. STO(1,IFLAG+1)) GO TO 3 IFLAG=IFLAG-1 D0=STO(2,IFLAG+1) AM0=ISTO(1,IFLAG+1) GO TO 2 3 IF (AN .GE. N) GO TO 1 IFLAG=IFLAG+1 IF (IFLAG .GE. 499) GO TO 6 C C CHECK TO SEE IF REFORMATTING IS NECESSARY. CURRENTLY, THE C PROGRAM CAN STORE ONLY 500 NODES, AT WHICH TIME BOTH COMPUTING AND C PRINTING WILL STOP. C STO(1,IFLAG+1)=AN STO(2,IFLAG+1)=D1 ISTO(1,IFLAG+1)=M ISTO(2,IFLAG+1)=INT D0=D1 AM0=AM GO TO 1 4 AN=AM*(S-D1)/(STO0-D1) IF (AN .LT. AM) AN=AM IF (AN .GT. N) GO TO 1 STO(1,2)=AN STO(2,2)=D1 ISTO(1,2)=M ISTO(2,2)=INT D0=D1 AM0=AM IFLAG=1 GO TO 1 5 IF (PATH .EQ. 2.0D0) THEN QP=100.0D0*QP QM=100.0D0*QM QO=100.0D0*QO WRITE(*,1000) S0,S1,A0,A1,R0,R1,QP,QM,QO,STO(1,2),PRIN1,PRIN2, $ PRIN3,STO0 ELSE WRITE(*,1010) STO(1,2),PRIN1,PRIN2,PRIN3,STO0 END IF 1010 FORMAT(10(/),' DETERMINATION OF OPTIMAL SAMPLING PLAN',///, & ' LOT SIZE N SAMPLE REJECT IF # EXPECTED',/ & ' SIZE FAILURES > COST =',// & ' 0.0 COST =',// $ ' 0.0 0, USING C ALGORITHM 291, COMMUNICATIONS OF THE A.C.M., VOLUME 9, NUMBER 9, C SEPTEMBER, 1966. C DOUBLE PRECISION FUNCTION LGAMMA(X) REAL*8 F,X,Y,Z Y=X IF (Y .LT. 7.0D0) THEN F=1.0D0 Z=Y-1.0D0 10 Z=Z+1.0D0 IF ( Z .LT. 7.0D0) THEN Y=Z F=F*Z GO TO 10 END IF Y=Y+1.0D0 F=-DLOG(F) ELSE F=0.0D0 END IF Z=1.0D0/Y**2 LGAMMA=F+(Y-.50)*DLOG(Y)-Y+.918938533204673 + $ (((-.000595238095238D0*Z + .000793650793651D0)*Z- $ .002777777777778D0)*Z+.083333333333333D0)/Y RETURN END