C @(#)distri.for 17.1.1.1 (ES0-DMD) 01/25/02 17:12:45 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C REAL FUNCTION GAMMAQ(A,X) REAL X,A REAL GLN,GAMSER,GAMMCF IF (X.LT.A+1)THEN CALL GSER(GAMSER,A,X,GLN) GAMMAQ=1-GAMSER ELSE CALL GCF(GAMMCF,A,X,GLN) GAMMAQ=GAMMCF ENDIF RETURN END REAL FUNCTION GAMMAP(A,X) REAL X,A REAL GLN,GAMSER,GAMMCF IF (X.LT.A+.1)THEN CALL GSER(GAMSER,A,X,GLN) GAMMAP=GAMSER ELSE CALL GCF(GAMMCF,A,X,GLN) GAMMAP=1.-GAMMCF ENDIF RETURN END SUBROUTINE GSER(GAMSER,A,X,GLN) REAL GAMSER,X,GLN,SUM,DEL,EPS,A REAL GAMMLN,AP INTEGER ITMAX,N PARAMETER (ITMAX=100,EPS=3.E-7) GLN=GAMMLN(A) IF (X.LE.0)THEN IF(X.LT.0)RETURN GAMSER=0 RETURN ENDIF AP = A SUM = 1./A DEL = SUM DO 10 N=1,ITMAX AP = AP+1 DEL = DEL*X/AP SUM = SUM+DEL IF (ABS(DEL).LT.ABS(SUM)*EPS)GOTO 20 10 CONTINUE 20 GAMSER=SUM*EXP(-X+A*LOG(X)-GLN) RETURN END SUBROUTINE GCF(GAMMCF,A,X,GLN) REAL GAMMCF,X,GLN,GOLD,A0,A1,B0,B1,FAC,EPS,ANA,G,A REAL GAMMLN,ANF INTEGER N,ITMAX PARAMETER (ITMAX=100,EPS=3.E-7) GLN=GAMMLN(A) A0 = 1. A1 = X B0 = 0 B1 = 1 FAC = 1 DO 10 N = 1,ITMAX ANA = N-A A0 = (A1+A0*ANA)*FAC B0 = (B1+B0*ANA)*FAC ANF = N*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)GOTO 20 GOLD=G ENDIF 10 CONTINUE 20 GAMMCF = EXP(-X+A*ALOG(X)-GLN)*G RETURN END REAL FUNCTION GAMMLN(XX) REAL XX INTEGER J DOUBLE PRECISION COF(6),STP,HALF,ONE,FPF,X,TMP,SER DATA COF/76.18009173,-86.50532033,24.01409822, 1 -1.231739516,0.120858003D-2,-0.536382D-5/ DATA STP/2.50662827465/ DATA HALF,ONE,FPF/0.5,1.,5.5/ X = XX - ONE TMP = X + FPF TMP = (X+HALF)*LOG(TMP)-TMP SER = ONE DO 10 J=1,6 X = X + ONE SER = SER + COF(J)/X 10 CONTINUE GAMMLN = TMP + LOG(STP*SER) RETURN END REAL FUNCTION ERF(X) REAL X,GAMMAP c returns the error function erf(X) IF(X.LT.0) THEN ERF = -1*GAMMAP(0.5,X**2) ELSE ERF = GAMMAP(0.5,X**2) ENDIF RETURN END REAL FUNCTION BETAI(A,B,X) C returns the imcomplete beta function REAL*4 A,B,X,GAMMLN,BT,BETACF IF (X.EQ.0.OR.X.EQ.1.)THEN BT = 0. ELSE BT = EXP(GAMMLN(A+B)-GAMMLN(A)-GAMMLN(B)+ 1 A*LOG(X)+B*LOG(1.-X)) ENDIF IF (X.LT.(A+1.)/(A+B+2.))THEN BETAI = BT*BETACF(A,B,X)/A ELSE BETAI = 1.-1*(BT*BETACF(B,A,1.-X)/B) ENDIF END REAL FUNCTION BETACF(A,B,X) REAL A,B,X,TEM,APP,BPP,AOLD REAL EPS,AM,BM,AZ,QAB,QAM,QAP,BZ,EM,D,AP,BP INTEGER ITMAX,M PARAMETER (ITMAX=100,EPS=3.E-7) AM = 1. BM = 1. AZ = 1. QAB = A+B QAP = A+1. QAM = A-1. BZ = 1.-QAB*X/QAP DO 10 M = 1,ITMAX EM = M TEM = EM+EM D = EM*(B-M)*X/((QAM+TEM)*(A+TEM)) AP = AZ+D*AM BP = BZ+D*BM D = -(A+EM)*(QAB+EM)*X/((A+TEM)*(QAP+TEM)) APP = AP+D*AZ BPP = BP+D*BZ AOLD = AZ AM = AP/BPP BM = BP/BPP AZ = APP/BPP BZ = 1. IF(ABS(AZ-AOLD).LT.EPS*ABS(AZ))GOTO 20 10 CONTINUE 20 BETACF = AZ RETURN END REAL FUNCTION ERFCC(X) REAL X,Z,T Z = ABS(X) T = 1./(1.+0.5*Z) ERFCC = T*EXP(-Z*Z-1.26551223+T*(1.00002368+T*(0.37409196+ 1 T*(0.09678418+T*(-0.18628806+T*(0.27886807+T*(-1.13520398+ 2 T*(1.48851587+T*(-0.82215223+T*0.17087277))))))))) IF (X.LT.0.) ERFCC = 2.-ERFCC RETURN END REAL FUNCTION UNIF(X,PARA) REAL X,PARA(2),A,B A = PARA(1) B = PARA(2) IF (X.LT.B.AND.X.GT.A)THEN UNIF = (B-X)/(B-A) ELSE UNIF = 0 ENDIF END REAL FUNCTION BINOM(NOBS,N,PROB) REAL FACT INTEGER NOBS,N,I REAL PROB BINOM = 0. DO 10 I = 0,NOBS BINOM = BINOM+FACT(N)/FACT(I)/FACT(N-I)*(PROB**I)* 1 ((1-PROB)**(N-I)) 10 CONTINUE RETURN END REAL FUNCTION FACT(N) REAL GAMMLN INTEGER N,J FACT = 1. IF (N.LT.0) THEN RETURN ELSE IF (N.LT.32) THEN DO 11 J=1,N FACT = FACT*J 11 CONTINUE ELSE FACT = EXP(GAMMLN(N+1.)) ENDIF RETURN END