PROGRAM GAGERR C BY ANDY CHIANG date 1-9-98 Friday C DECLARE VARIABLES INTEGER D(4),CHOICE CHARACTER OK CHARACTER*20 OUTFIL,INFIL REAL S1,S2,S3,S4,M,I,J,SS(4) REAL DF(4),R2,R3,R4,P1,P3,P4 REAL RE,VRE,LRE,URE,SERE REAL OV,VOV,LOV,UOV,SEOV REAL PO,VPO,LPO,UPO,WUPO,WLPO,SEPO REAL PA,VPA,LPA,UPA,WUPA,WLPA,SEPA REAL RED,VRED,LRED,URED,SERED REAL OVD,VOVD,LOVD,UOVD,SEOVD REAL POD,VPOD,LPOD,UPOD,SEPOD REAL PAD,VPAD,LPAD,UPAD,SEPAD REAL U1,U12,U13,U14,VU1,U2,U22,U23,U24,VU2 REAL U3,U31,U32,U33,U34,VU3 REAL V1,V2,v3 REAL XX,GS,A,FU3,FU2 REAL LA,UA,LB,UB,LAB,UAB,LE,UE,CONFID,ZERO REAL UA30,UB30,UC30,D3,E3,LA30,LB30,LC30 REAL UA31,UB31,UC31,LA31,LB31,LC31 REAL UA1,UB1,UC1,D1,E1,LA1,LB1,LC1 REAL CL1,CL2,CL3,CU1,CU2,CU3,SW1,SW2 COMMON /DEG/ DF,A COMMON /RS/ R2,R3,R4 COMMON /MS/ S1,S2,S3,S4 COMMON /PS/ P1,P3,P4 ZERO=0. C GIVE OPTION 388 CONTINUE 1 PRINT 2 2 FORMAT(/////,79('-'),//,25X,'GAGE PROGRAM',//,79('-')) PRINT *,'1. ENTER INPUT VALUES FROM KEYBOARD:' PRINT *,'2. READ VALUES FROM EXTERNAL FILE:' 111 PRINT 94 94 FORMAT(/,'#SELECT 1, 2 OR TO QUIT ANYTIME:',/) READ *,CHOICE PRINT 99 C KEY-IN DATA IF(CHOICE.EQ.1) THEN 3 PRINT *,'#ENTER SUM OF SQUARES: SSA, SSB, SSAB, SSERROR: (>0)' READ *,(SS(K), K = 1,4) PRINT 97 SMALL=MIN(SS(1),SS(2),SS(3),SS(4)) IF(SMALL.LE.0.)THEN PRINT *,'NOT ALL SS>0, PLEASE RE-KEY' PRINT 97 GOTO 3 ENDIF 4 PRINT *,'#ENTER I,J,M ( INTEGERS > 1):' READ *,I,J,M PRINT 97 SMALL=MIN(I,J,M) IF(SMALL.LT.2)THEN PRINT *,'NOT ALL VALUES > 1, PLEASE RE-ENTER' PRINT 97 GOTO 4 ENDIF PRINT 99 C READ DATA FROM EXTERNAL FILE ELSEIF(CHOICE.EQ.2) THEN PRINT *,'ENSURE THAT DATA LISTED (IN FILE) AS:' PRINT *,'--------------------------------' PRINT *,' LINE 1: SSA' PRINT *,' LINE 2: SSB' PRINT *,' LINE 3: SSAB' PRINT *,' LINE 4: SSE' PRINT *,' LINE 5: I' PRINT *,' LINE 6: J' PRINT *,' LINE 7: M' PRINT *,'--------------------------------' PRINT 97 C 388 PRINT *,'#ENTER EXTERNAL FILE NAME:' READ(*,'(A)') INFIL PRINT 97 OPEN(2,FILE=INFIL,STATUS='OLD',ERR=899) READ(2,*) SS(1),SS(2),SS(3),SS(4),I,J,M CLOSE(2) ELSE GOTO 111 ENDIF C ENTER CONFIDENCE LEVEL 5 PRINT *,'#ENTER DESIRED CONFIDENCE LEVEL (50-99):' READ *,CONFID PRINT 97 IF ((CONFID.GE.99.9999).OR.(CONFID.LE.49.)) THEN PRINT *,'CONFID LEVEL SHOULD BE BETWEEN 50 & 99, PLEASE RE-KEY' PRINT 97 GOTO 5 ENDIF C C COMPUTE DEGREES OF FREEDOM 13 DF(1) = I - 1. DF(2) = J - 1. DF(3) = (I - 1.)*(J - 1.) DF(4) = (M - 1.)*I*J C COMPUTE MEAN SQUARES S1 = SS(1) / DF(1) S2 = SS(2) / DF(2) S3 = SS(3) / DF(3) S4 = SS(4) / DF(4) C ANOVA TABLE PRINT 98 PRINT 25,'SOURCE','DF','SUM OF SQUARES','MEAN SQUARES' 25 FORMAT(A14,T15,3(2X,A14)) PRINT 27,'A',DF(1),SS(1),S1 PRINT 27,'B',DF(2),SS(2),S2 PRINT 27,'AB',DF(3),SS(3),S3 PRINT 27,'ERROR',DF(4),SS(4),S4 27 FORMAT(A14,T15,(2X, F14.0),3(2X,F14.8)) PRINT 99 PRINT *,'# IS TABLE CORRECT? (Y/N)' READ (*,'(A)') OK IF ((OK.EQ.'N').OR.(OK.EQ.'n')) GOTO 1 PRINT 99 C OPEN OUTPUT FILE PRINT *,'#ENTER THE NAME OF OUTPUT FILE (MAX.8 CHARACTER):' READ(*,'(A)') OUTFIL PRINT 99 OPEN(1,FILE=OUTFIL,STATUS='UNKNOWN') C F-TEST FOR INTERACTION TERM FTAB=S3/S4 PVALUE=BETAI(DF(4)/2.,DF(3)/2.,(DF(4)/(DF(4)+DF(3)*FTAB))) IF(PVALUE.GT..10)THEN PRINT 705 PRINT 77,FTAB,PVALUE 77 FORMAT(1X,'F = MSAB/MSE = ',G12.5,3X,'HAS P-VALUE = ',G12.5) PRINT 700 PAUSE 'HIT RETURN-KEY TO CONTINUE' ENDIF C C COMPUTE SPECIAL COEFFICIENTS R2 = 1. / (M*I) R3 = (I - 1.) / (M*I) R4 = 1. / M P1 = 1. / (M*J) P3 = (J - 1.) / (M*J) P4 = R4 A = .5 - CONFID/200 LA = S1*(1 - 3*SQRT(2 / DF(1))) UA = S1*(1 + 3*SQRT(2 / DF(1))) LB = S2*(1 - 3*SQRT(2 / DF(2))) UB = S2*(1 + 3*SQRT(2 / DF(2))) LAB = S3*(1 - 3*SQRT(2 / DF(3))) UAB = S3*(1 + 3*SQRT(2 / DF(3))) LE = S4*(1 - 3*SQRT(2 / DF(4))) UE = S4*(1 + 3*SQRT(2 / DF(4))) C C************************* C COMPUTE POINT ESTIMATES C************************* C SIG - REPEABILITY SQUARES RE = S4 C SIG - REPRODUCIBILITY SQUARES PO = MAX(0.,R2*S2 + R3*S3 - R4*S4) C SIG - OVERAL SQUARES OV = R2*S2 + R3*S3 + (1-R4)*S4 C SIG - PARTS SQUARES PA = MAX(0.,P1*S1 + P3*S3 - P4*S4) C SIG - REPREATABILITY RED = SQRT(RE) C SIG - REPRODUCIBILITY POD = SQRT(PO) C SIG - OVERAL OVD = SQRT(OV) C SIG - PARTS PAD = SQRT(PA) C SIG - REPEA SQUARES / SIG - OVERAL SQUARES U1 = RE / OV V1 = SQRT(U1) C SIG - REPRO SQUARES / SIG - OVERAL SQUARES U2 = PO / OV V2 = SQRT(U2) C SIG - PARTS SQUARES / SIG - OVERAVLL SQUARES U3 = PA / OV V3 = SQRT(U3) C**************************** C COMPUTE CONFIDENCE LIMITS C**************************** C SIG - REPEABILITY SQUARES LRE = DF(4)*S4 / XX(A,DF(4)) URE = DF(4)*S4 / XX(1 - A,DF(4)) C SIG - REPREATABILITY LRED = SQRT(LRE) URED = SQRT(URE) C SIG - OVERAL SQUARES LOV = (R2*S2*G(2))**2+(R3*S3*G(3))**2 + ((1-R4)*S4*G(4))**2 UOV = (R2*S2*H(2))**2+(R3*S3*H(3))**2 + ((1-R4)*S4*H(4))**2 LOV = MAX(0., OV - SQRT(LOV)) UOV = OV + SQRT(UOV) C SIG - OVERAL LOVD = SQRT(LOV) UOVD = SQRT(UOV) C SIG - REPRODUCIBILITY SQUARES WUPO = (R2*S2*H(2))**2 + (R3*S3*H(3))**2 + (R4*S4*G(4))**2 WUPO = WUPO + R2*S2*R4*S4*HH(2,4) + R3*S3*R4*S4*HH(3,4) WLPO = (R2*S2*G(2))**2 + (R3*S3*G(3))**2 + (R4*S4*H(4))**2 WLPO = WLPO + R2*S2*R4*S4*GG(2,4) + R3*S3*R4*S4*GG(3,4) WLPO = WLPO + R2*S2*R3*S3*GS(2,3) LPO = MAX(0., PO - SQRT(WLPO)) UPO = PO + SQRT(WUPO) C SIG - REPRODUCIBILITY LPOD = SQRT(LPO) UPOD = SQRT(UPO) C SIG - PARTS SQUARES WUPA = (P1*S1*H(1))**2 + (P3*S3*H(3))**2 + (P4*S4*G(4))**2 WUPA = WUPA + P1*S1*P4*S4*HH(1,4) + P3*S3*P4*S4*HH(3,4) WLPA = (P1*S1*G(1))**2 + (P3*S3*G(3))**2 + (P4*S4*H(4))**2 WLPA = WLPA + P1*S1*P4*S4*GG(1,4) + P3*S3*P4*S4*GG(3,4) WLPA = WLPA + P1*S1*P3*S3*GS(1,3) LPA = MAX(0., PA - SQRT(WLPA)) UPA = PA + SQRT(WUPA) C SIG - PARTS LPAD = SQRT(LPA) UPAD = SQRT(UPA) C SIG - REPEA SQUARES / SIG - OVERAL SQUARES UA1 =(H(4)*S4*(1-R4))**2 + (G(2)*R2*S2)**2 + (G(3)*R3*S3)**2 + - HH(4,2)*S4*S2*(1-R4)*R2 - HH(4,3)*S4*S3*(1-R4)*R3 UB1 = -2*(H(4)*S4)**2*(1-R4) + HH(4,2)*S4*S2*R2 + + HH(4,3)*S4*S3*R3 UC1 = (H(4)*S4)**2 D1 = - R2*S2 - R3*S3 - (1-R4)*S4 E1 = S4 LA1 = (G(4)*S4*(1-R4))**2 + (H(2)*R2*S2)**2 + (H(3)*R3*S3)**2 + - GG(4,2)*S4*S2*(1-R4)*R2 - GG(4,3)*S4*S3*(1-R4)*R3 LB1 = -2*(G(4)*S4)**2*(1-R4) + GG(4,2)*S4*S2*R2 + + GG(4,3)*S4*S3*R3 LC1 = (G(4)*S4)**2 CU1 = UQUAD(UA1,UB1,UC1,D1,E1) CU1 = MIN(1.,CU1) CL1 = LQUAD(LA1,LB1,LC1,D1,E1) CL1 = MAX(0.,CL1) CL2 = 1 - CU1 CU2 = 1 - CL1 C SIG - PARTS SQUARES / SIG - OVERAL SQUARES UA30 = (G(2)*R2*S2)**2+(G(3)*R3*S3)**2 +(G(4)*(1-R4)*S4)**2 UB30 = - 2*(G(3)*S3)**2*R3*P3 + 2*(G(4)*S4)**2*(1-R4)*P4 + + HH(1,2)*P1*R2*S1*S2 + HH(1,3)*P1*R3*S1*S3 + + HH(1,4)*P1*(1-R4)*S1*S4 UC30 = (H(1)*P1*S1)**2 + (G(3)*P3*S3)**2 +(G(4)*P4*S4)**2 - + HH(1,3)*P1*P3*S1*S3 + HH(1,4)*P1*P4*S1*S4 D3 = D1 E30 = P1*S1 + P3*S3 - P4*S4 LA30 = (H(2)*R2*S2)**2 + (H(3)*R3*S3)**2 + (H(4)*(1-R4)*S4)**2 LB30 = - 2*(H(3)*S3)**2*R3*P3 + 2*(H(4)*S4)**2*(1-R4)*P4 + + GG(1,2)*P1*R2*S1*S2 + GG(1,3)*P1*R3*S1*S3 + + GG(1,4)*P1*(1-R4)*S1*S4 LC30 = (G(1)*P1*S1)**2 + (H(3)*P3*S3)**2 + (H(4)*P4*S4)**2 - + GG(1,3)*P1*P3*S1*S3 + GG(1,4)*P1*P4*S1*S4 CU3 = UQUAD(UA30,UB30,UC30,D3,E3) CL3 = LQUAD(LA30,LB30,LC30,D3,E3) CL3 = MAX(0.,CL3) IF(CU3.LT.(P3/R3)) THEN UA31 = (H(3)*S3*R3)**2 +(G(2)*S2*R2)**2 + (G(4)*(1-R4)*S4)**2 + - HH(3,2)*S3*S2*R3*R2 - HH(3,4)*S3*S4*R3*(1-R4) UB31 = -2*(H(3)*S3)**2 *P3*R3 + 2*(G(4)*S4)**2*P4*(1-R4) + + HH(1,2)*S1*S2*P1*R2 + HH(1,4)*S1*S4*P1*(1-R4) + + HH(3,2)*S3*S2*P3*R2 + HH(3,4)*S3*S4*(P3*(1-R4)-P4*R3) UC31 = (H(1)*P1*S1)**2 + (H(3)*P3*S3)**2 + (G(4)*P4*S4)**2 + + HH(1,4)*S1*S4*P1*P4 + HH(3,4)*S3*S4*P3*P4 LA31 = (G(3)*S3*R3)**2+(H(2)*S2*R2)**2 + (H(4)*(1-R4)*S4)**2 + - GG(3,2)*S3*S2*R3*R2 - GG(3,4)*S3*S4*R3*(1-R4) LB31 = -2*(G(3)*S3)**2 *P3*R3 + 2*(H(4)*S4)**2*P4*(1-R4) + + GG(1,2)*S1*S2*P1*R2 + GG(1,4)*S1*S4*P1*(1-R4) + + GG(3,2)*S3*S2*P3*R2 + GG(3,4)*S3*S4*(P3*(1-R4)-P4*R3) LC31 = (G(1)*P1*S1)**2 + (G(3)*P3*S3)**2 + (H(4)*P4*S4)**2+ + GG(1,4)*S1*S4*P1*P4 + GG(3,4)*S3*S4*P3*P4 CU3 = UQUAD(UA31,UB31,UC31,D3,E3) CL3 = LQUAD(LA31,LB31,LC31,D3,E3) CL3 = MAX(0.,CL3) ENDIF C********************************************* C COMPUTE STANDARD ERRORS OF POINT ESTIMATES C********************************************* C SIG-REPEABILITY SQUARES VRE = 2*S4**2 / DF(4) SERE = SQRT(VRE) C SIG-REPEABILITY VRED = VRE/(4*RE) SERED = SQRT(VRED) C SIG-OVERAL SQUARES VOV = 2*((R2*S2)**2 / DF(2) + (R3*S3)**2 / DF(3) + + ((1-R4)*S4)**2 / DF(4)) SEOV = SQRT(VOV) C SIG-OVERAL VOVD = VOV/(4*OV) SEOVD = SQRT(VOVD) C SIG-REPRODUCIBILITY SQUARES VPO = 2*((R2*S2)**2 / DF(2) + (R3*S3)**2 / DF(3) + + (R4*S4)**2 / DF(4)) SEPO = SQRT(VPO) C SIG-REPRODUCIBILITY IF (PO.GT.0.) THEN C BY TAYLOR APPROX VPOD = VPO/(4*PO) SEPOD = SQRT(VPOD) ELSE VPOD = 0. SEPOD = 0. PRINT 500 PAUSE 'HIT RETURN TO CONTINUE' ENDIF C SIG-PARTS SQUARES VPA = 2*((P1*S1)**2 / DF(1) + (P3*S3)**2 / DF(3) + + (P4*S4)**2 / DF(4)) SEPA = SQRT(VPA) C SIG-PARTS IF (PA.GT.0.)THEN C BY TAYLOR APPROX VPAD = VPA / (4*PA) SEPAD = SQRT(VPAD) ELSE VPAD = 0. SEPAD = 0. PRINT 600 PAUSE 'HIT RETURN TO CONTINUE' ENDIF C SIG-REPEA SQ / SIG-OVERAL SQ U12 = R2*S4 / OV**2 U13 = R3*S4 / OV**2 U14 = (R2*S2 + R3*S3) / OV**2 VU1 = 2*((U12*S2)**2 / DF(2) + (U13*S3)**2 / DF(3) + + (U14*S4)**2 / DF(4)) C SIG-REPEA / SIG-OVERAL VV1 = VU1/(4*U1) C SIG-REPRODUCIBILITY SQ / SIG-OVERAL SQ C BY TAYLOR APPROX U22 = R2*S4 / OV**2 U23 = R3*S4 / OV**2 U24 = (R4*RE - PO) / OV**2 VU2 = 2*((U22*S2)**2 / DF(2) + (U23*S3)**2 / DF(3) + + (U24*S4)**2 / DF(4)) IF(VU2.EQ.0.)THEN C BY AVERAGE SLOPE (BOX,HUNTER & HUNTER) SW1=1 PRINT 550 PAUSE 'HIT RETURN TO CONTINUE' U22 = (FU2(UB,S3,S4) - FU2(LB,S3,S4)) / (UB - LB) U23 = (FU2(S2,UAB,S4) - FU2(S2,LAB,S4)) / (UAB - LAB) U24 = (FU2(S2,S3,LE) - FU2(S2,S3,UE)) / (UE - LE) VU2 = 2*((U22*S2)**2 / DF(2) + (U23*S3)**2 / DF(3) + + (U24*S4)**2 / DF(4)) ENDIF C SIG-REPRODUCIBILITY / SIG-OVERAL C BY TAYLOR APPROX ONLY IF(PO.GT.0.)THEN VV2 = VU2/(4*U2) ELSE VV2 = 0. ENDIF C SIG-PARTS SQ / SIG-OVERAL SQ C BY TAYLOR APPROX U31 = P1/OV U32 = -R2*PA/OV**2 U33 = (P3*OV - R3*PA) / OV**2 U34 = ( - P4*OV - (1-R4)*PA) / OV**2 VU3 = 2*((U31*S1)**2 / DF(1) + (U32*S2)**2 / DF(2) + + (U33*S3)**2 / DF(3) + (U34*S4)**2 / DF(4)) IF(VU3.EQ.0.)THEN C BY AVERAGE SLOPE (BOX,HUNTER & HUNTER) SW2=1 PRINT 650 PAUSE 'HIT RETURN TO CONTINUE' U31 = (FU3(UA,S2,S3,S4) - FU3(LA,S2,S3,S4)) / (UA - LA) U32 = (FU3(S1,UB,S3,S4) - FU3(S1,LB,S3,S4)) / (UB - LB) U33 = (FU3(S1,S2,UAB,S4) - FU3(S1,S2,LAB,S4)) /(UAB - LAB) U34 = (FU3(S1,S2,S3,UE) - FU3(S1,S2,S3,LE)) / (UE - LE) VU3 = 2*((U31*S1)**2 / DF(1) + (U32*S2)**2 / DF(2) + + (U33*S3)**2 / DF(3) + (U34*S4)**2 / DF(4)) ENDIF C SIG-PARTS / SIG-OVERAL C BY TAYLOR APPROX ONLY IF(PA.GT.0.)THEN VV3 = VU3/(4*U3) ELSE VV3 = 0. ENDIF C C OUTPUT RESULT TO FILE C PRINT ANOVA TABLE C DO 33 I = 1,4 C D(I) = INT(DF(I)) C 33 CONTINUE PRINT 34,CONFID PRINT 35 PRINT 36 34 FORMAT(52X,F5.2,' % CONFIDENCE INTERVAL') 35 FORMAT(4X,' PARAMETER',9X,'ESTIMATE',7X,'STD ERROR', + 5X,'LOWER LIMIT',7X,'UPP LIMIT') 36 FORMAT(5X,9('-'),9X,8('-'), 7X,9('-'), 5X,11('-'), 7X,9('-')) PRINT 40,'SIG2_REPEA' ,RE,SERE,LRE,URE PRINT 40,'SIG2_REPRO' ,PO,SEPO,LPO,UPO PRINT 40,'SIG2_OVERAL' ,OV,SEOV,LOV,UOV PRINT 40,'SIG2_PARTS' ,PA,SEPA,LPA,UPA PRINT 97 PRINT 40,'SIG_REPEA' ,RED,SERED,LRED,URED PRINT 40,'SIG_REPRO' ,POD,SEPOD,LPOD,UPOD PRINT 40,'SIG_OVERAL' ,OVD,SEOVD,LOVD,UOVD PRINT 40,'SIG_PARTS' ,PAD,SEPAD,LPAD,UPAD 40 FORMAT(A14,1X,4(2X, F14.8)) PRINT 99 PAUSE 'HIT RETURN-KEY TO CONTINUE' PRINT 43,CONFID PRINT 44 PRINT 45 43 FORMAT( 55X,F5.2,' % CONFID INTERVAL') 44 FORMAT(8X,' PARAMETER',14X,'ESTIMATE',4X,'STD ERROR',2X, + 'LOWER LIMIT',4X,'UPP LIMIT') 45 FORMAT(9X,9('-'),14X,8('-'),4X,9('-'),2X,11('-'),4X,9('-')) PRINT 50,'SIG2_REPEA/SIG2_OVERAL' ,U1,SQRT(VU1),CL1,CU1 PRINT 50,'SIG2_REPRO/SIG2_OVERAL' ,U2,SQRT(VU2),CL2,CU2 PRINT 50,'SIG2_PARTS/SIG2_OVERAL' ,U3,SQRT(VU3),CL3,CU3 PRINT 97 PRINT 50,'SIG_REPEA/SIG_OVERAL',V1,SQRT(VV1),SQRT(CL1), + SQRT(CU1) PRINT 50,'SIG_REPRO/SIG_OVERAL',V2,SQRT(VV2),SQRT(CL2), + SQRT(CU2) PRINT 50,'SIG_PARTS/SIG_OVERAL',V3,SQRT(VV3),SQRT(CL3), + SQRT(CU3) 50 FORMAT(A26,1X,2(2X, F11.7),2(2X,F11.5)) PRINT 98 PAUSE 'HIT RETURN-KEY TO CONTINUE' WRITE(1,*)'FILE: ',OUTFIL WRITE(1,98) WRITE(1,51) 51 FORMAT(25X,'GAGE PROGRAM OUTPUT') WRITE(1,98) WRITE (1,25)'SOURCE','DEG. OF FREEDOM','SUM OF SQUARES', + 'MEAN SQUARES' WRITE(1,27) 'A',DF(1),SS(1),S1 WRITE(1,27) 'B',DF(2),SS(2),S2 WRITE(1,27) 'AB',DF(3),SS(3),S3 WRITE(1,27) 'ERROR',DF(4),SS(4),S4 WRITE(1,97) WRITE(1,34)CONFID WRITE(1,35) WRITE(1,36) WRITE(1,40)'SIG2_REPEA' ,RE,SERE,LRE,URE WRITE(1,40)'SIG2_REPRO' ,PO,SEPO,LPO,UPO WRITE(1,40)'SIG2_OVERAL' ,OV,SEOV,LOV,UOV WRITE(1,40)'SIG2_PARTS' ,PA,SEPA,LPA,UPA WRITE(1,97) WRITE(1,40)'SIG_REPEA' ,RED,SERED,LRED,URED WRITE(1,40)'SIG_REPRO' ,POD,SEPOD,LPOD,UPOD WRITE(1,40)'SIG_OVERAL' ,OVD,SEOVD,LOVD,UOVD WRITE(1,40)'SIG_PARTS' ,PAD,SEPAD,LPAD,UPAD WRITE(1,97) WRITE(1,43)CONFID WRITE(1,44) WRITE(1,45) WRITE(1,50)'SIG2_REPEA/SIG2_OVERAL' ,U1,SQRT(VU1),CL1,CU1 WRITE(1,50)'SIG2_REPRO/SIG2_OVERAL' ,U2,SQRT(VU2),CL2,CU2 WRITE(1,50)'SIG2_PARTS/SIG2_OVERAL' ,U3,SQRT(VU3),CL3,CU3 WRITE(1,97) WRITE(1,50)'SIG_REPEA/SIG_OVERAL',V1,SQRT(VV1),SQRT(CL1), + SQRT(CU1) WRITE(1,50)'SIG_REPRO/SIG_OVERAL',V2,SQRT(VV2),SQRT(CL2), + SQRT(CU2) WRITE(1,50)'SIG_PARTS/SIG_OVERAL',V3,SQRT(VV3),SQRT(CL3), + SQRT(CU3) WRITE(1,98) IF(PVALUE.GT..1) THEN WRITE(1,705) WRITE(1,77) FTAB, PVALUE WRITE(1,700) ENDIF IF(PA.LE.0.)WRITE(1,600) IF(PO.LE.0.)WRITE(1,500) IF(SW1.EQ.1)WRITE(1,550) IF(SW2.EQ.1)WRITE(1,650) CLOSE(1) PRINT 87 87 FORMAT(//,'#RUN PROGRAM AGAIN? [''N'' TO END]') READ (*,'(A)') OK IF((OK.EQ.'N').OR.(OK.EQ.'n'))GOTO99999 GOTO 1 97 FORMAT(/) 98 FORMAT(/,79('*'),/) 99 FORMAT(/,79('-'),/) 500 FORMAT(//, + //,' NOTE: SIG_REPRO HAS ZERO ESTIMATE AND SO', + /,12X,'STD ERROR OF [SIG_REPROD/SIG_OVERAL] AND ', + /,12X,'STD ERROR OF [SIG_REPROD] ARE NOT ESTIMATED.',/) 550 FORMAT(//, + //,7X,'* STD ERROR OF [SIG2_REPRO SQ/SIG2_OVERAL] IS EST BY', + /,11X,'THE BOX-HUNTER-HUNTER PROCEDURE IN WHICH EACH 1ST ORDER', + /,11X,'DERIVATIVE IN TAYLOR APPROX HAS BEEN REPLACED BY', + /,11X,'A 3-SIGMA AVERAGE SLOPE IN THE DIRECTION OF THE MS.',/) 600 FORMAT(//, + //,' NOTE: SIG_PARTS HAS ZERO ESTIMATE AND SO', + /,12X,'STD ERROR OF [SIG_PARTS/SIG_OVERAL] AND ', + /,12X,'STD ERROR OF [SIG_PARTS] ARE NOT ESTIMATED.',/) 650 FORMAT(//, + //,7X,'* STD ERROR OF [SIG2_PARTS SQ/SIG2_OVERAL] IS EST BY', + /,11X,'THE BOX-HUNTER-HUNTER PROCEDURE IN WHICH EACH 1ST ORDER', + /,11X,'DERIVATIVE IN TAYLOR APPROX HAS BEEN REPLACED BY', + /,11X,'A 3-SIGMA AVERAGE SLOPE IN THE DIRECTION OF THE MS.',/) 700 FORMAT(//,' WARNING: THE F-TEST SUGGESTS THAT THE AB INTERACTION', + /,12x,'IS NOT SIGNIFICANT AT .10 LEVEL.', + //,10x,'THIS PROGRAM PRODUCES ESTIMATES THAT ARE FOR A MODEL', + /,12x,'WITH SIGNIFICANT INTERACTION TERM AND SO MAY NOT BE', + /,12x,' APPROPRIATE FOR OTHER REDUCED MODELS.',//,3(/)) 705 FORMAT(/,65('-'),/) 899 PRINT 900 GOTO 388 900 FORMAT(//,25('='),/,' ERROR READING FILE!!! ',/,25('='),///) 99999 END C C C C***************************** C* FUNCTIONS & SUBROUTINES * C***************************** C C CHI - SQUARE PERCENTILE NUMERIC FUNC C COMPUTES 100Pth UPPER PERCENTAGE POINT OF CHI-SQ DIST WITH Y DF C XINV USES MODIFIED PROGRAMS FROM THE NUMERICAL RECIPES BOOK, SEE SECTION BELOW C ALTERNATIVELY, YOU MAY REPLACE XINV WITH CHIIN WHICH IS AN IMSL SUBROUTINE C IF AVAILABLE FUNCTION XX(P,Y) REAL P,Y,Q Q = 1. - P C XX = CHIIN(Q,Y) XX = XINV(Y,P) RETURN END C C F - DIST PERCENTILE NUMERIC FUNCT C COMPUTES\ 100Pth UPPER PERCENTAGE POINT OF F DIST WITH (X,Y) DF C FINV USES MODIFIED PROGRAMS FROM THE NUMERICAL RECIPES BOOK, SEE SECTION BELOW C ALTERNATIVELY, YOU MAY REPLACE FINV WITH FIN WHICH IS AN IMSL SUBROUTINE C IF AVAILABLE C P = UPPER TAIL AREA C X = 1ST DEG OF FREEDOM C Y = 2ND DEG OF FREEDOM FUNCTION F(P,X,Y) REAL P,X,Y,Q Q = 1. - P C F = FIN(Q,X,Y) F = FINV(X,Y,P) RETURN END C G - FUNCTION INDEXED FUNCTION C I IS THE FACTOR INDEX FUNCTION G(I) REAL DF(4),P INTEGER I COMMON /DEG/ DF,P G = 1 - DF(I) / XX(P,DF(I)) RETURN END C GG - FUNCTION INDEXED FUNCTION FUNCTION GG(I,J) REAL DF(4),P INTEGER I,J COMMON /DEG/ DF,P GG =(F(P,DF(I),DF(J))-1)**2-(G(I)*F(P,DF(I),DF(J)))**2-H(J)**2 GG =GG / F(P,DF(I),DF(J)) RETURN END C H - FUNCTION INDEXED FUNCTION C I IS THE FACTOR INDEX FUNCTION H(I) REAL DF(4),P INTEGER I COMMON /DEG/ DF,P H = DF(I) / XX(1 - P,DF(I)) - 1 RETURN END C HH - FUNCTION INDEXED FUNCTION FUNCTION HH(I,J) REAL DF(4),P INTEGER I,J COMMON /DEG/ DF,P HH = (1. - F(1 - P,DF(I),DF(J)))**2 - + (H(I)*F(1 - P,DF(I),DF(J)))**2 - G(J)**2 HH = HH / F(1 - P,DF(I),DF(J)) RETURN END C GS - FUNCTION INDEXED FUNCTION FUNCTION GS(I,J) REAL DF(4),DFZ INTEGER I,J COMMON /DEG/ DF,P DFZ = DF(I) + DF(J) GS = ((1. - DFZ / XX(P,DFZ))**2)*DFZ*DFZ / (DF(I)*DF(J)) GS = GS - (G(I)**2)*DF(I) / DF(J) - (G(J)**2)*DF(J) / DF(I) RETURN END C FUNCTIONS FOR BOX-HUNTER-HUNTER AVERAGE SLOPE PROCEDURE C FU2,FU3,UQUAD,LQUAD C FU2 FUNCTION FUNCTION FU2(B,AB,E) REAL R2,R3,R4 COMMON /RS/ R2,R3,R4 FU2 = (R2*B + R3*AB - R4*E) / (R2*B + R3*AB + (1-R4)*E) RETURN END C FU3 FUNCTION FUNCTION FU3(A,B,AB,E) REAL R2,R3,R4,P1,P3,P4 COMMON /RS/ R2,R3,R4 COMMON /PS/ P1,P3,P4 FU3 = (P1*A + P3*AB - P4*E) / (R2*B + R3*AB + (1-R4)*E) RETURN END C UQUAD FUNCTION FUNCTION UQUAD(A,B,C,D,E) REAL A,B,C,D,E,AA,BB,CC,Q AA = D*D - A BB = 2.*D*E - B CC = E*E - C Q = BB*BB - 4.*AA*CC Q = MAX(0.,Q) UQUAD = ( - BB + SQRT(Q)) / (2*AA) RETURN END C LQUAD FUNCTION FUNCTION LQUAD(A,B,C,D,E) REAL A,B,C,D,E,AA,BB,CC,Q AA = D*D - A BB = 2.*D*E - B CC = E*E - C Q = BB*BB - 4.*AA*CC Q = MAX(0.,Q) LQUAD = ( - BB - SQRT(Q)) / (2*AA) RETURN END C FUNCTIONS FOR CHI-SQUARES/F DIST PERCENTAGE POINTS C **BETAI, BETACF, BETA, GAMMLIN, GAMMQ, GCF, GSER C ARE COPIED FROM CHAPTER 6, Numerical Recipes in Fortran, C Cambridge University Press, ISBN 0 521 43064-X C THESE 7 FUNCTIONS ARE USED TO COMPUTE THE CHI-SQ & F C UPPER TAIL PROB GIVEN A PERCENTILE VALUE C **XINV, FINV ARE MODIFIED FROM RTBIS C COPIED FROM CHAPTER 9, Numerical Recipes in Fortran, C Cambridge University Press, ISBN 0 521 43064-X C THESE 2 FUNCTIONS ARE USED TO SOLVE FOR INTEGRAL EQUATION C BY BISECTION METHOD (AVOIDING THE NEWTON-RAPHSON METHOD C BECAUSE OF DERIVATIVE PROBLEM) FUNCTION BETAI(A,B,X) REAL BETAI,A,B,X REAL BT,BETACF,GAMMLN IF(X.LT.0..OR.X.GT.1.)PAUSE 'BAD ARGUMENT X IN BETAI' IF(X.EQ.0..OR.X.EQ.1.)THEN BT=0. ELSE BT=EXP(GAMMLN(A+B)-GAMMLN(A)-GAMMLN(B)+A*LOG(X)+B*LOG(1.-X)) ENDIF IF(X.LT.(A+1.)/(A+B+2.))THEN BETAI=BT*BETACF(A,B,X)/A RETURN ELSE BETAI=1.-BT*BETACF(B,A,1.-X)/B RETURN ENDIF END FUNCTION BETA(Z,W) REAL BETA,W,Z REAL GAMMLN BETA=EXP(GAMMLN(Z)+GAMMLN(W)-GAMMLN(Z+W)) RETURN END C FUNCTION BETACF(A,B,X) INTEGER MAXIT REAL BETACF,A,B,X,EPS,FPMIN PARAMETER (MAXIT=1000,EPS=3.E-7,FPMIN=1.E-30) INTEGER M,M2 REAL AA,C,D,DEL,H,QAB,QAM,QAP QAB=A+B QAP=A+1. QAM=A-1. C=1. D=1.-QAB*X/QAP IF(ABS(D).LT.FPMIN)D=FPMIN D=1./D H=D DO 11 M=1,MAXIT M2=2*M AA=M*(B-M)*X/((QAM+M2)*(A+M2)) D=1.+AA*D IF(ABS(D).LT.FPMIN)D=FPMIN C=1.+AA/C IF(ABS(C).LT.FPMIN)C=FPMIN D=1./D H=H*D*C AA=-(A+M)*(QAB+M)*X/((A+M2)*(QAP+M2)) D=1.+AA*D IF(ABS(D).LT.FPMIN)D=FPMIN C=1.+AA/C IF(ABS(C).LT.FPMIN)C=FPMIN D=1./D DEL=D*C H=H*DEL IF(ABS(DEL-1.).LT.EPS)GOTO 1 11 CONTINUE PAUSE 'A OR B TOO BIG, OR MAXIT TOO SMALL IN BETACF' 1 BETACF=H RETURN END C FUNCTION FPROB(X,V1,V2,A) REAL X,V1,V2,A FPROB=BETAI(V2/2,V1/2,(V2/(V2+V1*X))) - A RETURN END C FUNCTION FINV(V1,V2,A) INTEGER JMAX REAL FINV,X1,X2,XACC,V1,V2,A PARAMETER (JMAX=1000,X1=1.E-20,X2=1.E9,XACC=1.E-5) EXTERNAL FPROB INTEGER J REAL DX,F,FMID,XMID FMID=FPROB(X2,V1,V2,A) F=FPROB(X1,V1,V2,A) IF(F*FMID.GE.0.) PAUSE 'ROOT MUST BE BRACKETED IN FINV' IF(F.LT.0.)THEN FINV=X1 DX=X2-X1 ELSE FINV=X2 DX=X1-X2 ENDIF DO 11 J=1,JMAX DX=DX*.5 XMID=FINV+DX FMID=FPROB(XMID,V1,V2,A) IF(FMID.LE.0.)FINV=XMID IF(ABS(DX).LT.XACC .OR. FMID.EQ.0.) RETURN 11 CONTINUE PAUSE 'TOO MANY BISECTIONS IN FINV' END C FUNCTION GAMMLN(XX) REAL GAMMLN,XX INTEGER J DOUBLE PRECISION SER,STP,TMP,X,Y,COF(6) SAVE COF,STP DATA COF,STP/76.18009172947146D0,-86.50532032941677D0, *24.01409824083091D0,-1.231739572450155D0,.1208650973866179D-2, *-.5395239384953D-5,2.5066282746310005D0/ X=XX Y=X TMP=X+5.5D0 TMP=(X+0.5D0)*LOG(TMP)-TMP SER=1.000000000190015D0 DO 11 J=1,6 Y=Y+1.D0 SER=SER+COF(J)/Y 11 CONTINUE GAMMLN=TMP+LOG(STP*SER/X) RETURN END C FUNCTION GAMMQ(A,X) REAL A,GAMMQ,X REAL GAMMCF,GAS4R,GLN IF(X.LT.0..OR.A.LE.0.)PAUSE 'BAD ARGUMENTS IN GAMMQ' IF(X.LT.A+1.)THEN CALL GSER(GAS4R,A,X,GLN) GAMMQ=1.-GAS4R ELSE CALL GCF(GAMMCF,A,X,GLN) GAMMQ=GAMMCF ENDIF RETURN END C SUBROUTINE GCF(GAMMCF,A,X,GLN) INTEGER ITMAX REAL A,GAMMCF,GLN,X,EPS,FPMIN PARAMETER (ITMAX=1000,EPS=3.E-7,FPMIN=1.E-30) INTEGER I REAL AN,B,C,D,DEL,H,GAMMLN GLN=GAMMLN(A) B=X+1.-A C=1./FPMIN D=1./B H=D DO 11 I=1,ITMAX AN=-I*(I-A) B=B+2. D=AN*D+B IF(ABS(D).LT.FPMIN)D=FPMIN C=B+AN/C IF(ABS(C).LT.FPMIN)C=FPMIN D=1./D DEL=D*C H=H*DEL IF(ABS(DEL-1.).LT.EPS)GOTO 1 11 CONTINUE PAUSE 'A TOO LARGE, ITMAX TOO SMALL IN GCF' 1 GAMMCF=EXP(-X+A*LOG(X)-GLN)*H RETURN END C SUBROUTINE GSER(GAS4R,A,X,GLN) INTEGER ITMAX REAL A,GAS4R,GLN,X,EPS PARAMETER (ITMAX=1000,EPS=3.E-7) INTEGER N REAL AP,DEL,SUM,GAMMLN GLN=GAMMLN(A) IF(X.LE.0.)THEN IF(X.LT.0.)PAUSE 'X < 0 IN GSER' GAS4R=0. RETURN ENDIF AP=A SUM=1./A DEL=SUM DO 11 N=1,ITMAX AP=AP+1. DEL=DEL*X/AP SUM=SUM+DEL IF(ABS(DEL).LT.ABS(SUM)*EPS)GOTO 1 11 CONTINUE PAUSE 'A TOO LARGE, ITMAX TOO SMALL IN GSER' 1 GAS4R=SUM*EXP(-X+A*LOG(X)-GLN) RETURN END C FUNCTION XPROB(X,V,A) REAL A,V,X EXTERNAL GAMMQ XPROB=GAMMQ(V/2,X/2) - A RETURN END C FUNCTION XINV(V,A) INTEGER JMAX REAL XINV,X1,X2,XACC,V,A PARAMETER (JMAX=1000,X1=1.E-20,X2=1.E9,XACC=1.E-5) EXTERNAL XPROB INTEGER J REAL DX,F,FMID,XMID FMID=XPROB(X2,V,A) F=XPROB(X1,V,A) IF(F*FMID.GE.0.) PAUSE 'ROOT MUST BE BRACKETED IN XINV' IF(F.LT.0.)THEN XINV=X1 DX=X2-X1 ELSE XINV=X2 DX=X1-X2 ENDIF DO 11 J=1,JMAX DX=DX*.5 XMID=XINV+DX FMID=XPROB(XMID,V,A) IF(FMID.LE.0.)XINV=XMID IF(ABS(DX).LT.XACC .OR. FMID.EQ.0.) RETURN 11 CONTINUE PAUSE 'TOO MANY BISECTIONS IN XINV' END