SUBROUTINE CHIBNDRY C This subroutine estimates the uncertainties in selected parameters C for a night's binary star data as fitted with BINFIT. The approach C taken here is to increment one parameter, then optimize the fit by C allowing the other selected parameters to vary (while holding fixed C the parameters that were not varied in the original fit). This C process is continued until the reduced chi-square has increased by C some amount. The amount of chi-square increase desired depends on C the number of degrees of freedom and the confidence level desired. C J.T. Armstrong 12 Feb. 1991 c INCLUDE 'BINFIT.INC' CHARACTER*1 CHOICE REAL*8 ALAMDA, PROBCHI, OCHISQ, DELCHISQ, NCHISQ REAL*8 FREE, COVAR(MXPARM,MXPARM), ALPHA(MXPARM,MXPARM) REAL*8 DP1D(MXPARM),T,A0,A1,B1,B2,XP,EXPCHISQ REAL*8 LOWVALUE, HIGHVALUE INTEGER*4 ITST, LOOPCT, ALLVARY, KK, LISTASV(MXPARM) INTEGER*4 NVPARSV, FXM CHARACTER*64 OUTNAM LOGICAL*1 REPEAT, INDEPT(MXPARM), TOOBIG c The following four lines are borrowed from clisterr.f (c.hummel) CHARACTER*80 HISTRY CHARACTER*1 BLANK INTEGER*4 LEN1 DATA BLANK/' '/ DATA A0,A1,B1,B2/2.30753,0.27061,0.99229,0.04481/ METHOD = 'Marquardt-Levenberg' NVPARSV=NVPARTOT C Fill the array of amounts to vary the parameters by ALAMDA = -1. CALL BINMRQMI(COVAR,ALPHA,ALAMDA,REPEAT) ALAMDA = 0. CALL BINMRQMI(COVAR,ALPHA,ALAMDA,REPEAT) DO I = 1,MXPARM DP1D(I) = 0. END DO DO I = 1,MXPARM IF (ALPHA(I,I).NE.0.) THEN DP1D(LISTA(I)) = 1./SQRT(ALPHA(I,I)) END IF END DO C WRITE(OUTC,*) ' ALPHA array:' C DO I = 1,NVPARTOT C WRITE(OUTC,'(8(G8.2,2X))') (ALPHA(I,J),J=1,NVPARTOT) C END DO C Find parameters to be varied that depend on NONE of the other C parameters to be varied. DO I = 1,NVPARTOT INDEPT(LISTA(I)) = .TRUE. DO J = 1,I-1 IF(ABS(ALPHA(I,J)).GT.1D-6*ABS(ALPHA(I,I))) THEN INDEPT(LISTA(I)) = .FALSE. END IF END DO DO J = I+1,NVPARTOT IF(ABS(ALPHA(I,J)).GT.1D-6*ABS(ALPHA(I,I))) THEN INDEPT(LISTA(I)) = .FALSE. END IF END DO END DO D WRITE(OUTC,*) ' List of parameters not dependent on others:' D DO I = 1,NVPARTOT D IF(INDEPT(LISTA(I))) WRITE(OUTC,*) LISTA(I) D END DO C Get the increase in CHISQ desired WRITE(OUTC,*) ' Desired confidence level (> 50 percent): ' 8 READ(INC,*) CONF IF (CONF.GT.1.) CONF = CONF/100. ! if CONF given in percent IF (CONF.LT.0.5) THEN WRITE(OUTC,*) ' Confidence level must be > 50 percent: ' GO TO 8 END IF C Calculate CHISQ corresponding to indicated confidence level. C See Abramowitz & Stegun 26.2.22 and 26.4.17. T = DSQRT(-2.*LOG(1.-CONF)) XP = T - (A0+A1*T)/(1.+B1*T+B2*T*T) FREE = NZTOT-NVPARTOT NCHISQ = (1.-(2./(9.*FREE))+XP*DSQRT(2./(9.*FREE)))**3 WRITE(OUTC,'(A,F7.4)') ' Corresponding reduced chi squared is ', 1 NCHISQ C Calculate CHISQ expected for this number of degrees of freedom. C (Same calculation as above, with CONF = 0.5 => XP = 0.) EXPCHISQ = (1.-(2./(9.*FREE)))**3 DELCHISQ = NCHISQ-EXPCHISQ IF(TCHISQ.GT.EXPCHISQ) THEN NCHISQ = NCHISQ*TCHISQ/EXPCHISQ WRITE(OUTC,'(A,F7.4)') 1 ' Chi squared at confidence level rescaled to ',NCHISQ END IF C Calculate the confidence intervals from the diagonal elements of the C covariance matrix (see, e.g., Numerical Recipes, p. 536, eq. 14.5.4). C In the approximation of chi^2 forming a parabolic hypersurface in C parameter space, this will give the uncertainty in each parameter C when all the others are free to adjust to the perturbation. DO I=1,NVPARTOT PUNCERT3D(LISTA(I),0) 1 = SQRT(FREE*DELCHISQ*COVAR(LISTA(I),LISTA(I))) END DO C For each I, check to see if the corresponding parameter was varied. C If yes, change it a little, then hold it constant while varying all C the other variable parameters. WRITE(OUTC,'(A,F6.3,A)') 1' Fractional change in CHISQ for convergence (suggested value: ', 2 1./REAL(FREE),'): ' READ(INC,*) DELCHISQ 5 WRITE(OUTC,'(A,A)') ' Enter filename for listing points at ', 1 'chi^2 boundary: ' READ(INC,'(A)') OUTNAM IF(OUTNAM.NE.' ') THEN OPEN(UNIT=OUTDAT,FILE=OUTNAM,STATUS='NEW',FORM='FORMATTED', 1 ERR=5) END IF CALL MKHIST('ERRFIT','4 Mar 1991',HISTRY) L = LEN1(HISTRY) WRITE(OUTDAT,1000) HISTRY(1:L) WRITE (OUTDAT,1000) ' Input data file(s): ' DO I = 1, NFILES L = LEN1(INFILE(I)) WRITE(OUTDAT,1000) ' ',INFILE(I)(1:L) END DO L = LEN1(MODDSN) WRITE (OUTDAT,1000) ' Final model file: ',MODDSN(1:L) WRITE (OUTDAT,1001) STARNAME,NDAY,NYEAR WRITE (OUTDAT,1000) ' Filter Width Reduced Degrees of' WRITE (OUTDAT,1000) ' Chi^2 Freedom' DO IG=1,NFILT WRITE (OUTDAT,1002) LAMBDA0(IG),DLAMBDA(IG),FCHISQ(IG), 1 NZVIS(IG)-NVPAR(IG) END DO WRITE (OUTDAT,1003) 'Reduced chi squared for all data: ',TCHISQ WRITE (OUTDAT,1004) 'Number of data points: ',NZTOT WRITE (OUTDAT,1004) 'Degrees of freedom: ', 1 NZTOT-NVPARTOT WRITE (OUTDAT,1000) ' Final model:' CALL WRITEHD(OUTDAT,BLANK) DO IC = 1, NCOMP DO IG = 1,NFILT CALL CWRITEMD(OUTDAT,BLANK,PARM,IC,IG) END DO END DO 1000 FORMAT (A,A,A,A) 1001 FORMAT (' Source: ',A,' Day: ',I4,' Year: ',I5) 1002 FORMAT (' ',3X,F5.1,' nm',2X,F3.0,' nm',3X,F8.3,3X,I4) 1003 FORMAT (' ',(A), F8.3 ) 1004 FORMAT (' ',(A), I4 ) 1005 FORMAT (' ',(A), F10.0 ) DO I = 1,MXPARM NVPARTOT=NVPARSV-1 IF (VARY1D(I)) THEN D WRITE(OUTC,'(A,I3)') ' Varying parameter ',I VARY1D(I) = .FALSE. C Make list of parameters to fit in this step K = 0 DO J = 1,MXPARM IF(VARY1D(J).EQ..TRUE.) THEN C --if(I not a color or J not a mag. diff) then C --include J in list of variable parameters FXM=MXFILT*MXMOD IF( .NOT.(I.GT.5*FXM .AND. I.LE.6*FXM) 1 .OR. 2 .NOT.(J.GT.4*FXM .AND. J.LE.5*FXM) ) THEN K=K+1 LISTA(K)=J D WRITE(OUTC,'(A,I3,A,I3)') D 1 ' LISTA(',K,') = ',LISTA(K) ELSE NVPARTOT=NVPARTOT-1 END IF END IF C ---These lines are the old code C IF(VARY1D(J).EQ..TRUE.) THEN C K=K+1 C LISTA(K)=J C END IF END DO C Vary PARM by +DELTAP (when J=1) or by -DELTAP (when J=-1) DO J=1,-1,-2 DO K = 1,MXPARM PARM1D(K) = PARMST1D(K) END DO C WRITE(OUTC,'(A,G11.4,A,G11.4)') C 1 ' Starting at ',PARM1D(I),'; increment by ', C 1 DP1D(I) LOWVALUE = PARM1D(I) TOOBIG = .FALSE. PARM1D(I) = PARM1D(I) + J*DP1D(I) C Get best fit of other NVPARSV-1 parameters 10 IF (NVPARTOT.GT.0 .AND. .NOT.INDEPT(I)) THEN ALAMDA = -1. ! Return here for repeat or for next step ITER = 0 ITST = 0 C WRITE(OUTC,*) 'Initializing BINMRQMI' CALL BINMRQMI(COVAR,ALPHA,ALAMDA,REPEAT) 100 OCHISQ = TCHISQ ! Loop to here if no convergence ITER = ITER + 1 C WRITE(OUTC,'(A,G11.4)') C 1 ' CHIBNDRY: OCHISQ = ',OCHISQ C Call BINMRQMI and observe change in CHISQ CALL BINMRQMI(COVAR,ALPHA,ALAMDA,REPEAT) C CALL CAFOUT(1) C WRITE(OUTC,'(A,I3,A,G11.4,A,F8.4)') C 1 ' For PARM1D(',I,') = ',PARM1D(I), C 2 ', TCHISQ = ',TCHISQ C --If chi^2 < NCHISQ, don't repeat BINMRQMI IF( (NCHISQ-TCHISQ)/NCHISQ.GT.DELCHISQ ) THEN ITST = 2 C WRITE(OUTC,*) 'Low chi^2; ITST set to 2' END IF C --If NCHISQ has been bracketed and bracketing C interval is small, don't repeat BINMRQMI IF( ABS(HIGHVALUE-LOWVALUE).LT..01*PUNCERT3D(I,0) 1 .AND. TOOBIG) THEN ITST = 2 C WRITE(OUTC,*) 'Close to chi^2 boundary;'// C 1 ' ITST set to 2' END IF C --If BINMRQMI reduced chi^2 and chi^2 is close C enough to the previous value, increment ITST IF( ABS(OCHISQ/TCHISQ-1.).LT.DELCHISQ 1 .AND. .NOT.REPEAT ) THEN ITST = ITST + 1 C WRITE(OUTC,'(A,G15.8,A,I2)')' OCHISQ/TCHISQ=', C 1 OCHISQ/TCHISQ,'; add one to ITST, now ',ITST END IF C --If ITST isn't big enough, repeat BINMRQMI; C otherwise, go on to test chi^2 against NCHISQ IF (ITST.LT.2) THEN C WRITE(OUTC,*) ' Try again:' GO TO 100 END IF ELSE CALL CADJMOD CALL BINVIS CALL CAGRFAC END IF C Test for convergence to desired TCHISQ value C Chi squared still too small IF ((NCHISQ-TCHISQ)/NCHISQ.GT.DELCHISQ 1 .AND. .NOT.TOOBIG) THEN C WRITE(OUTC,'(A,I3,A,G11.4,A,G11.4,A,G11.4)') C 1 ' CHISQ TOO SMALL; PARM1D(',I,')=',PARM1D(I), C 2 ' DP1D=',DP1D(I),' TCHISQ=',TCHISQ LOWVALUE = PARM1D(I) DP1D(I) = 2.*DP1D(I) PARM1D(I) = PARM1D(I) + J*DP1D(I)/2. GO TO 10 C Chi squared got too big ELSE IF ((NCHISQ-TCHISQ)/NCHISQ.LT.-DELCHISQ 1 .AND. .NOT.TOOBIG) THEN C WRITE(OUTC,'(A,I3,A,G11.4,A,G11.4,A,G11.4)') C 1 ' CHISQ TOO BIG; PARM1D(',I,')=',PARM1D(I), C 2 ' DP1D=',DP1D(I),' TCHISQ=',TCHISQ HIGHVALUE = PARM1D(I) PARM1D(I) = (LOWVALUE+HIGHVALUE)/2. TOOBIG = .TRUE. GO TO 10 C Chi squared too small ELSE IF ((NCHISQ-TCHISQ)/NCHISQ.GT.DELCHISQ 1 .AND. TOOBIG) THEN C WRITE(OUTC,'(A,I3,A,G11.4,A,G11.4)') C 1 ' CHISQ TOO SMALL; PARM1D(',I,')=',PARM1D(I), C 2 ' TCHISQ=',TCHISQ LOWVALUE = PARM1D(I) PARM1D(I) = (LOWVALUE+HIGHVALUE)/2. IF(ABS(HIGHVALUE-LOWVALUE).GT..01*PUNCERT3D(I,0)) 1 THEN GO TO 10 END IF C Chi squared got too big ELSE IF ((NCHISQ-TCHISQ)/NCHISQ.LT.-DELCHISQ 1 .AND. TOOBIG) THEN C WRITE(OUTC,'(A,I3,A,G11.4,A,G11.4)') C 1 ' CHISQ TOO BIG; PARM1D(',I,')=',PARM1D(I), C 2 ' TCHISQ=',TCHISQ HIGHVALUE = PARM1D(I) PARM1D(I) = (LOWVALUE+HIGHVALUE)/2. IF(ABS(HIGHVALUE-LOWVALUE).GT..01*PUNCERT3D(I,0)) 1 THEN GO TO 10 END IF END IF C Chi squared is juuuuuuuust right PUNCERT3D(I,J) = PARM1D(I) - PARMST1D(I) DP1D(I) = PARM1D(I) - PARMST1D(I) WRITE(OUTC,'(A,I3,A,G11.4)') 1 ' DONE WITH THIS VARIABLE; PARM1D(',I,')= ', 2 PARM1D(I) WRITE(OUTDAT,'(A,I3,A,G11.4)') 1 ' DONE WITH THIS VARIABLE; PARM1D(',I,')= ', 2 PARM1D(I) WRITE(OUTDAT,'(A,F8.4,A)') 1 ' CHISQ = ',NCHISQ,' at these values: ' CALL CWRITEMD(OUTDAT,' ',PARM,0,0) END DO VARY1D(I) = .TRUE. END IF END DO CLOSE(OUTDAT) RETURN END