' ******************************************************************
' Probability functions related to the incomplete Beta function
' ******************************************************************

' ------------------------------------------------------------------
' Error codes
' ------------------------------------------------------------------

#DEFINE FOk 0
' No error

#DEFINE FDomain -1
' Argument domain error

' ------------------------------------------------------------------
' Global variable
' ------------------------------------------------------------------

COMMON SHARED ErrCode AS INTEGER
' Error code from the latest function evaluation

' ------------------------------------------------------------------
' External function
' ------------------------------------------------------------------

DECLARE FUNCTION IBeta(A AS DOUBLE,_ 
                       B AS DOUBLE,_ 
                       X AS DOUBLE) AS DOUBLE
' Incomplete Beta function

' ******************************************************************

FUNCTION FBeta(A AS DOUBLE, B AS DOUBLE, X AS DOUBLE) AS DOUBLE
' Returns Prob(Z <= X), where Z follows the Beta distribution
' with parameters A and B

  FBeta = IBeta(A, B, X)
END FUNCTION

FUNCTION FBinom(N AS INTEGER, P AS DOUBLE, K AS INTEGER) AS DOUBLE
' Returns Prob(X <= K) where X follows the 
' binomial distribution B(N,P)

  IF P < 0 OR P > 1 OR N <= 0 OR N < K THEN
    ErrCode = FDomain
    FBinom = 0
  ELSEIF K = 0 THEN
    ErrCode = FOk
    FBinom = (1 - P) ^ N
  ELSEIF K = N THEN
    ErrCode = FOk
    FBinom = 1
  ELSE
    FBinom = 1 - IBeta(K + 1, N - K, P)
  END IF
END FUNCTION

FUNCTION FSnedecor(Nu1 AS INTEGER, Nu2 AS INTEGER,_
                   X AS DOUBLE) AS DOUBLE
' Returns Prob(F <= X), where F follows the Fisher-Snedecor
' distribution with Nu1 and Nu2 degrees of freedom

  IF Nu1 < 1 OR Nu2 < 1 OR X <= 0 THEN
    ErrCode = FDomain
    FSnedecor = 0
  ELSE
    FSnedecor = 1 - IBeta(0.5 * Nu2, 0.5 * Nu1, Nu2 / (Nu2 + Nu1 * X))
  END IF
END FUNCTION

FUNCTION PSnedecor(Nu1 AS INTEGER, Nu2 AS INTEGER,_
                   X AS DOUBLE) AS DOUBLE
' Returns Prob(F > X), where F follows the Fisher-Snedecor
' distribution with Nu1 and Nu2 degrees of freedom

  IF Nu1 < 1 OR Nu2 < 1 OR X <= 0 THEN
    ErrCode = FDomain
    PSnedecor = 0
  ELSE
    PSnedecor = IBeta(0.5 * Nu2, 0.5 * Nu1, Nu2 / (Nu2 + Nu1 * X))
  END IF
END FUNCTION

FUNCTION FStudent(Nu AS INTEGER, X AS DOUBLE) AS DOUBLE
' Returns Prob(t <= X), where t follows the Student distribution
' with Nu degrees of freedom

  DIM F AS DOUBLE

  IF Nu < 1 THEN
    ErrCode = FDomain
    FStudent = 0
  ELSEIF X = 0 THEN
    ErrCode = FOk
    FStudent = 0.5
  ELSE
    F = 0.5 * IBeta(0.5 * Nu, 0.5, Nu / (Nu + X * X))
    IF X < 0 THEN FStudent = F ELSE FStudent = 1 - F
  END IF
END FUNCTION

FUNCTION PStudent(Nu AS INTEGER, X AS DOUBLE) AS DOUBLE
' Returns Prob(|t| > X), where t follows the Student distribution
' with Nu degrees of freedom

  IF Nu < 1 THEN
    ErrCode = FDomain
    PStudent = 0
  ELSE
    PStudent = IBeta(0.5 * Nu, 0.5, Nu / (Nu + X * X))
  END IF
END FUNCTION

