SUBROUTINE ERREST(A,LISTA,MFIT,COVAR) C---------------------------------------------------------------------- C Subroutine called by SETELL to find the one-sigma points along the C parameters fit by ORBFIT. The program varies each parameter from its C best-fit value and finds the best-fit set of remaining parameters. C The chi-squared value at that point is used to determine whether the C one-sigma point in that parameter has been reached. The uncertainties C are returned in the array UNCERT(I,J), where J=-1 is the uncertainty C below and J=1 is the uncertainty above the best-fit value. By C J.T. Armstrong C 25 Apr 1991 J.T. Armstrong Revised to use error ellipse C parameters C---------------------------------------------------------------------- IMPLICIT UNDEFINED (A-Z) INCLUDE 'PARAMETR.INC' INCLUDE 'ELLIPSE.INC' CHARACTER*1 CHOICE CHARACTER*80 INS REAL*8 ALAMDA, RMOD, THMOD, RAMOD, DECMOD REAL*8 PROBCHI, T, XP, A0, A1, B1, B2 REAL*8 OCHISQ, ORIGCHISQ, DELCHISQ, NCHISQ, CONF, LOOPFACT REAL*8 EXPCHISQ REAL*8 TORB, FREE REAL*8 A(7), ASAVE(7), DELA(7), DYDA(7,2) REAL*8 LOWVALUE, HIGHVALUE REAL*8 COVAR(7,7), ALPHA(7,7), UNCERT(7,-1:1) REAL*8 DA_TRUE, DE_TRUE, DPERIOD, DEPOCH, DINCLIN REAL*8 DASC_NODE, DARG_PERI INTEGER*4 NDATA, I, J, K, KK, LISTA(7), ITST, NCA, MFIT INTEGER*4 LISTB(7), EFIT, INITMRQ, LOOPCT LOGICAL*1 TOOBIG PARAMETER (NCA = 7) DATA A0,A1,B1,B2/2.30753,0.27061,0.99229,0.04481/ J = 1 DO I = 1, NPTS IF(.NOT.DELETE(I)) J = J + 1 END DO NDATA = 2 * (NPTS-DEL_COUNT) EFIT = MFIT - 1 C GET THE INCREASE IN CHISQR DESIRED ORIGCHISQ = CHISQ CALL LOCATE (0,17) CALL WRITE_STRING('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 CALL LOCATE(0,17) CALL WRITE_STRING('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 = DBLE(NDATA-MFIT) NCHISQ = (1.-(2./(9.*FREE))+XP*DSQRT(2./(9.*FREE)))**3 NCHISQ = NCHISQ*(NDATA-MFIT) 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 EXPCHISQ = EXPCHISQ*FREE DELCHISQ = NCHISQ-EXPCHISQ WRITE(OUTC,'(A,I3,A,F6.3,A,F6.3,A,F4.1,A)') 1 ' Expect Chi^2 (',INT(FREE),' d.o.f) = ',EXPCHISQ,' (min) or ', 2 NCHISQ,' (',CONF*100.,'% confidence).' C IMPLEMENT THIS "IF" STATEMENT IF CHISQ/FREE < 1 BY A SIGNIFICANT C AMOUNT BECAUSE OF SYSTEMATICS (AS OPPOSED TO BECAUSE OF MIS-SCALING C OF UNCERTAINTIES OF DATA POINTS): C IF (CHISQ.GT.EXPCHISQ) THEN NCHISQ = NCHISQ*CHISQ/EXPCHISQ C END IF WRITE(OUTC,'(A,F6.3,A,F6.3,A,F4.1,A)') 1 ' Observe Chi^2(min) = ',CHISQ,'; seek Chi^2 = ',NCHISQ, 2 ' at ',CONF*100.,'% confidence.' C FILL THE ARRAY OF BEST-FIT ORBITAL ELEMENTS DO I = 1,MA ASAVE(I) = A(I) END DO C FILL THE ARRAY OF AMOUNTS TO VARY THE PARAMETERS BY DO I = 1,MA IF(COVAR(I,I).NE.0.) THEN DELA(I) = SQRT(DELCHISQ*COVAR(I,I)) UNCERT(I,0) = SQRT(DELCHISQ*COVAR(I,I)) D WRITE(OUTC,'(A,I1,A,G11.4,A)') D 1 ' UNCERT(',I,') = ',UNCERT(I,0),' from COVAR.' ELSE DELA(I) = 0.0 END IF END DO C READ THE LIST OF FITTED PARAMETERS LISTA() INTO A WORKING LIST C LISTB(). FOR EACH I, FIT THE OTHER MFIT-1 PARAMETERS C WRITE(PR,'(1X,A,7I3)') 'ERREST1: LISTA= ',(LISTA(J),J=1,7) CALL LOCATE(0,13) CALL WRITE_STRING('FRACTIONAL CHANGE IN CHISQ FOR CONVERGENCE: ') READ(INC,*) DELCHISQ DO I = 1,MFIT DO J = 1,MA A(J) = ASAVE(J) LISTB(J) = 0 END DO DO J = 1,MFIT LISTB(J) = LISTA(J) END DO C WRITE(PR,'(1X,A,7I3)') 'ERREST2: LISTB= ',(LISTB(J),J=1,7) LISTB(I) = 0 C MAKE A LIST OF THE MFIT-1 PARAMETERS TO FIT DO J = 1,MFIT IF (LISTB(J).EQ.0.AND.J.LT.MA) THEN LISTB(J) = LISTB(J+1) LISTB(J+1)= 0 END IF END DO C WRITE(PR,'(1X,A,7I3)') 'ERREST3: LISTB= ',(LISTB(J),J=1,7) C VARY PARAMETER I BY +DELA (WHEN J=1) OR BY -DELA (WHEN J=-1) DO J=1,-1,-2 A(LISTA(I)) = ASAVE(LISTA(I)) LOWVALUE = A(LISTA(I)) TOOBIG = .FALSE. A(LISTA(I)) = A(LISTA(I)) + J*DELA(LISTA(I)) D WRITE(PR,'(1X,A,I1,A,I1,A,F20.10,F20.10)') D 1 'ERREST: A(',LISTA(I),'), DELA(',LISTA(I),')=', D 2 A(LISTA(I)),DELA(LISTA(I)) C GET BEST FIT OF OTHER MFIT-1 PARAMETERS IF(MFIT.GT.1) THEN IF(I.EQ.1.AND.J.EQ.1) THEN CALL LOCATE(1,14) CALL WRITE_STRING 1 ('NEWLINENEWLINENEWLINENEWLINENEWLINENEWL ') CALL LOCATE(1,14) CALL WRITE_STRING 1 ('NEWLINENEWLINENEWLINENEWLINENEWLINENEWL ') CALL LOCATE(1,14) CALL WRITE_STRING 1 ('INIT. VAL. OF MARQUARDT COEFF. ALAMDA: ') READ(INC,*) ALAMDA CALL LOCATE(0,15) CALL WRITE_STRING 1 ('NEWLINENEWLINENEWLINENEWLINENEWLINENEWLINENEW ') CALL LOCATE(0,15) CALL WRITE_STRING 1 ('NEWLINENEWLINENEWLINENEWLINENEWLINENEWLINENEW ') CALL LOCATE(0,15) CALL WRITE_STRING 1 ('NUMBER TO MULT. ALAMDA BY AS MARQUARDT METHOD') 11 CALL LOCATE(0,16) CALL WRITE_STRING 1 ('NEWLINENEWLINENEWLINENEWLINENEWLINENEWL ') CALL LOCATE(0,16) CALL WRITE_STRING 1 ('NEWLINENEWLINENEWLINENEWLINENEWLINENEWL ') CALL LOCATE(0,16) CALL WRITE_STRING 1 ('GOES OVER TO PARABOLIC CHISQ FIT: [<1] ') READ(INC,*) LOOPFACT IF (LOOPFACT.GE.1.) GO TO 11 END IF 10 INITMRQ = -1 ! RETURN HERE FOR REPEAT ! OR FOR NEXT STEP C WRITE(OUTC,*) ' Going to MRQMIN' CALL MRQMIN(NDATA,A,LISTB,EFIT,COVAR,ALPHA, 1 NCA,ALAMDA,LOOPFACT,INITMRQ) C WRITE(OUTC,*) ' Back from MRQMIN' K = 1 ITST = 0 100 K = K+1 ! LOOP TO HERE IF NO CONVERGENCE OCHISQ = CHISQ KK = MOD(K,9) C CALL MRQMIN AND OBSERVE CHANGE IN CHISQ C WRITE(OUTC,*) ' Going to MRQMIN' CALL MRQMIN(NDATA,A,LISTB,EFIT,COVAR,ALPHA, 1 NCA,ALAMDA,LOOPFACT,INITMRQ) C WRITE(OUTC,*) ' Back from MRQMIN' IF(CHISQ.LT.OCHISQ) THEN C WRITE(PR,1001) A(1),A(2),A(3),A(4),A(5)/RPDEG,A(6)/RPDEG, C 1 A(7)/RPDEG,CHISQ OCHISQ = CHISQ ELSE IF (CHISQ.GT.OCHISQ) THEN ITST = 0 ELSE IF (ABS(OCHISQ/CHISQ-1.).LT.DELCHISQ) THEN ITST = ITST + 1 END IF END IF IF (ITST.LT.2) GO TO 100 C CALL LOCATE(0,17+KK) C CALL WRITE_STRING (' REPEAT? [N]') C READ(INC,*) COMND C CALL EGA_RESTORE_DEFAULT C IF (COMND.EQ.'Y'.OR.COMND.EQ.'y') GO TO 10 ELSE CHISQ = 0. DO K=1,NPTS CALL TRUE2APP(TOBS(K),A,RMOD,THMOD,RAMOD,DECMOD,DYDA,PR) R_EST(K) = RMOD THETA_EST(K)= THMOD END DO CALL CHISQR(R_EST,ROBS,THETA_EST,THETAOBS,DRA,DDEC, 1 PHI_ERR,NPTS,DELETE,CHISQ,DELRA,DELDEC) C WRITE(PR,1001) A(1),A(2),A(3),A(4),A(5)/RPDEG,A(6)/RPDEG, C 1 A(7)/RPDEG,CHISQ END IF C TEST FOR CONVERGENCE TO DESIRED CHISQR VALUE C --Chi squared still too small IF ((NCHISQ-CHISQ)/NCHISQ.GT.DELCHISQ 1 .AND. .NOT.TOOBIG) THEN D WRITE(OUTC,'(A,I3,A,G11.4,A,G11.4,A,G11.4)') D 1 ' CHISQ TOO SMALL; A(',LISTA(I),')=',A(LISTA(I)), D 2 ' DELA=',DELA(LISTA(I)),' CHISQ=',CHISQ LOWVALUE = A(LISTA(I)) DELA(LISTA(I))= 2.*DELA(LISTA(I)) A(LISTA(I)) = A(LISTA(I)) + J*DELA(LISTA(I))/2. GO TO 10 C --Chi squared got too big ELSE IF ((NCHISQ-CHISQ)/NCHISQ.LT.-DELCHISQ 1 .AND. .NOT.TOOBIG) THEN D WRITE(OUTC,'(A,I3,A,G11.4,A,G11.4,A,G11.4)') D 1 ' CHISQ TOO BIG; A(',LISTA(I),')=',A(LISTA(I)), D 2 ' DELA=',DELA(LISTA(I)),' CHISQ=',CHISQ HIGHVALUE = A(LISTA(I)) A(LISTA(I)) = (LOWVALUE+HIGHVALUE)/2. TOOBIG = .TRUE. GO TO 10 C --Chi squared too small ELSE IF ((NCHISQ-CHISQ)/NCHISQ.GT.DELCHISQ 1 .AND. TOOBIG) THEN D WRITE(OUTC,'(A,I3,A,G11.4,A,G11.4)') D 1 ' CHISQ TOO SMALL; A(',LISTA(I),')=',A(LISTA(I)), D 2 ' CHISQ=',CHISQ LOWVALUE = A(LISTA(I)) A(LISTA(I)) = (LOWVALUE+HIGHVALUE)/2. IF(ABS(HIGHVALUE-LOWVALUE).GT..01*UNCERT(LISTA(I),0)) 1 THEN GO TO 10 END IF C --Chi squared got too big ELSE IF ((NCHISQ-CHISQ)/NCHISQ.LT.-DELCHISQ 1 .AND. TOOBIG) THEN D WRITE(OUTC,'(A,I3,A,G11.4,A,G11.4)') D 1 ' CHISQ TOO BIG; A(',LISTA(I),')=',A(LISTA(I)), D 2 ' CHISQ=',CHISQ HIGHVALUE = A(LISTA(I)) A(LISTA(I)) = (LOWVALUE+HIGHVALUE)/2. IF(ABS(HIGHVALUE-LOWVALUE).GT..01*UNCERT(LISTA(I),0)) 1 THEN GO TO 10 END IF END IF C --Chi squared is juuuuuuuust right UNCERT(LISTA(I),J)= A(LISTA(I)) - ASAVE(LISTA(I)) DELA(LISTA(I)) = A(LISTA(I)) - ASAVE(LISTA(I)) WRITE(OUTC,'(1X,A,I1,A,F20.10)') 1 'DONE WITH THIS VARIABLE; DEL A(',LISTA(I),')= ', 2 A(LISTA(I)) - ASAVE(LISTA(I)) WRITE(PR,1001) A(1),A(2),A(3),A(4),A(5)/RPDEG,A(6)/RPDEG, 1 A(7)/RPDEG,CHISQ END DO IF (LISTA(I).LE.4) THEN WRITE(PR,1000) LISTA(I), ASAVE(LISTA(I)), UNCERT(LISTA(I),1), 1 UNCERT(LISTA(I),-1) CALL LOCATE(0,2+I) WRITE(INS,1000)LISTA(I), ASAVE(LISTA(I)), UNCERT(LISTA(I),1), 1 UNCERT(LISTA(I),-1) CALL WRITE_STRING ( INS ) ELSE WRITE(PR,1000) LISTA(I), ASAVE(LISTA(I))/RPDEG, 1 UNCERT(LISTA(I),1)/RPDEG, 2 UNCERT(LISTA(I),-1)/RPDEG CALL LOCATE(0,2+I) WRITE(INS,1000)LISTA(I), ASAVE(LISTA(I))/RPDEG, 1 UNCERT(LISTA(I),1)/RPDEG, 2 UNCERT(LISTA(I),-1)/RPDEG CALL WRITE_STRING ( INS ) END IF END DO CALL PGASK(.TRUE.) CALL PGADVANCE CALL PGASK(.FALSE.) RETURN 1000 FORMAT(1X,'ERREST: A(',I1,')= ',1X,F20.10,2X,F20.10,2X,F20.10) 1001 FORMAT(1X,'ERREST: ELEMENTS = ',F7.3,1X,F6.4,1X,F9.4,1X,F12.4,1X, 1 3(F8.3,1X),' CHISQ = ',F8.2) END C===================================================================== SUBROUTINE GAMMQ(GAMMQ,A,X) C FROM NUMERICAL RECIPES, P. 161 REAL*8 GAMMQ,A,X,GAMSER,GAMMCF,GLN IF(X.LT.0..OR.A.LE.0.)PAUSE IF(X.LT.A+1.)THEN CALL GSER(GAMSER,A,X,GLN) GAMMQ=1.-GAMSER ELSE CALL GCF(GAMMCF,A,X,GLN) GAMMQ=GAMMCF ENDIF RETURN END C===================================================================== SUBROUTINE GSER(GAMSER,A,X,GLN) C FROM NUMERICAL RECIPES, P. 162 INTEGER*4 ITMAX, N REAL*8 GAMSER,A,X,GLN,AP,SUM,DEL,EPS PARAMETER (ITMAX=100,EPS=3.D-7) CALL GAMMLN(GLN,A) IF(X.LE.0.)THEN IF(X.LT.0.)PAUSE GAMSER=0. RETURN ENDIF AP=A SUM=1./A DEL=SUM DO 11 N=1,ITMAX AP=AP+1. DEL=DEL*X/AP SUM=SUM+DEL IF(ABS(DEL).LT.ABS(SUM)*EPS)GO TO 1 11 CONTINUE PAUSE 'A too large, ITMAX too small' 1 GAMSER=SUM*DEXP(-X+A*DLOG(X)-GLN) RETURN END C===================================================================== SUBROUTINE GCF(GAMMCF,A,X,GLN) C FROM NUMERICAL RECIPES, P. 162 INTEGER*4 ITMAX,N REAL*8 GAMMCF,A,X,GLN,EPS,GOLD,A0,A1,B0,B1,FAC,AN,ANA,ANF,G PARAMETER (ITMAX=100,EPS=3.E-7) CALL GAMMLN(GLN,A) GOLD=0. A0=1. A1=X B0=0. B1=1. FAC=1. DO 11 N=1,ITMAX AN=FLOAT(N) ANA=AN-A A0=(A1+A0*ANA)*FAC B0=(B1+B0*ANA)*FAC ANF=AN*FAC A1=X*A0+ANF*A1 B1=X*B0+ANF*B1 IF(A1.NE.0.)THEN FAC=1./A1 G=B1*FAC IF(ABS((G-GOLD)/G).LT.EPS)GO TO 1 GOLD=G ENDIF 11 CONTINUE PAUSE 'A too large, ITMAX too small' 1 GAMMCF=DEXP(-X+A*DLOG(X)-GLN)*G RETURN END C===================================================================== SUBROUTINE GAMMLN(GAMLN,XX) C FROM NUMERICAL RECIPES, P. 157 REAL*8 COF(6),STP,HALF,ONE,FPF,X,TMP,SER,XX,HALF,ONE,FPF,GAMLN DATA COF,STP/76.18009173D0,-86.50532033D0,24.01409822D0, * -1.231739516D0,.120858003D-2,-.536382D-5,2.50662827465D0/ DATA HALF,ONE,FPF/0.5D0,1.0D0,5.5D0/ X=XX-ONE TMP=X+FPF TMP=(X+HALF)*DLOG(TMP)-TMP SER=ONE DO 11 J=1,6 X=X+ONE SER=SER+COF(J)/X 11 CONTINUE GAMLN=TMP+DLOG(STP*SER) RETURN END