SUBROUTINE tutest(data1,n1,data2,n2,t,prob) INTEGER n1,n2 REAL prob,t,data1(n1),data2(n2) CU USES avevar,betai REAL ave1,ave2,df,var1,var2,betai call avevar(data1,n1,ave1,var1) call avevar(data2,n2,ave2,var2) t=(ave1-ave2)/sqrt(var1/n1+var2/n2) df=(var1/n1+var2/n2)**2/((var1/n1)**2/(n1-1)+(var2/n2)**2/(n2-1)) prob=betai(0.5*df,0.5,df/(df+t**2)) return END C (C) Copr. 1986-92 Numerical Recipes Software "W#(1&. c************************************************ C (C) Copr. 1986-92 Numerical Recipes Software "W#(1&. FUNCTION BETAI(A,B,X) IF(X.LT.0..OR.X.GT.1.)PAUSE 'bad argument X in BETAI' IF(X.EQ.0..OR.X.EQ.1.)THEN BT=0. ELSE BT=EXP(GAMMLN(A+B)-GAMMLN(A)-GAMMLN(B) + +A*ALOG(X)+B*ALOG(1.-X)) ENDIF IF(X.LT.(A+1.)/(A+B+2.))THEN BETAI=BT*BETACF(A,B,X)/A RETURN ELSE BETAI=1.-BT*BETACF(B,A,1.-X)/B RETURN ENDIF print*,'betai ',a,b,x END FUNCTION GAMMLN(XX) REAL*8 COF(6),STP,HALF,ONE,FPF,X,TMP,SER 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)*LOG(TMP)-TMP SER=ONE DO 11 J=1,6 X=X+ONE SER=SER+COF(J)/X 11 CONTINUE GAMMLN=TMP+LOG(STP*SER) RETURN END FUNCTION BETACF(A,B,X) 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 11 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))GO TO 1 11 CONTINUE PAUSE 'A or B too big, or ITMAX too small' 1 BETACF=AZ RETURN END