      program T1mprobadr
      parameter (n=3)

      real data1(n),data2(n)

      data data1/290,311,284/
      data data2/271,304,260/
      call tptest(data1,data2,n,t,prob)

      write(*,*)t,prob

      stop
      end

      SUBROUTINE tptest(data1,data2,n,t,prob)
      INTEGER n
      REAL prob,t,data1(n),data2(n)
C     USES avevar,betai
C     Given the paired arrays data1(1:n) and data2(1:n),
C     this routine returns Student's t for
C     paired data as t, and its signicance as prob, small
C     values of prob indicating a signicant
C     difference of means.
      INTEGER j
      REAL ave1,ave2,cov,df,sd,var1,var2,betai
      call avevar(data1,n,ave1,var1)
      call avevar(data2,n,ave2,var2)
      cov=0.
      do j=1,n
      cov=cov+(data1(j)-ave1)*(data2(j)-ave2)
      enddo
      df=n-1
      cov=cov/df
      sd=sqrt((var1+var2-2.*cov)/n)
      t=(ave1-ave2)/sd
      prob=betai(0.5*df,0.5,df/(df+t**2))
      return
      END

      SUBROUTINE avevar(data,n,ave,var)
      INTEGER n
      REAL ave,var,data(n)
C      Given array data(1:n), returns its mean as ave and its variance
C     as var.
      INTEGER j
      REAL s,ep
      ave=0.0
      do j=1,n
      ave=ave+data(j)
      enddo
      ave=ave/n
      var=0.0
      ep=0.0
      do j=1,n
      s=data(j)-ave
      ep=ep+s
      var=var+s*s
      enddo
      var=(var-ep**2/n)/(n-1)
      return
      END

      FUNCTION betai(a,b,x)
      REAL betai,a,b,x
C     USES betacf,gammln
C     Returns the incomplete beta function Ix(a,b).
      REAL bt,betacf,gammln
      if(x.lt.0..or.x.gt.1.)pause
C      'bad argument x in betai'
      if(x.lt.0..or.x.eq.1.) then
	    bt=0.
      else
C     Factors in front of the continued fraction.
	    bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.-x))
      endif
      if(x.lt.(a+1.)/(a+b+2.))then
C     Use continued fraction directly.
	    betai=bt*betacf(a,b,x)/a
	    return
      else
      betai=1.-bt*betacf(b,a,1.-x)/b
C     Use continued fraction after making the symmetry transformation.
      return
      endif
      END


      FUNCTION betacf(a,b,x)
      INTEGER MAXIT
      REAL betacf,a,b,x,EPS,FPMIN
      PARAMETER (MAXIT=100,EPS=3.e-7,FPMIN=1.e-30)
C     Used by betai: Evaluates continued fraction for
C     incomplete beta function by modied
C     Lentz's method ( x 5.2).
      INTEGER m,m2
      REAL aa,c,d,del,h,qab,qam,qap
      qab=a+b
C     These q's will be used in factors that occur in the coecients
C     (6.4.6).
      qap=a+1.
      qam=a-1.
      c=1.
C     First step of Lentz's method.
      d=1.-qab*x/qap
      if(abs(d).lt.FPMIN)d=FPMIN
      d=1./d
      h=d
      do m=1,MAXIT
      m2=2*m
      aa=m*(b-m)*x/((qam+m2)*(a+m2))
      d=1.+aa*d
C     One step (the even one) of the recurrence.
      if(abs(d).lt.FPMIN)d=FPMIN
      c=1.+aa/c
      if(abs(c).lt.FPMIN)c=FPMIN
      d=1./d
      h=h*d*c
      aa=-(a+m)*(qab+m)*x/((a+m2)*(qap+m2))
      d=1.+aa*d
C     Next step of the recurrence (the odd one).
      if(abs(d).lt.FPMIN)d=FPMIN
      c=1.+aa/c
      if(abs(c).lt.FPMIN)c=FPMIN
      d=1./d
      del=d*c
      h=h*del
      if(abs(del-1.).lt.EPS)goto 1
C     Are we done?
      enddo
      pause
C     'a or b too big, or MAXIT too small in betacf'
    1 betacf=h
      return
      END

      FUNCTION gammln(xx)
      REAL gammln,xx
C     Returns the value ln((xx)) for xx > 0.
      INTEGER j
      DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
C     Internal arithmetic will be done in double precision, a nicety
C     that you can omit if five-figure accuracy is good enough.
      SAVE cof,stp
      DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,
     X24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
     X-.5395239384953d-5,2.5066282746310005d0/
      x=xx
      y=x
      tmp=x+5.5d0
      tmp=(x+0.5d0)*log(tmp)-tmp
      ser=1.000000000190015d0
      do j=1,6
	    y=y+1.d0
	    ser=ser+cof(j)/y
      enddo
      gammln=tmp+log(stp*ser/x)
      return
      END

