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