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