SUBROUTINE MRQSHELL C This subroutine estimates selected parameters for a night's binary C star data using the Levenberg-Marquardt method. C J.T. Armstrong 4 Mar. 1991 INCLUDE 'BINFIT.INC' REAL*8 ALAMDA, OCHISQ, DELCHISQ, NCHISQ REAL*8 FREE, COVAR(MXPARM,MXPARM), ALPHA(MXPARM,MXPARM) REAL*8 DP1D(MXPARM),T,A0,A1,B1,B2,XP,EXPCHISQ,DCHISQ REAL*8 DUMMY(MXPARM),COVARI,COVARJ INTEGER*4 ITST, LOOPCT, ALLVARY, KK, LISTASV(MXPARM), WHERE INTEGER*4 IERR CHARACTER*12 COVFILE LOGICAL*1 REPEAT DATA A0,A1,B1,B2/2.30753,0.27061,0.99229,0.04481/ METHOD = 'Marquardt-Levenberg' WRITE(OUTC,*) ' Fractional change in CHISQ for convergence: ' READ(INC,*) DELCHISQ C Get best fit to NVPARTOT parameters C WRITE(OUTC,'(A,I3)') ' MRQSHELL: NVPARTOT = ',NVPARTOT 10 IF (NVPARTOT.GT.1) THEN ALAMDA = -1. ! Initialization step ITER = 0 ITST = 0 C WRITE(OUTC,*) ' MRQSHELL: starting model:' C CALL CAFOUT(1) CALL BINMRQMI(COVAR,ALPHA,ALAMDA,REPEAT) 100 OCHISQ = TCHISQ ! Loop to here if no convergence C WRITE(OUTC,*) 'Try again:' ITER = ITER + 1 C WRITE(OUTC,'(A,G11.4)') ' MRQSHELL: OCHISQ = ',OCHISQ C Call BINMRQMI and observe change in CHISQ CALL BINMRQMI(COVAR,ALPHA,ALAMDA,REPEAT) CALL CAFOUT(1) C WRITE(OUTC,'(A,G11.4,3X,G11.4,3X,I2)') C 1 ' MRQSHELL: TCHISQ, OCHISQ, ITST = ',TCHISQ,OCHISQ,ITST IF( .NOT.REPEAT .AND. OCHISQ/TCHISQ-1..GT.0. 1 .AND. OCHISQ/TCHISQ-1..LT.DELCHISQ ) THEN C WRITE(OUTC,'(A,G15.8,A)') C 1 ' OCHISQ/TCHISQ=',OCHISQ/TCHISQ,'; add one to ITST' ITST = ITST + 1 END IF IF (ITST.LT.2) GO TO 100 ELSE CALL CADJMOD CALL BINVIS CALL CAGRFAC END IF ALAMDA = 0. ! To get final COVAR and ALPHA CALL BINMRQMI(COVAR,ALPHA,ALAMDA,REPEAT) C Get the confidence level desired WRITE(OUTC,*) ' Desired confidence level (> 50 percent): ' 200 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 200 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 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 DCHISQ = NCHISQ-EXPCHISQ WRITE(OUTC,'(A,I3,A,F5.3,A,F5.3,A,F4.1,A)') 1 ' Expect Chi^2 (',FREE,' d.o.f.) = ',EXPCHISQ,' (min) or ', 2 NCHISQ,' (',CONF*100.,'% confidence).' cExpect Chi^2 (123 d.o.f.) = 0.988 (min) or 1.234 (99.9% confidence). WRITE(OUTC,'(A,F6.3,A,F6.3,A,F4.1,A)') 1 ' Observe Chi^2(min) = ',TCHISQ,'; seek Chi^2 = ',NCHISQ, 2 ' at ',CONF*100.,'% confidence.' cObserve Chi^2(min) = 1.234; seek Chi^2 = 2.345 at 99.9% confidence. IF(TCHISQ.GT.EXPCHISQ) THEN NCHISQ = NCHISQ*TCHISQ/EXPCHISQ WRITE(OUTC,'(A,G11.4)') ' Scaled NCHISQ = ',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*DCHISQ*COVAR(LISTA(I),LISTA(I))) WRITE(OUTC,'(A,A,A,G11.4)') 1 ' Uncertainty in ',PARMNAME(LISTA(I)),' = ', 2 PUNCERT3D(LISTA(I),0) END DO C Repeat with smaller DELCHISQ if desired WRITE(OUTC,'(A,F7.3,A)') ' Convergence is slow (',TCHISQ,')' IF(QUERY(' Continue (Y or N)?')) THEN DELCHISQ = DELCHISQ/5. GO TO 10 END IF C Write the covariance matrix WRITE(OUTC,*) 1 'Covariance matrix to screen (S), file (F), or both (B)? ' READ(INC,'(A1)') ANSWER CALL CAPS(ANSWER) IF(ANSWER.EQ.'S' .OR. ANSWER.EQ.'B') THEN WRITE(OUTC,1000) STARNAME(1:8),DATE WRITE(OUTC,1001) (PARMNAME(LISTA(I)),I=1,NVPARTOT) DO I=1,NVPARTOT WRITE(OUTC,1002) PARMNAME(LISTA(I)), 1 (COVAR(LISTA(I),LISTA(J)),J=1,NVPARTOT) END DO WRITE(OUTC,1003) WRITE(OUTC,1001) (PARMNAME(LISTA(I)),I=1,NVPARTOT) DO I=1,NVPARTOT COVARI=SQRT(COVAR(LISTA(I),LISTA(I))) DO J=1,NVPARTOT COVARJ=SQRT(COVAR(LISTA(J),LISTA(J))) DUMMY(J)=COVAR(LISTA(I),LISTA(J))/(COVARI*COVARJ) END DO WRITE(OUTC,1002) PARMNAME(LISTA(I)), 1 (DUMMY(J),J=1,NVPARTOT) END DO END IF IF(ANSWER.EQ.'F' .OR. ANSWER.EQ.'B') THEN 300 WRITE(OUTC,*) 'Enter file for covariance matrix output: ' READ(INC,'(A12)') COVFILE OPEN(UNIT=OUTDAT,FILE=COVFILE,STATUS='NEW',IOSTAT=IERR) IF(IERR.NE.0) THEN WRITE(OUTC,*) 'Cannot open file ',COVFILE GO TO 300 END IF WRITE(OUTDAT,1000) STARNAME(1:8),DATE WRITE(OUTDAT,1010) WRITE(OUTDAT,1011) (PARMNAME(LISTA(I)),I=1,NVPARTOT) WRITE(OUTDAT,1010) DO I=1,NVPARTOT WRITE(OUTDAT,1012) PARMNAME(LISTA(I)), 1 (COVAR(LISTA(I),LISTA(J)),J=1,NVPARTOT) WRITE(OUTDAT,1010) END DO WRITE(OUTDAT,1010) WRITE(OUTDAT,1010) WRITE(OUTDAT,1003) WRITE(OUTDAT,1010) WRITE(OUTDAT,1011) (PARMNAME(LISTA(I)),I=1,NVPARTOT) WRITE(OUTDAT,1010) DO I=1,NVPARTOT COVARI=SQRT(COVAR(LISTA(I),LISTA(I))) DO J=1,NVPARTOT COVARJ=SQRT(COVAR(LISTA(J),LISTA(J))) DUMMY(J)=COVAR(LISTA(I),LISTA(J))/(COVARI*COVARJ) END DO WRITE(OUTDAT,1012) PARMNAME(LISTA(I)), 1 (DUMMY(J),J=1,NVPARTOT) WRITE(OUTDAT,1010) END DO CLOSE(OUTDAT) END IF IF(ANSWER.EQ.'S') THEN WRITE(OUTC,*) 'Type to continue.' READ(INC,'(A1)') ANSWER END IF RETURN 1000 FORMAT(' Covariance matrix for ',A8,' on ',I6) 1001 FORMAT(1X,5X, 20(2X,A3,4X)) 1002 FORMAT(1X,A3,2X,20(F7.3,2X) ) 1003 FORMAT(' Correlation matrix') 1010 FORMAT(' ') 1011 FORMAT(1X,5X, 20(3X,A3,4X)) 1012 FORMAT(1X,A3,2X,20(F8.4,2X) ) END