B4J Code Snippet Gamma and Beta functions for computing Cumulative Distribution Functions

Sub AppStart (Args() As String)
    Dim markTime As Long = DateTime.Now
    Log("2-tailed z = 3.291   p=" & (NumberFormat2((1-cdfZ(3.291)), 1,4,4,False)))
    Log("right-tail chiSq(14) = 36.1   p=" & NumberFormat2((1-cdfC(36.1, 14)),1,4,4,False))
    Log("right-tail F(3,5) = 33.2025   p=" & NumberFormat2((1-cdfF(33.2025, 3, 5)), 1,4,4,False))
    Log("2-tailed Pearson R(14) = 0.742    p=" & NumberFormat2((1-cdfR(0.742, 14)), 1,4,4,False))
    Log("2-tailed t-Value(13) = 4.221   p=" & NumberFormat2((1-cdfT(4.221, 13)), 1,4,4,False))
    Log((DateTime.now - markTime))   '~15 milliseconds in release mode
End Sub

Sub cdfZ(x As Double) As Double
    If x>=0 Then Return (gammaNormal(.5, x * x / 2)) Else Return (1 - cdfZ(-x))
End Sub

Sub cdfC(x As Double, v As Double) As Double
    If x>=0 Then Return (gammaNormal(.5 * v, x / 2)) Else Return (1 - cdfC(-x, v))
End Sub

Sub cdfT(x As Double, v As Double) As Double
    If x>=0 Then Return (1 - betaNormal(v / 2, .5, v / (v + x*x))) Else Return (1 - cdfT(-x, v))
End Sub

Sub cdfF(x As Double, v1 As Double, v2 As Double) As Double
    If x>=0 Then Return (1 - betaNormal(v2 / 2, v1 / 2, v2 / (v2 + v1*x))) Else Return (1 - cdfF(-x, v1, v2))
End Sub

Sub cdfR(x As Double, v As Double) As Double
    If x>=0 Then Return cdfT(x * Sqrt(v)/Sqrt(1.0 - x*x), v) Else Return cdfR(-x, v)
End Sub

Sub gammaNormal(a As Double, x As Double) As Double
    Dim sum, delta, aStep, result As Double
    result = 0.0
    If (x = 0.0) Then
        result = .5
    Else
        aStep = a: delta = 1.0 / a
        sum = delta
        Dim n As Int = 0
        Do While (Abs(delta) > Abs(sum) * .0000001)
            n = n + 1
            If n>100 Then Return 0.0
            aStep = aStep + 1
            delta = delta * x / aStep
            sum = sum + delta
        Loop
        result = sum * Power(cE, -x + a * Logarithm(x,cE) - logGamma(a))
    End If
    Return result
End Sub

Sub betaNormal(a As Double, b As Double, x As Double) As Double
    If (x < 0.0 Or x > 1.0) Then Return 0.0
    Dim result As Double
    If (x = 0.0 Or x = 1.0) Then
        result = 0.0
    Else
        result = Power(cE, logGamma(a + b) - logGamma(a) - logGamma(b) + a * Logarithm(x, cE) + b * Logarithm(1.0 - x, cE))
    End If
    If (x < (a + 1.0)/(a + b + 2.0)) Then
        result = result * continuousFractions(a, b, x) / a
    Else
        result =  1.0 - result * continuousFractions(b, a, 1.0-x) / b
    End If
    Return result
End Sub

Sub continuousFractions(a As Double, b As Double, x As Double) As Double
    Dim m, m2 As Int
    Dim aa, c, d, delta, result As Double
    c = 1.0
    d = 1.0 / Abs(1.0 - (a +  b) * x / (a + 1))
    result = d
    m = 0: m2 = 0
    Do While Abs(delta-1.0) > .0000001
        m = m + 1
        If m>100 Then Return 0.0
        m2 = m2 + 2
        aa = m * (b - m) * x / ((a + m2) * (a - 1 + m2))
        c = (1.0 + aa/c)
        d = 1.0 / (1.0 + aa*d)
        result = result * c * d
    'alternating two different fractional terms
        aa = -(a + m) * (a + b + m) * x / ((a + m2) * (a + 1 + m2))
        c = (1.0 + aa/c)
        d = 1.0 / (1.0 + aa*d)
        delta = c * d
        result = result * delta
    Loop
    Return result
End Sub

Sub logGamma(z As Double) As Double
    Dim  zStep, denominator, series, Lanczos() As Double   'Lanzos coefficients g=5 n=7
    Lanczos = Array As Double( 1.000000000190015, _
        76.18009172947146, -86.50532032941677, 24.01409824083091, _
        -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5)
    denominator = (z + 5.5) - (z + 0.5) * Logarithm((z + 5.5), cE)
    series = Lanczos(0)
    zStep = z
    For j = 1 To 6
        zStep = zStep + 1
        series = series + Lanczos(j) / zStep
    Next
    Return Logarithm(Sqrt(2*cPI) * series / z, cE) - denominator
End Sub