' ******************************************************************
' Probability functions related to the Gamma function
' ******************************************************************

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

#DEFINE FOk 0
' No error

#DEFINE FDomain -1
' Argument domain error

#DEFINE FSing -2
' Singularity

#DEFINE FUnderflow -4
' Underflow range error

' ------------------------------------------------------------------
' Machine-dependent constant
' ------------------------------------------------------------------

#DEFINE MaxNum 1.797693134862315D+308   
' Max. floating point number: 2^1024

#DEFINE MinLog -708.3964185322641#       
' Min. argument for Exp

' ------------------------------------------------------------------
' Mathematical constant
' ------------------------------------------------------------------

#DEFINE Pi 3.141592653589793#
' Pi

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

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

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

DECLARE FUNCTION LnGamma (X AS DOUBLE) AS DOUBLE
' Ln of Gamma function

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

FUNCTION DBeta(A AS DOUBLE, B AS DOUBLE, X AS DOUBLE) AS DOUBLE
' Density of Beta distribution with parameters A and B

  DIM AS DOUBLE L

  IF A <= 0 OR B <= 0 OR X < 0 OR X > 1 THEN
    ErrCode = FDomain
    DBeta = 0
    EXIT FUNCTION
  END IF

  IF X = 0 THEN
    IF A < 1 THEN
      ErrCode = FSing
      DBeta = MaxNum
    ELSE
      ErrCode = FOk
      DBeta = 0
    END IF
    EXIT FUNCTION
  END IF

  IF X = 1 THEN
    IF B < 1 THEN
      ErrCode = FSing
      DBeta = MaxNum
    ELSE
      ErrCode = FOk
      DBeta = 0
    END IF
    EXIT FUNCTION
  END IF

  L = LnGamma(A + B) - LnGamma(A) - LnGamma(B)_
      + (A - 1) * LOG(X) + (B - 1) * LOG(1 - X)

  IF L < MinLog THEN
    ErrCode = FUnderflow
    DBeta = 0
  ELSE
    ErrCode = FOk
    DBeta = EXP(L)
  END IF
END FUNCTION

FUNCTION DGamma(A AS DOUBLE, B AS DOUBLE, X AS DOUBLE) AS DOUBLE
' Density of Gamma distribution with parameters A and B

  DIM AS DOUBLE L

  IF A <= 0 OR B <= 0 OR X < 0 THEN
    ErrCode = FDomain
    DGamma = 0
    EXIT FUNCTION
  END IF

  IF X = 0 THEN
    IF A < 1 THEN
      ErrCode = FSing
      DGamma = MaxNum
    ELSEIF A = 1 THEN
      ErrCode = FOk
      DGamma = B
    ELSE
      ErrCode = FOk
      DGamma = 0
    END IF
    EXIT FUNCTION
  END IF

  L = A * LOG(B) - LnGamma(A) + (A - 1) * LOG(X) - B * X

  IF L < MinLog THEN
    ErrCode = FUnderflow
    DGamma = 0
  ELSE
    ErrCode = FOk
    DGamma = EXP(L)
  END IF
END FUNCTION

FUNCTION DKhi2(Nu AS INTEGER, X AS DOUBLE) AS DOUBLE
' Density of Khi-2 distribution with Nu d.o.f.

  DKhi2 = DGamma(0.5 * Nu, 0.5, X)
END FUNCTION

FUNCTION DSnedecor(Nu1 AS INTEGER, Nu2 AS INTEGER,_
                   X AS DOUBLE) AS DOUBLE
' Density of Fisher-Snedecor distribution with Nu1 and Nu2 d.o.f.

  DIM AS DOUBLE L, P1, P2, R, S

  IF Nu1 < 1 OR Nu2 < 1 OR X <= 0 THEN
    ErrCode = FDomain
    DSnedecor = 0
    EXIT FUNCTION
  END IF

  R = Nu1 / Nu2
  P1 = 0.5 * Nu1
  P2 = 0.5 * Nu2
  S = P1 + P2
  L = LnGamma(S) - LnGamma(P1) - LnGamma(P2)_
      + P1 * LOG(R) + (P1 - 1) * LOG(X) - S * LOG(1 + R * X)

  IF L < MinLog THEN
    ErrCode = FUnderflow
    DSnedecor = 0
  ELSE
    ErrCode = FOk
    DSnedecor = EXP(L)
  END IF
END FUNCTION

FUNCTION DStudent(Nu AS INTEGER, X AS DOUBLE) AS DOUBLE
' Density of Student distribution with Nu d.o.f.

  DIM AS DOUBLE L, P, Q

  IF Nu < 1 THEN
    ErrCode = FDomain
    DStudent = 0
    EXIT FUNCTION
  END IF

  P = 0.5 * (Nu + 1)
  Q = 0.5 * Nu
  L = LnGamma(P) - LnGamma(Q) - 0.5 * LOG(Nu * Pi)_
      - P * LOG(1 + X * X / Nu)

  IF L < MinLog THEN
    ErrCode = FUnderflow
    DStudent = 0
  ELSE
    ErrCode = FOk
    DStudent = EXP(L)
  END IF
END FUNCTION

