C************************************************************* C **** TAXFIT - FOR LS FORTRAN 15-NOV-1996 **** C **** COPYRIGHT Michael Birnbaum, all rights reserved C **** MODIFIED 3-MAY-1998 **** C *** No warranties or guarantees. C *** This program may be used for noncommercial, scholarly C *** purposes without charge. Publications should cite following: C C Birnbaum, M. H., & Chavez, A. (1997). Tests of Theories of C Decision Making: Violations of Branch Independence C and Distribution Independence. Organizational Behavior C and Human Decision Processes, 71(2), 161-194. C Birnbaum, M. H., & Stegner, S. E. (1979). Source credibility C in social judgment: Bias, expertise, and the judge's C point of view. Journal of Personality and Social Psychology, C 37, 48-74. C C*** Note: This version includes Chandler's STEPIT Subroutine. C*** If you use STEPIT, also cite Chandler (1969). C************************************************************* C C THIS PROGRAM WILL FIT SIMPLIFIED BIRNBAUM-STEGNER MODEL TO EXPT C WITH Judged Strength of pref TO GET ONE GAMBLE OR ANOTHER C FOR CHOICES BETWEEN TWO, THREE, OR FOUR OUTCOME GAMBLES C MODEL ALLOWS ONE CONFIGURAL WEIGHT PARAMETER, HERE CALLED C C THIS VERSION FITS SINGLE S DATA AS A FUNCTION OF UTILITY DIFF C THERE ARE TWO INDICES OF FIT: C 1. -LOG LIKELIHOOD OF CHOICES, GIVEN MODEL C 2. SUM OF SQUARED DEVIATIONS between Data and Model C DATA FORMAT (Each CARD is a line of 80 columns ) C CARD 1 2I3, F5.0 # OF TRIALS TO BE FIT, MODEL ID, 1-h C 1 - h is the weight of -Log(likelihood| Model) C h is the weight of sum of squared deviations C MODEL ID = 1 JUDGMENTS, MODEL ID = 2 is not active here C 2 I3 No. of Subjects (I3) C 3 5F5.2 NEXT CARDS (5F5.2) INITIAL ESTIMATES, C C Parameters are GAMMA, BETA, A, B, and DELTA (here called C) C Warning: Your design may not allow separation of parameters C In that case, you should fix one or more parameters C A good parameter to fix is beta = 1, which has been C found as a reasonable approximation for small outcomes. C You should also constrain delta to avoid negative weights C UNLIKE Birnbaum & Chavez, B is the scaling parameter for choice C A is the scaling parameter for judgment C C 4-6 5F5.2 MASKS, MIN, MAX-These constrain parameter values C EACH SET STARTS ON A NEW CARD-Format is 5F5.2 C MASK=0 is free, MASK=1 is fixed to initial estimate C XMIN =0 XMAX =0 for free parameters, C 7 FORMAT for mean judgments in this file C 8 FORMAT for raw data of subjects in file = SUBJECTS.DAT C 9 FORMAT for output of predictions and residuals C 10 2I3 # OF OUTCOMES LEFT, # OF OUTCOMES RIGHT C 11 (8F5.0) PROBABILITIES OF LEFT AND RIGHT OUTCOMES C 12 (8F5.0) OUTCOMES OF LEFT AND RIGHT GAMBLES C NOTE: lowest outcomes are on the left, highest on right C 13 (1F7.0) DATA VALUE (E.G., MEAN JUDGMENT OF STR preference) C C CARDS 10-13 REPEATED FOR # OF PROBLEMS C Note: This file MUST have the problems in the same order C as the order of judgments in the subjects data array. C Data coded: negative numbers = chose gamble on the left C positive numbers = chose gamble on the right C C ************************************************************** C PROGRAM TAXFIT IMPLICIT DOUBLEPRECISION(A-H,O-Z) COMMON/STPIT/CHISQ,X(50),XMAX(50),XMIN(50),DELTAX(50),ERR(50,50), +DELMIN(50),MASK(50),NV,NTRACE,MATRIX COMMON/TOFUNK/D(120,300),PD(120,300),RESD(120,300), 1GAMMA(120),BETA(120),A(120),B(120),NPAIRS,WWLIK, 2PX(300,5),PY(300,5),XO(300,5),YO(300,5),NOL(300),NOR(300), 3W(5),C(110),CWT(5,5),UTL(300),UTR(300),MID,ISUB,CHISS,CHILK COMMON/FASTER/JVARY COMMON/ICALL/ICALLS,IMAX DIMENSION FMTDAT(20),FMT2(20),FMT3(20) DIMENSION XS(50),D1(300) C C *** FILES TO BE USED ********************************** FILES C OPEN (5, FILE = 'TAXFIT.DAT', STATUS = 'OLD') OPEN (6, FILE = 'TAXFIT.OUT', STATUS = 'NEW') OPEN (7, FILE = 'SUBJECTS.DAT', STATUS = 'OLD') OPEN (8, FILE = 'TAXFIT.PAR', STATUS = 'NEW') OPEN (9, FILE = 'TAXFIT.RES', STATUS = 'NEW') OPEN (10, FILE = 'TAXFIT.PRE', STATUS = 'NEW') C C READ IN THE NUMBER OF TRIALS, GAMBLE PAIRS = NPAIRS READ(5,99)NPAIRS,MID,WWLIK READ(5,99)NSUBS NSP1 = NSUBS +1 99 FORMAT(2I3,F5.0) C NV IS THE NO.OF ESTIMATED PARAMETERS NV=5 C NTRACE=1 MEANS PRINT OUT OF MINIMIZATION C TO MAKE STEPIT HAPPY: NTRACE=0 MATRIX=110 C READ IN INITIAL ESTIMATES OF THE X VALUES C TRANSFER TABLE C GAMMA, BETA, A, B, C = DELTA C X(1) X(2) X(3) X(4) X(5) READ(5,151)(XS(J),J=1,NV) 151 FORMAT(10F5.2) WRITE(6,155)(XS(J),J=1,NV) 155 FORMAT(/,' INITIAL ESTIMATES',8(/,' ',7F9.2)) C SOME PARAMETERS ARE FIXED OR CONSTRAINED C TO FIX A PARAMETER, SET MASKS TO 1 --TO CONSTRAIN USE XMIN AND XMAX C MASK = 0 FOR FREE PARAMETERS XMIN = XMAX = 0 = UNCONSTRAINED C MASK = 1 FOR FIXED PARAMETERS C IF YOU WANT EVERYTHING FREE, LEAVE ALL ZEROS READ(5,161) (MASK(J),J=1,NV) 161 FORMAT(5I5) READ(5,151) (XMIN(J),J=1,NV) READ(5,151) (XMAX(J),J=1,NV) READ(5,167) FMTDAT READ(5,167) FMT2 READ(5,167) FMT3 167 FORMAT(20A4) C *************************** START DATA LOOP HERE ********* DO 110 IP=1,NPAIRS READ(5,99)NOL(IP),NOR(IP) C NOL = NO OF OUTCOMES LEFT, NOR = NO OF OUTCOMES ON RIGHT C READ IN PROBABILITIES NL = NOL(IP) NR = NOR(IP) READ(5,101)(PX(IP,J),J=1,NL),(PY(IP,K),K=1,NR) 101 FORMAT(8F5.0) READ(5,101)(XO(IP,J),J=1,NL),(YO(IP,K),K=1,NR) READ(5,FMTDAT) D1(IP) C READ IN THE DATA AVERAGED DATA, D1(IP), USING SPECIFIED FORMAT 102 FORMAT(F7.0) 110 CONTINUE TOTSS=0.0 TOTLK=0.0 C ******************************+++++++++++++++************ SUBJECT LOOP C START SUBJECT LOOP HERE C ********************************************************* DO 990 ISUB=1,NSUBS C READ SUBJECTS DATA READ(7,FMT2) (D(ISUB,IP),IP=1,NPAIRS),ISID C SET INITIAL VALUES FOR EACH SUBJECT DO 120 I =1,NV X(I) = XS(I) 120 CONTINUE C ICALLS IS NO. OF CALLS TO FUNK, IMAX IS THE UPPER LIMIT FOR ICALLS ICALLS=0 IMAX=2500 C SET VALUES OF DELMIN AND DELTAX DO 140I=1,NV DELMIN(I)=0.0 DELTAX(I)=X(I)*(.01) 140 CONTINUE C ********************************************************** WRITE(6,141) ISUB 141 FORMAT(' SUBJECT=',I5) C NOW CALL STEPIT CALL STEPIT C ********************************************************** C PRINT OUT RESULTS C PRINTOUTS PER SUBJECT C FILE OF PARAMETERS AND INDICES OF FIT WRITE(8,703) GAMMA(ISUB),BETA(ISUB),A(ISUB),B(ISUB),C(ISUB), 1 CHISS,CHILK,ISID 703 FORMAT(5F8.3,2F12.3,2X,I5) C FIND WEIGHTS FOR EACH SUBJECT OF 4 EQUALLY LIKELY OUTCOMES C CALCULATE W(P) TTT = 4*(.25)**GAMMA(ISUB) V1=(.25)**GAMMA(ISUB) V2=C(ISUB) W1 = V1*(1.0-(.6)*V2)/TTT W2 = (V1-(.2)*C(ISUB)*V1)/TTT W3 = (V1+(.2)*C(ISUB)*V1)/TTT W4 = (V1+(.6)*C(ISUB)*V1)/TTT C PRINT WEIGHTS FOR EACH SUBJECT WRITE(8,723) W1,W2,W3,W4,ISID 723 FORMAT(4F8.3,34X,I5) C CUMULATE TOTAL INDICES OF FIT TOTSS=TOTSS+CHISS TOTLK=TOTLK+CHILK C PRINT FILE OF RESIDUALS WRITE(9,FMT3) (RESD(ISUB,IP),IP=1,NPAIRS),ISID WRITE(10,FMT3) (PD(ISUB,IP),IP=1,NPAIRS),ISID C ****************************************************** 990 CONTINUE C C ********************************* SUBJECT LOOP ******* C END PRINTOUT PER SUBJECT HERE ************************* SUBJECT LOOP ENDS C WRITE(6,803) NSUBS,TOTSS,TOTLK 803 FORMAT(' FIT OF BIRNBAUM-STEGNER MODEL FOR',I5,2X,' SUBJECTS', 1 /, ' TOTAL SUM OF SQUARED DEVIATIONS =',F20.2, 2 /, ' TOTAL NEGATIVE LOG LIKELIHOOD =',F20.2) C FIND AVERAGE DEVIATION AVSS = TOTSS/(NSUBS*NPAIRS) AVLK = TOTLK/(NSUBS*NPAIRS) RMS = AVSS**.5 WRITE (6,825) AVSS,RMS,AVLK 825 FORMAT(/,10X,' AVERAGE SQUARED DEVIATION=',F10.3, 1 /,10X, ' ROOT MEAN SQUARED DEVIATION=',F10.3, 2 /,10X, ' AVERAGE NEG. LOG LIKELIHOOD=',F10.3) C C ************************************************************* C NOW FIT THE AVERAGED DATA --These results will appear as the last Subject NSP1 = NSUBS + 1 ISUB = NSP1 DO 901 IP = 1,NPAIRS D(ISUB,IP)=D1(IP) 901 CONTINUE C**************************----------- C SET INITIAL VALUES DO 920 I =1,NV X(I) = XS(I) 920 CONTINUE C ICALLS IS NO. OF CALLS TO FUNK, IMAX IS THE UPPER LIMIT FOR ICALLS ICALLS=0 IMAX=2500 C SET VALUES OF DELMIN AND DELTAX DO 940I=1,NV DELMIN(I)=0.0 DELTAX(I)=X(I)*(.01) 940 CONTINUE C ********************************************************** WRITE(6,941) 941 FORMAT(' FIT FOR AVERAGED DATA') C NOW CALL STEPIT CALL STEPIT C*************************---------------- WRITE(6,801) 801 FORMAT('1',T36,'FIT OF BIRNBAUM-STEGNER TAX MODEL-ONE OMEGA = C') C PRINT OUT GAMBLES FIT WRITE(6,804) 804 FORMAT(T15,/,T15, ' GAMBLE ON THE LEFT GAMBLE ON THE RIGHT') C----------BEGIN LOOP OF PAIRS HERE-------------- DO 810 IP=1,NPAIRS NL = NOL(IP) NR = NOR(IP) WRITE(6,821) IP,NOL(IP),NOR(IP) 821 FORMAT(/,T5,' TRIAL=',I4,' NO OUTS LEFT=',I4,8X, + ' NO OUTS RIGHT=',I4) WRITE(6,805) (PX(IP,J),J=1,NL),(PY(IP,K),K=1,NR) 805 FORMAT(,T10, ' PROBS=',2X,8F8.2) WRITE(6,806) (XO(IP,J),J=1,NL),(YO(IP,K),K=1,NR) 806 FORMAT(,T10,' OUTS =',2X,8F8.2) WRITE(6,822) UTL(IP),UTR(IP) 822 FORMAT(,T15, ' CE EQUIV LEFT=', F7.2,5X,' CE EQUIV RIGHT =',F7.2) WRITE(6,807) D(NSP1,IP),PD(NSP1,IP),RESD(NSP1,IP) 807 FORMAT(T5,' DATA =',3X,F10.2,3X,' PRED =',3X,F10.2,3X, + ' RESID=',3X,F10.2) 808 FORMAT(/,T20,' PRED =',10X,F10.2) 809 FORMAT(/T20,' RESID=',10X,F10.2) 810 CONTINUE WRITE(8,703) GAMMA(ISUB),BETA(ISUB),A(ISUB),B(ISUB),C(ISUB), 1 CHISS,CHILK,NSP1 C FIND WEIGHTS FOR GROUP DATA OF 4 EQUALLY LIKELY OUTCOMES C CALCULATE W(P) TTT = 4*(.25)**GAMMA(ISUB) V1=(.25)**GAMMA(ISUB) V2=C(ISUB) W1 = V1*(1.0-(.6)*V2)/TTT W2 = (V1-(.2)*C(ISUB)*V1)/TTT W3 = (V1+(.2)*C(ISUB)*V1)/TTT W4 = (V1+(.6)*C(ISUB)*V1)/TTT WRITE(8,723) W1,W2,W3,W4,ISUB WRITE(9,FMT3) (RESD(ISUB,IP),IP=1,NPAIRS),ISUB WRITE(10,FMT3) (PD(ISUB,IP),IP=1,NPAIRS),ISUB C PRINT WEIGHTS FOR GROUP DATA INTO SUBJECT FILE SUBJECT C C C IF WE GET THIS MESSAGE, ALL IS WELL WRITE(6,999) 999 FORMAT(///,T15,'NORMAL END') STOP END C C ******************************************************************** C **** **** C **** HERE IS SUBROUTINE ' FUNK ' **** C **** Copyright Michael Birnbaum **** C ******************************************************************** C C HERE IS FUNK--WHICH FITS BIRNBAUM-STEGNER MODEL WITH ONE PARAMETER C TO JUDGMENTS of Strength of Preferences between pairs of Gambles C Negative Numbers mean preference for the Gamble on the Left C Positive Numbers mean preference for the Gamble on the Right C 0 = No preference. C SUBROUTINE FUNK C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT DOUBLEPRECISION(A-H,O-Z) COMMON/STPIT/CHISQ,X(50),XMAX(50),XMIN(50),DELTAX(50),ERR(50,50), +DELMIN(50),MASK(50),NV,NTRACE,MATRIX COMMON/TOFUNK/D(120,300),PD(120,300),RESD(120,300), 1GAMMA(120),BETA(120),A(120),B(120),NPAIRS,WWLIK, 2PX(300,5),PY(300,5),XO(300,5),YO(300,5),NOL(300),NOR(300), 3W(5),C(110),CWT(5,5),UTL(300),UTR(300),MID,ISUB,CHISS,CHILK COMMON/FASTER/JVARY COMMON/ICALL/ICALLS,IMAX C DECODE THE TRANSFER TABLE GAMMA(ISUB) = X(1) BETA(ISUB) = X(2) A(ISUB) = X(3) B(ISUB) = X(4) C(ISUB) = X(5) C FIND SUM OF SQUARED RESIDUALS, CHISS, AND NEGATIVE LOG LIKELIHOOD, CHILK CHISQ = 0.0 CHISS = 0.0 CHILK = 0.0 C ************************************************************ C--------------IN THIS SECTION, WE CALCULATE C THE PREDICTION FOR EACH DATA POINT, FIND THE RESIDUAL C SQUARE THE RESIDUALS, AND ADD THEM UP INTO CHISS C--------------WE ALSO CALCULATE HERE NEGATIVE LOG LIKELIHOOD C OF CHOICES GIVEN THE PARAMETERS AND MODEL C ************************************************************* DO 550 IP =1,NPAIRS C COMPUTE UTILITIES OF THE GAMBLES ON THE LEFT AND THE RIGHT NL = NOL(IP) NR = NOR(IP) SUML = 0.0 TTL = 0.0 C GAMBLE ON THE LEFT ***************************** LEFT C FIRST LOOP== FIND SUMS OF WEIGHTS AND WEIGHT-UTILITY PRODS DO 520, I = 1,NL W(I) = PX(IP,I)**GAMMA(ISUB) SUML = SUML + W(I)*XO(IP,I)**BETA(ISUB) TTL = TTL + W(I) 520 CONTINUE C SECOND SET OF LOOPS==FIND SUM OF CONFIGURAL TERMS SUMC = 0.0 DO 525 I = 2,NL JL = I - 1 DO 525 J =1,JL TN = NL+1 XDIF = (XO(IP,I)**BETA(ISUB)) - (XO(IP,J)**BETA(ISUB)) IF (C(ISUB) .LE. 0.0) THEN SUMC = SUMC + C(ISUB)*(XDIF/TN)*(PX(IP,I)**GAMMA(ISUB)) ELSE SUMC = SUMC + C(ISUB)*(XDIF/TN)*(PX(IP,J)**GAMMA(ISUB)) END IF 525 CONTINUE UTL(IP) = (SUML + SUMC)/TTL C C GAMBLE ON THE RIGHT ***************************** RIGHT C FIRST LOOP== FIND SUMS OF WEIGHTS AND WEIGHT-UTILITY PRODS SUMR = 0.0 TTR = 0.0 DO 530, I = 1,NR W(I) = PY(IP,I)**GAMMA(ISUB) SUMR = SUMR + W(I)*YO(IP,I)**BETA(ISUB) TTR = TTR + W(I) 530 CONTINUE C SECOND SET OF LOOPS==FIND SUM OF CONFIGURAL TERMS SUMC = 0.0 DO 535 I = 2,NR JL = I - 1 DO 535 J =1,JL TN = NR+1 YDIF = YO(IP,I)**BETA(ISUB) - YO(IP,J)**BETA(ISUB) IF (C(ISUB) .LE. 0.0) THEN SUMC = SUMC + C(ISUB)*(YDIF/TN)*(PY(IP,I)**GAMMA(ISUB)) ELSE SUMC = SUMC + C(ISUB)*(YDIF/TN)*(PY(IP,J)**GAMMA(ISUB)) END IF 535 CONTINUE UTR(IP) = (SUMR + SUMC)/TTR C ******************************************* CALCULATE PREDICTED RESPONSE PD(ISUB,IP) = A(ISUB)*(UTR(IP) - UTL(IP)) C C Here WE CALCULATE SUM OF SQUARED DEVIATIONS AND SUM OF -LOG(PROBS) RESD(ISUB,IP) = D(ISUB,IP) - PD(ISUB,IP) CHISS = CHISS + RESD(ISUB,IP)**2 C PPP = B(ISUB)*(UTR(IP) - UTL(IP)) PP = (1.0/(1.0+EXP(PPP))) IF (D(ISUB,IP) .LE. 0.0) THEN CHILK = CHILK - LOG(PP) ELSE PP = 1.0 - PP CHILK = CHILK - LOG(PP) END IF C 550 CONTINUE CHISQ = (WWLIK)*CHILK +(1.0-WWLIK)*CHISS ICALLS=ICALLS+1 RETURN END C C ****************************************************************** C **** **** C **** THIS IS SUBROUTINE STEPIT **** C **** **** C ****************************************************************** SUBROUTINE STEPIT C COPYRIGHT 1965 -- J. P. CHANDLER, PHYSICS DEPT., INDIANA UNIVERSITY. C C STEPIT 5.2 AUGUST 1967 C AVAILABLE FROM.... QUANTUM CHEMISTRY PROGRAM EXCHANGE C I.U. CHEMISTRY DEPT., BLOOMINGTON, INDIANA. C C MODIFIED FOR FORTRAN G OS/360 MVT C C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT DOUBLEPRECISION(A-H,O-Z) DIMENSION VEC(50),TRIAL(50),XSAVE(50),CHI(50),DX(50),SECOND(2,2) DIMENSION OLDVEC(50),SALVO(50),XOSC(50,15),CHIOSC(15),NFLAT(50) C COMMON/STPIT/CHISQ,X(50),XMAX(50),XMIN(50),DELTAX(50),ERR(50,50) +,DELMIN(50),MASK(50),NV,NTRACE,MATRIX COMMON/FASTER/JVARY COMMON/ICALL/ICALLS,IMAX NVMAX=50 MOSQUE=15 KW=61 RATIO=10.0 COLIN=0.99 NCOMP=5 ACK=2.0 SIGNIF=2.E8 HUGE=1.E37 C JVARY=0 40 IF(NV)290,290,50 50 NACTIV=0 DO 150I=1,NV IF(MASK(I))150,60,150 60 IF(SIGNIF*DABS(DELTAX(I))-DABS(X(I)))70,70,100 70 IF(X(I))90,80,90 80 DELTAX(I)=0.01 GOTO 100 90 DELTAX(I)=0.01*X(I) 100 IF(DELMIN(I))120,110,120 110 DELMIN(I)=DELTAX(I)/SIGNIF 120 IF(XMAX(I)-XMIN(I))130,130,140 130 XMAX(I)=HUGE XMIN(I)=-HUGE 140 NACTIV=NACTIV+1 VACD=DMIN1(XMAX(I),X(I)) X(I)=DMAX1(XMIN(I),VACD) 150 CONTINUE COMPAR=0.0 IF(NACTIV-1)160,190,180 160 DO 170J=1,NV 170 MASK(J)=0 GOTO 50 180 A=NACTIV SUB=2.0/(A-1.0) P=2.0*(1.0/DSQRT(A)/(1.0-0.5**SUB)-1.0) COMPAR=DMIN1(9.99D-1,DABS((1.0-(1.0-COLIN)**SUB)*(1.0+P*(1.- +COLIN)))) 190 IF(NTRACE)280,200,200 200 WRITE(6,210) 210 FORMAT( +'1ENTER SUBROUTINE STEPIT. COPYRIGHT 1965 J. P. CHANDLER, PHYSICS + DEPT., INDIANA UNIVERSITY.'//' INITIAL VALUES....'/) WRITE(6,220)(MASK(J),J=1,NV) 220 FORMAT(/' MASK = '10(I6,6X)/(4X10I12)) WRITE(6,230)(X(J),J=1,NV) 230 FORMAT(/' X = '10E12.4/(10X10E12.4)) WRITE(6,240)(XMAX(J),J=1,NV) 240 FORMAT(/' XMAX = '10E12.4/(10X10E12.4)) WRITE(6,250)(XMIN(J),J=1,NV) 250 FORMAT(/' XMIN = '10E12.4/(10X10E12.4)) WRITE(6,260)(DELTAX(J),J=1,NV) 260 FORMAT(/' DELTAX = '10E12.4/(10X10E12.4)) WRITE(6,270)(DELMIN(J),J=1,NV) 270 FORMAT(/' DELMIN = '10E12.4/(10X10E12.4)) 280 CALL FUNK NF=1 JOCK=1 290 IF(NTRACE)320,300,300 300 WRITE(6,310)NV,NACTIV,MATRIX,NCOMP,RATIO,ACK,COLIN,COMPAR,CHISQ 310 FORMAT(//' 'I3,' VARIABLES,'I3,' ACTIVE.'10X8HMATRIX =I4, +10X7H NCOMP=I2//' RATIO ='F5.1,10X5HACK =F5.1,10X7HCOLIN =F6.3, +10X8HCOMPAR =F6.3///' CHISQ ='E15.8) 320 IF(NV)2150,2150,330 330 IF(NTRACE)360,360,340 340 WRITE(6,350) 350 FORMAT(//60(1X1H*)//10X37HTRACE MAP OF THE MINIMIZATION PROCESS//) C 360 DO 370I=1,NV DX(I)=DELTAX(I) VEC(I)=0. DO 370J=1,20 370 ERR(I,J)=0. CHIOLD=CHISQ NOSC=0 C 380 NCIRC=0 NZIP=0 C C MAIN DO LOOP FOR CYCLING THROUGH THE VARIABLES. C FIRST TRIAL STEP WITH EACH VARIABLE IS SEPARATE. 390 NACK=0 DO 1350I=1,NV OLDVEC(I)=VEC(I) VEC(I)=0.0 TRIAL(I)=0.0 IF(MASK(I))400,410,400 400 VEC(I)=-0.0 NFLAT(I)=1 GOTO 1350 410 NACK=NACK+1 SAVE=X(I) IF(SIGNIF*DABS(DX(I))-DABS(X(I)))580,580,420 420 X(I)=SAVE+DX(I) JVARY=0 IF(JOCK)440,440,430 430 JOCK=0 JVARY=I 440 NFLAG=1 IF(X(I)-XMIN(I))460,460,450 450 IF(X(I)-XMAX(I))470,460,460 460 NFLAG=NFLAG+3 GOTO 490 470 CALL FUNK IF(ICALLS.LE.IMAX)GOTO10 RETURN 10 CONTINUE NF=NF+1 JVARY=I CHIME=CHISQ IF(CHISQ-CHIOLD)620,480,490 480 NFLAG=NFLAG+1 490 X(I)=SAVE-DX(I) IF(X(I)-XMIN(I))590,590,500 500 IF(X(I)-XMAX(I))510,590,590 510 CALL FUNK NF=NF+1 JVARY=I IF(CHISQ-CHIOLD)610,520,530 520 NFLAG=NFLAG+1 530 IF(NFLAG-3)540,580,590 540 IF((CHISQ-CHIME)*(CHIME-2.*CHIOLD+CHISQ))550,590,550 550 TRIAL(I)=.5*DX(I)*(CHISQ-CHIME)/(CHIME-2.*CHIOLD+CHISQ) VEC(I)=TRIAL(I)/DABS(DX(I)) NFLAT(I)=0 X(I)=SAVE+TRIAL(I) CALL FUNK NF=NF+1 IF(CHISQ-CHIOLD)560,570,570 560 CHIOLD=CHISQ JOCK=1 GOTO 600 570 TRIAL(I)=0.0 VEC(I)=0.0 GOTO 590 580 VEC(I)=-0.0 NFLAT(I)=1 590 X(I)=SAVE 600 NCIRC=NCIRC+1 IF(NCIRC-NACTIV)690,1430,1430 610 DX(I)=-DX(I) C C A LOWER VALUE HAS BEEN FOUND. HENCE THIS VARIABLE WILL CHANGE. C 620 NCIRC=0 DEL=DX(I) 630 CHIME=CHIOLD CHIOLD=CHISQ VEC(I)=VEC(I)+DEL/DABS(DX(I)) NFLAT(I)=0 TRIAL(I)=TRIAL(I)+DEL DEL=ACK*DEL SAVE=X(I) X(I)=SAVE+DEL IF(X(I)-XMIN(I))680,680,640 640 IF(X(I)-XMAX(I))650,680,680 650 CALL FUNK NF=NF+1 IF(CHISQ-CHIOLD)630,660,660 660 ZZZ=ACK*CHIME-(ACK+1.0)*CHIOLD+CHISQ IF(ZZZ.EQ.0.)GOTO661 CINDER=(0.5/ACK)*(ACK**2*CHIME-(ACK**2-1.0)*CHIOLD-CHISQ)/ZZZ GOTO 662 661 CINDER=0.0 662 X(I)=SAVE+CINDER*DEL CALL FUNK NF=NF+1 IF(CHISQ-CHIOLD)670,680,680 670 CHIOLD=CHISQ TRIAL(I)=TRIAL(I)+CINDER*DEL VEC(I)=VEC(I)+CINDER*DEL/DABS(DX(I)) GOTO 690 680 X(I)=SAVE 690 IF(NZIP-1)1340,700,700 700 IF(DABS(VEC(I))-ACK)750,710,710 710 DX(I)=ACK*DABS(DX(I)) VEC(I)=VEC(I)/ACK OLDVEC(I)=OLDVEC(I)/ACK DO 720J=1,MOSQUE 720 ERR(I,J)=ERR(I,J)/ACK IF(NTRACE)750,750,730 730 WRITE(6,740)I,DX(I) 740 FORMAT(' STEP SIZE'I3,' INCREASED TO 'E12.5) 750 SUMO=0.0 SUMV=0.0 DO 760J=1,NV SUMO=SUMO+OLDVEC(J)**2 760 SUMV=SUMV+VEC(J)**2 IF(SUMO*SUMV)1340,1340,770 770 SUMO=DSQRT(SUMO) SUMV=DSQRT(SUMV) COSINE=0.0 DO 780J=1,NV 780 COSINE=COSINE+(OLDVEC(J)/SUMO)*(VEC(J)/SUMV) IF(NZIP-1)1340,790,800 790 IF(NACK-NACTIV)1340,820,820 800 IF(NACK-NACTIV)820,810,810 810 IF(NZIP-NCOMP)820,830,830 820 IF(COSINE-COMPAR)1340,830,830 C C SIMON SAYS, TAKE AS MANY GIANT STEPS AS POSSIBLE... C 830 IF(NTRACE)860,860,840 840 WRITE(6,850)CHIOLD,(VEC(J),J=1,I) 850 FORMAT(' CHISQ ='E15.8,5X14HNO. OF STEPS =10F9.2/(42X10F9.2)) 860 NGIANT=0 NTRY=0 NRETRY=0 KL=1 NOSC=NOSC+1 IF(NOSC-MOSQUE)890,890,870 870 NOSC=MOSQUE DO 880K=2,MOSQUE CHIOSC(K-1)=CHIOSC(K) DO 880J=1,NV XOSC(J,K-1)=XOSC(J,K) 880 ERR(J,K-1)=ERR(J,K) 890 DO 900J=1,NV XOSC(J,NOSC)=X(J) 900 ERR(J,NOSC)=VEC(J)/SUMV CHIOSC(NOSC)=CHIOLD IF(NOSC-3)960,910,910 C C SEARCH FOR A PREVIOUS SUCCESSFUL GIANT STEP IN A DIRECTION MORE C NEARLY PARALLEL TO THE DIRECTION OF THE PROPOSED STEP THAN WAS THE C IMMEDIATELY PREVIOUS ONE. C 910 COXCOM=0.0 DO 920J=1,NV 920 COXCOM=COXCOM+ERR(J,NOSC)*ERR(J,NOSC-1) NAH=NOSC-2 930 NTRY=0 DO 950K=KL,NAH NRETRY=NAH-K COSINE=0.0 DO 940J=1,NV 940 COSINE=COSINE+ERR(J,NOSC)*ERR(J,K) IF(COSINE-COXCOM)950,950,970 950 CONTINUE 960 CHIBAK=CHI(I) GOTO 1020 970 NTRY=1 KL=K+1 IF(NTRACE)1000,1000,980 980 NT=NOSC-K WRITE(6,990)NT 990 FORMAT(/1X8H********5X34HPOSSIBLEOSCILLATIONWITHPERIODI2, +' DETECTED.') 1000 DO 1010J=1,NV SALVO(J)=TRIAL(J) 1010 TRIAL(J)=(X(J)-XOSC(J,K))/ACK CHIBAK=CHIOLD+(CHIOSC(K)-CHIOLD)/ACK C 1020 DO 1040J=1,NV XSAVE(J)=X(J) TRIAL(J)=ACK*TRIAL(J) IF(MASK(J))1040,1030,1040 1030 TRAAL=X(J)+TRIAL(J) TRAIL=DMIN1(TRAAL,XMAX(J)) X(J)=DMAX1(TRAIL,XMIN(J)) 1040 CONTINUE JOCK=0 JVARY=0 CALL FUNK NF=NF+1 IF(CHISQ-CHIOLD)1050,1080,1080 1050 CHIBAK=CHIOLD CHIOLD=CHISQ NGIANT=NGIANT+1 IF(NTRACE)1020,1020,1060 1060 WRITE(6,1070)CHISQ,(X(J),J=1,NV) 1070 FORMAT(' CHISQ ='E15.8/' X(I)....'/(10(1XE12.5))) GOTO 1020 C 1080 IF(NRETRY)1100,1100,1090 1090 IF(NGIANT)1150,1150,1100 1100 CINDER=(0.5/ACK)*(ACK**2*CHIBAK-(ACK**2-1.0)*CHIOLD-CHISQ)/(ACK* +CHIBAK-(ACK+1.0)*CHIOLD+CHISQ) DO 1120J=1,NV IF(MASK(J))1120,1110,1120 1110 CANDER=XSAVE(J)+CINDER*TRIAL(J) XIJ=DMIN1(CANDER,XMAX(J)) X(J)=DMAX1(XIJ,XMIN(J)) 1120 CONTINUE JOCK=0 JVARY=0 CALL FUNK NF=NF+1 IF(CHISQ-CHIOLD)1280,1130,1130 1130 IF(NGIANT)1170,1140,1170 1140 IF(NTRY)1150,1170,1150 1150 DO 1160J=1,NV TRIAL(J)=SALVO(J) 1160 X(J)=XSAVE(J) GOTO 1190 1170 DO 1180J=1,NV TRIAL(J)=TRIAL(J)/ACK 1180 X(J)=XSAVE(J) 1190 IF(NTRACE)1240,1240,1200 1200 WRITE(6,1210)CHIOLD,NGIANT 1210 FORMAT(/' CHISQ ='E15.8,' AFTER'I3,' GIANT STEPS.') WRITE(6,1220)(X(J),J=1,NV) 1220 FORMAT(' X(I)....'/(10(1XE12.5))) WRITE(6,1230) 1230 FORMAT(/) 1240 IF(NGIANT)1250,1250,1310 1250 IF(NRETRY)1260,1260,930 1260 IF(NTRY)1270,1330,1270 1270 NTRY=0 GOTO 960 C 1280 CHIOLD=CHISQ JOCK=1 IF(NTRACE)1310,1310,1290 1290 STEPS=FLOAT(NGIANT)+CINDER WRITE(6,1300)CHIOLD,STEPS 1300 FORMAT(/' CHISQ ='E15.8,' AFTER'F6.1,' GIANT STEPS.') WRITE(6,1220)(X(J),J=1,NV) WRITE(6,1230) 1310 IF(NTRY)1320,380,1320 1320 NOSC=0 GOTO 380 1330 NOSC=MAX0(NOSC-1,0) 1340 CHI(I)=CHIOLD 1350 CONTINUE C C ANOTHER CYCLE THROUGH THE VARIABLES HAS BEEN COMPLETED. C PRINT ANOTHER LINE OF TRACES. C IF(NTRACE)1370,1370,1360 1360 WRITE(6,850)CHIOLD,(VEC(J),J=1,NV) 1370 CONTINUE 1380 IF(NZIP)1420,1390,1420 1390 IF(NTRACE)1420,1420,1400 1400 WRITE(6,1220)(X(J),J=1,NV) WRITE(6,1410) 1410 FORMAT(' ') 1420 NZIP=NZIP+1 GOTO 390 C C A MINIMUM HAS BEEN FOUND. PRINT THE REMAINING TRACES. C 1430 IF(NTRACE)1450,1450,1440 1440 WRITE(6,850)CHIOLD,(VEC(J),J=1,I) 1450 IF(NTRACE)1470,1470,1460 1460 WRITE(6,1220)(X(J),J=1,NV) WRITE(6,1230) C C DECREASE THE SIZE OF THE STEPS FOR ALL VARIABLES. C 1470 CONTINUE 1480 NOSC=0 NGATE=1 DO 1520J=1,NV IF(MASK(J))1520,1490,1520 1490 IF(NFLAT(J))1500,1500,1520 1500 IF(DABS(DX(J))-DABS(DELMIN(J)))1520,1520,1510 1510 NGATE=0 1520 DX(J)=DX(J)/RATIO IF(NGATE)1530,1530,1600 1530 IF(NTRACE)1570,1570,1540 1540 WRITE(6,1550)(DX(J),J=1,NV) 1550 FORMAT(60(1X1H*)//' STEP SIZES REDUCED TO....'//(10(1XE12.5))) WRITE(6,1560) 1560 FORMAT(//) 1570 GOTO 380 1580 WRITE(6,1590)(DX(J),J=1,NV) 1590 FORMAT(' OPERAT. TERM.'//10(1X,E12.5)) C 1600 CHISQ=CHIOLD IF(NTRACE)1630,1610,1610 1610 WRITE(6,1620)NF 1620 FORMAT(//1X,I5,' FUNCTION COMPUTATIONS ') 1630 CONTINUE 1640 IF(IABS(MATRIX-100)-50)1650,1650,2160 1650 IF(NACTIV-NV)2160,1660,2160 C C COMPUTE THE STANDARD ERRORS AND THE CORRELATIONS. C 1660 FAC=RATIO**(MATRIX-100) DO 1680I=1,NV DX(I)=DABS(FAC*DX(I)) XSAVE(I)=X(I) JVARY=0 DO 1670J=1,2 X(I)=XSAVE(I)+DX(I) CALL FUNK NF=NF+1 JVARY=I SECOND(1,J)=CHISQ 1670 DX(I)=-DX(I) X(I)=XSAVE(I) 1680 ERR(I,I)=(SECOND(1,1)-2.0*CHIOLD+SECOND(1,2))/DX(I)**2 DO 1710I=2,NV IM=I-1 DO 1710J=1,IM DO 1700K=1,2 X(I)=XSAVE(I)+DX(I) JVARY=0 DO 1690L=1,2 X(J)=XSAVE(J)+DX(J) CALL FUNK NF=NF+1 JVARY=J SECOND(K,L)=CHISQ X(J)=XSAVE(J) 1690 DX(J)=-DX(J) X(I)=XSAVE(I) 1700 DX(I)=-DX(I) ERR(I,J)=0.25*(SECOND(1,1)-SECOND(1,2)-SECOND(2,1)+SECOND(2,2))/ +DABS(DX(I)*DX(J)) 1710 ERR(J,I)=ERR(I,J) IF(NTRACE)1770,1720,1720 1720 WRITE(6,1730) 1730 FORMAT('1SIZES OF INCREMENTS TO BE USED BELOW....') WRITE(6,1740)(DX(J),J=1,NV) 1740 FORMAT(/(10(1XE12.5))) WRITE(6,1750) 1750 FORMAT(/////' MATRIX OF THE SECOND PARTIAL DERIVATIVES....'/) DO 1760I=1,NV 1760 WRITE(6,1740)(ERR(I,J),J=1,I) 1770 DO 1780I=1,NV DO 1780J=1,I IF(ERR(I,J))1780,1790,1780 1780 CONTINUE GOTO 1810 1790 WRITE(6,1800) 1800 FORMAT(//// +' THE ABOVE MATRIX CONTAINS ONE OR MORE ZEROES. A LARGER VALUE OF + -MATRIX- SHOULD BE TRIED, TO SEE IF THEY ARE LEGITIMATE. ') C C INVERT THE MATRIX USING SYMINV2 (COMM. OF THE A.C.M. 6, P. 67). C 1810 DET=1.0 DETLOG=0. DO 1820J=1,NV 1820 SALVO(J)=1.0 DO 1970I=1,NV BIGAJJ=0.0 DO 1850J=1,NV IF(SALVO(J))1830,1850,1830 1830 IF(DABS(ERR(J,J))-BIGAJJ)1850,1850,1840 1840 BIGAJJ=DABS(ERR(J,J)) K=J 1850 CONTINUE IF(BIGAJJ)1870,1860,1870 1860 DET=0.0 GOTO 1980 1870 SALVO(K)=0.0 IF(DET.GT.1.E35.OR.DET.LT.1.E-35)GOTO1871 GOTO 1875 1871 WRITE(6,1872) 1872 FORMAT(///,'DETERMINANT COMPUTATION WILL CAUSE OVERFLOW.') GOTO 2160 1875 DET=DET*ERR(K,K) DETLOG=DETLOG+DLOG(DABS(ERR(K,K)))/2.303 TRIAL(K)=1.0/ERR(K,K) ERR(K,K)=0.0 XSAVE(K)=1.0 M=K-1 IF(M)1910,1910,1880 1880 DO 1900J=1,M XSAVE(J)=ERR(K,J) TRIAL(J)=ERR(K,J)*TRIAL(K) IF(SALVO(J))1860,1900,1890 1890 TRIAL(J)=-TRIAL(J) 1900 ERR(K,J)=0.0 1910 M=K+1 IF(M-NV)1920,1920,1960 1920 DO 1950J=M,NV XSAVE(J)=ERR(J,K) IF(SALVO(J))1860,1930,1940 1930 XSAVE(J)=-XSAVE(J) 1940 TRIAL(J)=-ERR(J,K)*TRIAL(K) 1950 ERR(J,K)=0.0 1960 DO 1970J=1,NV DO 1970K=J,NV 1970 ERR(K,J)=ERR(K,J)+XSAVE(J)*TRIAL(K) IF(DET)2000,1980,2020 1980 WRITE(6,1990) 1990 FORMAT(//// +' ERROR MATRIX IS SINGULAR. -MATRIX- SHOULD PROBABLY BE INCREASED +. '/////) GOTO 2150 2000 WRITE(6,2010) 2010 FORMAT(//// +' ERROR MATRIX IS NEGATIVE DEFINITE. -MATRIX- SHOULD PROBABLY BE +DECREASED.') 2020 IF(NTRACE)2050,2030,2030 2030 WRITE(6,2040)DET,DETLOG 2040 FORMAT(////' DETERMINANT OF ABOVE MATRIX = 'E12.5, E12.5) C +10X14HLOG10F(DET)=E12.5) C ****** THE ABOVE LINE WAS COMMENTED OUT TO MAKE THE C ****** COMPILIER HAPPY. IT BELONGS WITH LINE 2040 C *************************************************** C 2050 DO 2090I=1,NV DO 2060J=1,I ERR(I,J)=ERR(I,J)*2.0 2060 ERR(J,I)=ERR(I,J) IF(ERR(I,I))2070,2070,2090 2070 WRITE(6,2080)ERR(I,I) 2080 FORMAT(///' NEGATIVE OR ZERO MEAN SQUARE ERROR ENCOUNTERED...', +3X,E15.8/' -MATRIX- SHOULD PROBABLY BE DECREASED.'///) 2090 XSAVE(I)=DSIGN(DSQRT(DABS(ERR(I,I))),ERR(I,I)) IF(NTRACE)2160,2100,2100 2100 WRITE(6,2110) 2110 FORMAT(/////' STANDARD ERRORS....') WRITE(6,1740)(XSAVE(J),J=1,NV) WRITE(6,2120) 2120 FORMAT(/////' LOWER TRIANGLE OF THE CORRELATION MATRIX....'/) DO 2140I=2,NV IM=I-1 DO 2130J=1,IM 2130 TRIAL(J)=ERR(I,J)/DABS(XSAVE(I)*XSAVE(J)) 2140 WRITE(6,1740)(TRIAL(J),J=1,IM) 2150 WRITE(6,1620)NF C 2160 CONTINUE 2190 JVARY=0 CALL FUNK IF(NTRACE)2230,2200,2200 2200 WRITE(6,2210)(X(J),J=1,NV) 2210 FORMAT(///10X24HFINAL VALUES OF X(I)....//(7(1XE16.9))) WRITE(6,2220)CHISQ 2220 FORMAT(//' FINAL VALUE OF CHISQ = 'E15.8//) 2230 RETURN END