C*************************************************************
C ****
C **** CPTFIT - CPT Model FOR LS FORTRAN 29-oct-1996 ****
C **** by Michael H. Birnbaum, copyright
C all rights reserved.
C **** Updated Dec 15, 1997 Updated 5-3-98 ****
C*************************************************************
C
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 You will also cite Tversky & Kahneman (1992) CPT model
C
C*** Note: This version includes Chandler's STEPIT Subroutine.
C*** If you use STEPIT, also cite Chandler (1969).
C****************************************************************
C THIS PROGRAM WILL FIT CPT TO ANY EXPERIMENT
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 W(p) = F(p) 2 PARAMETER Weighting Function
C parameters are Gamma and 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
C DATA FORMAT (Each CARD is a line of 80 columns ) in
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 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 Some parameters should be constrained to be positive
C UNLIKE Birnbaum & Chavez,
C 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
values will be adjacent
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
C **************************************************************
C
PROGRAM CPTFIT
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 = 'CPTFIT.DAT', STATUS = 'OLD')
OPEN (6, FILE = 'CPTFIT.OUT', STATUS = 'NEW')
OPEN (7, FILE = 'SUBJECTS.DAT', STATUS = 'OLD')
OPEN (8, FILE = 'CPTFIT.PAR', STATUS = 'NEW')
OPEN (9, FILE = 'CPTFIT.RES', STATUS = 'NEW')
OPEN(10, FILE = 'CPTFIT.PRE', STATUS = 'NEW')
C
C READ IN THE NUMBER OF TRIALS, GAMBLE PAIRS = NPAIRS
READ(5,99)NPAIRS,MID,WWLIK
READ(5,99)NSUBS
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=103
C READ IN INITIAL ESTIMATES OF THE X VALUES
C TRANSFER TABLE
C GAMMA, BETA, A, B, C
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 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 LAST THREE CARDS BLANK=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)
C This section checks probabilities add to 1
SUMPX = 0.0
DO 107 J=1,NL
SUMPX=SUMPX+PX(IP,J)
107 CONTINUE
SUMPY=0.0
DO 108 K=1,NR
SUMPY=SUMPY+PY(IP,K)
108 CONTINUE
IF (SUMPY .NE. 1.0 .OR. SUMPX .NE. 1.0) GO TO 109
GO TO 110
109 WRITE(6,106) IP,SUMPX,SUMPY
106 FORMAT(' PROBS DO NOT ADD TO 1.0 IN PROBLEM ',I4,2X,2F5.3)
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)
W1 = (C(ISUB)*(.25)**GAMMA(ISUB))/(C(ISUB)*(.25)**GAMMA(ISUB) + (.75)**GAMMA(ISUB))
W2 = (C(ISUB)*(.50)**GAMMA(ISUB))/(C(ISUB)*(.50)**GAMMA(ISUB) + (.50)**GAMMA(ISUB)) -W1
W3 = (C(ISUB)*(.75)**GAMMA(ISUB))/(C(ISUB)*(.75)**GAMMA(ISUB) + (.25)**GAMMA(ISUB))
W3 = W3 - (C(ISUB)*(.50)**GAMMA(ISUB))/(C(ISUB)*(.50)**GAMMA(ISUB) + (.50)**GAMMA(ISUB))
W4 = 1 - W3 - W2 - W1
C PRINT WEIGHTS FOR EACH SUBJECT
WRITE(8,723) W4,W3,W2,W1,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 2 PARAMETER GENERAL CPT 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 ---This 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 **********************************************************
C NOW CALL STEPIT
CALL STEPIT
C*************************----------------
WRITE(6,941)
941 FORMAT(' FIT FOR AVERAGED DATA')
WRITE(6,801)
801 FORMAT(' FIT OF CUMULATIVE PROSPECT THEORY')
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(T20,' DATA =',F8.2,5X,' PRED =',F8.2,5X,' RESID=',5X,F8.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)
W1 = (C(ISUB)*(.25)**GAMMA(ISUB))/(C(ISUB)*(.25)**GAMMA(ISUB) + (.75)**GAMMA(ISUB))
W2 = (C(ISUB)*(.50)**GAMMA(ISUB))/(C(ISUB)*(.50)**GAMMA(ISUB) + (.50)**GAMMA(ISUB)) -W1
W3 = (C(ISUB)*(.75)**GAMMA(ISUB))/(C(ISUB)*(.75)**GAMMA(ISUB) + (.25)**GAMMA(ISUB))
W3 = W3 - (C(ISUB)*(.50)**GAMMA(ISUB))/(C(ISUB)*(.50)**GAMMA(ISUB) + (.50)**GAMMA(ISUB))
W4 = 1 - W3 - W2 - W1
C PRINT WEIGHTS FOR GROUP DATA INTO SUBJECT FILE SUBJECT
WRITE(8,723) W4,W3,W2,W1,NSP1
WRITE(9,FMT3) (RESD(ISUB,IP),IP=1,NPAIRS),NSP1
WRITE(10,FMT3) (PD(ISUB,IP),IP=1,NPAIRS),NSP1
C
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 **** ****
C ********************************************************************
C
C HERE IS FUNK--WHICH FITS GEN CUMULATIVE PROSPECT THEORY TO JUDGMENTS
C OF PREFERENCES BETWEEN PAIRS OF GAMBLES
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--------------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
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
SUMR = 0.0
C PC IS DECUMULATIVE PROBABILITY
C NOTE THAT IK WILL GO IN REVERSE ORDER
C GAMBLE ON THE LEFT *****************************LEFT
IKM = NL + 1
W(IKM) = 0
PC = 0.0
DO 520 I = 1,NL
WTEMP = (C(ISUB)*(PC**GAMMA(ISUB)))
WTEMP = WTEMP/(C(ISUB)*PC**GAMMA(ISUB)+(1.0-PC)**GAMMA(ISUB))
IK = NL+1-I
PC = PX(IP,IK) + PC
IKK = IK +1
IF (PC .GE. 1.0) GO TO 518
W(IK) =(C(ISUB)*(PC**GAMMA(ISUB)))
W(IK) = W(IK)/(C(ISUB)*PC**GAMMA(ISUB)+(1.0-PC)**GAMMA(ISUB))-WTEMP
GO TO 519
518 W(IK) = 1.0 - WTEMP
519 CONTINUE
SUML = SUML + W(IK)*XO(IP,IK)**BETA(ISUB)
UTL(IP) = SUML
520 CONTINUE
C GAMBLE ON THE RIGHT *****************************RIGHT
IKM = NR + 1
W(IKM) = 0
PC = 0.0
DO 530 I = 1,NR
WTEMP = (C(ISUB)*(PC**GAMMA(ISUB)))
WTEMP = WTEMP/(C(ISUB)*PC**GAMMA(ISUB)+(1.0-PC)**GAMMA(ISUB))
IK = NR+1-I
PC = PY(IP,IK) + PC
IKK = IK +1
IF (PC .GE. 1.0) GO TO 528
W(IK) =(C(ISUB)*(PC**GAMMA(ISUB)))
W(IK) = W(IK)/(C(ISUB)*PC**GAMMA(ISUB)+(1.0-PC)**GAMMA(ISUB))-WTEMP
GO TO 529
528 W(IK) = 1.0 - WTEMP
529 CONTINUE
SUMR = SUMR + W(IK)*YO(IP,IK)**BETA(ISUB)
UTR(IP) =SUMR
530 CONTINUE
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