' ******************************************************************
' Probability functions related to the incomplete Gamma function
' and Error function
' ******************************************************************

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

#DEFINE FOk 0
' No error

#DEFINE FDomain -1
' Argument domain error

#DEFINE FUnderflow -4
' Underflow range error

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

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

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

#DEFINE Sqrt2Div2 0.7071067811865476#
' Sqrt(2)/2

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

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

' ------------------------------------------------------------------
' External functions
' ------------------------------------------------------------------

DECLARE FUNCTION IGamma(A AS DOUBLE, X AS DOUBLE) AS DOUBLE
' Incomplete Gamma function

DECLARE FUNCTION Erf(X AS DOUBLE) AS DOUBLE
' Error function

DECLARE FUNCTION Erfc(X AS DOUBLE) AS DOUBLE
' Complement of error function

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

FUNCTION FPoisson(Mu AS DOUBLE, K AS INTEGER) AS DOUBLE
' Returns Prob(X <= K) for the Poisson distribution P(Mu)

  IF Mu <= 0 OR K < 0 THEN
    ErrCode = FDomain
    FPoisson = 0
  ELSEIF K = 0 THEN
    IF (-Mu) < MinLog THEN
      ErrCode = FUnderflow
      FPoisson = 0
    ELSE
      ErrCode = FOk
      FPoisson = EXP(-Mu)
    END IF
  ELSE
    FPoisson = 1 - IGamma(K + 1, Mu)
  END IF
END FUNCTION

FUNCTION FNorm(X AS DOUBLE) AS DOUBLE
' Returns Prob(U <= X), where U follows the standard normal
' distribution

  FNorm = 0.5 * (1 + Erf(X * Sqrt2div2))
END FUNCTION

FUNCTION PNorm(X AS DOUBLE) AS DOUBLE
' Returns Prob(|U| > X), where U follows the standard normal
' distribution

  DIM A AS DOUBLE

  A = ABS(X)

  IF A = 0 THEN
    ErrCode = FOk
    PNorm = 1
  ELSEIF A < 1 THEN
    PNorm = 1 - Erf(A * Sqrt2div2)
  ELSE
    PNorm = Erfc(A * Sqrt2div2)
  END IF
END FUNCTION

FUNCTION FKhi2(Nu AS INTEGER, X AS DOUBLE) AS DOUBLE
' Returns Prob(Khi2 <= X) for the Khi-2 distribution
' with Nu degrees of freedom

  IF Nu < 1 OR X <= 0 THEN
    ErrCode = FDomain
    FKhi2 = 0
  ELSE
    FKhi2 = IGamma(0.5 * Nu, 0.5 * X)
  END IF
END FUNCTION

FUNCTION PKhi2(Nu AS INTEGER, X AS DOUBLE) AS DOUBLE
' Returns Prob(Khi2 > X) for the Khi-2 distribution
' with Nu degrees of freedom

  IF Nu < 1 OR X <= 0 THEN
    ErrCode = FDomain
    PKhi2 = 0
  ELSE
    PKhi2 = 1 - IGamma(0.5 * Nu, 0.5 * X)
  END IF
END FUNCTION

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

  FGamma = IGamma(A, B * X)
END FUNCTION
