'*************************************************************
'*                STUDENT  T-PROBABILITY  LAW                *
'* --------------------------------------------------------- *
'* SAMPLE RUN:                                               *
'* (Calculate Student T-probability (lower tail and upper    *
'*  tail) for T=0.257).                                      *
'*                                                           *
'*  X= .257                                                  *
'*  PROB1= .6000293520025725                                 *
'*  ERROR= 0                                                 *
'*  X= .257                                                  *
'*  PROB2= .3999704738211087                                 *
'*  ERROR= 0                                                 *
'*  PROB1+PROB2= .9999998258236812                           *
'*                                                           *
'* --------------------------------------------------------- *
'* Ref.:"JOURNAL OF APPLIED STATISTICS (1968) VOL.17, P.189, *
'*       & VOL.19, NO.1".                                    *
'*                                                           *
'*                      Basic Release By J-P Moreau, Paris.  *
'*                               (www.jpmoreau.fr)           *
'*************************************************************
'PROGRAM STUDENT
DEFDBL A-H, O-Z
DEFINT I-N

  ZERO = 0#
  HALF = .5#
  ONE = 1#
  TWO = 2#
  FOUR = 4#
  PI = 3.1415926535#

' X,XNU,PROB1,PROB2: REAL
' NU, IFAULT: INTEGER

  X = .257#
  NU = 19

  GOSUB 1000 'call PROBST(X,NU,PROB1,IFAULT)
  CLS
  PRINT
  PRINT "  X="; X
  PRINT "  PROB1="; PROB1
  PRINT "  ERROR="; IFAULT

  X = .257#
  XNU = 19#

  GOSUB 2000 'call STUDNT(X,XNU,PROB2,IFAULT)
  PRINT "  X="; T
  PRINT "  PROB2="; PROB2
  PRINT "  ERROR="; IFAULT

  PRINT "  PROB1+PROB2="; PROB1 + PROB2
  PRINT

END 'of main program


1000 'SUBROUTINE PROBST(T, IDF, PROB1, IFAULT)
' ---------------------------------------------------------------------
'        ALGORITHM AS 3  APPL. STATIST. (1968) VOL.17, P.189
'
'        STUDENT T PROBABILITY (LOWER TAIL)
' ---------------------------------------------------------------------
       G10 = ONE / PI

'      A, B, C, F, S, FK: REAL
'      IM2,IOE,K,KS:Integer

      T = X: IDF = NU
      IFAULT = 1
      PROB1 = ZERO
      IF IDF < 1 THEN RETURN
      IFAULT = 0
      F = ONE * IDF
      A = T / SQR(F)
      B = F / (F + T * T)
      IM2 = IDF - 2
      IOE = IDF MOD 2
      S = ONE
      C = ONE
      F = ONE
      KS = 2 + IOE
      FK = KS
      IF IM2 < 2 THEN GOTO 20
      K = KS
      WHILE K <= IM2
	C = C * B * (FK - ONE) / FK
	S = S + C
	IF S = F THEN GOTO 20
	F = S
	FK = FK + TWO
	K = K + 2
      WEND
20    IF IOE = 1 THEN GOTO 30
      PROB1 = HALF + HALF * A * SQR(B) * S
      RETURN
30    IF IDF = 1 THEN S = ZERO
      PROB1 = HALF + (A * B * S + ATN(A)) * G10
RETURN


2000  'SUBROUTINE STUDNT (T, DOFF, PROB2, IFAULT)
' ----------------------------------------------------------------
'     ALGORITHM AS 27  APPL. STATIST. VOL.19, NO.1
'
'     Calculate the upper tail area under Student's t-distribution
'
'     Translated from Algol by Alan Miller
' ----------------------------------------------------------------
'     V, X, TT: REAL
'     A1, A2, A3, A4, A5, B1, B2, C1, C2, C3, C4, C5, D1, D2: REAL
'     E1, E2, E3, E4, E5, F1, F2, G1, G2, G3, G4, G5, H1, H2: REAL
'     XI1, XI2, XI3, XI4, XI5, XJ1, XJ2: REAL
'     IPOS: Integer (0 or 1)

      T = X: DOFF = XNU

      A1 = .09979441#: A2 = -.581821: A3 = 1.390993: A4 = -1.222452: A5 = 2.151185
      B1 = 5.537409: B2 = 11.42343
      C1 = .04431742#: C2 = -.2206018#: C3 = -.03317253#: C4 = 5.679969: C5 = -12.96519
      D1 = 5.166733: D2 = 13.49862
      E1 = 9.694901000000001D-03: E2 = -.1408854#: E3 = 1.88993: E4 = -12.75532: E5 = 25.77532
      F1 = 4.233736: F2 = 14.3963
      G1 = -9.187228E-05: G2 = .03789901#: G3 = -1.280346: G4 = 9.249528: G5 = -19.08115
      H1 = 2.777816: H2 = 16.46132
      XI1 = 5.79602E-04: XI2 = -.02763334#: XI3 = .4517029#: XI4 = -2.657697: XI5 = 5.127212
      XJ1 = .5657187#: XJ2 = 21.83269

'     Check that number of degrees of freedom > 4

      IF DOFF < TWO THEN
	IFAULT = 1
	PROB2 = -ONE
	RETURN
      END IF

      IF DOFF <= FOUR THEN
	IFAULT = INT(DOFF)
      ELSE
	IFAULT = 0
      END IF

'     Evaluate series

      V = ONE / DOFF
      IF T >= ZERO THEN
	IPOS = 1
      ELSE
	IPOS = 0
      END IF
      TT = ABS(T)
'     X =  ONE + TT * (((A1 + V * (A2 + V * (A3 + V * (A4 + V * A5)))) / (ONE - V * (B1 - V * B2))) +
'          TT * (((C1 + V * (C2 + V * (C3 + V * (C4 + V * C5)))) / (ONE - V * (D1 - V * D2))) +
'          TT * (((E1 + V * (E2 + V * (E3 + V * (E4 + V * E5)))) / (ONE - V * (F1 - V * F2))) +
'          TT * (((G1 + V * (G2 + V * (G3 + V * (G4 + V * G5)))) / (ONE - V * (H1 - V * H2))) +
'          TT * ((XI1 + V * (XI2 + V * (XI3 + V * (XI4 + V * XI5)))) / (ONE - V * (XJ1 - V * XJ2)))))))
      TMP1 = ((A1 + V * (A2 + V * (A3 + V * (A4 + V * A5)))) / (ONE - V * (B1 - V * B2)))
      TMP2 = ((C1 + V * (C2 + V * (C3 + V * (C4 + V * C5)))) / (ONE - V * (D1 - V * D2)))
      TMP3 = ((E1 + V * (E2 + V * (E3 + V * (E4 + V * E5)))) / (ONE - V * (F1 - V * F2)))
      TMP4 = ((G1 + V * (G2 + V * (G3 + V * (G4 + V * G5)))) / (ONE - V * (H1 - V * H2)))
      TMP5 = (XI1 + V * (XI2 + V * (XI3 + V * (XI4 + V * XI5))))
      TMP6 = (ONE - V * (XJ1 - V * XJ2))
      X = ONE + TT * (TMP1 + TT * (TMP2 + TT * (TMP3 + TT * (TMP4 + TT * (TMP5 / TMP6)))))

      X = HALF * (X ^ (-8))
      IF IPOS = 1 THEN
	PROB2 = X
      ELSE
	PROB2 = ONE - X
      END IF

RETURN

'end of file Student.bas