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

Discussion in 'B4J Code Snippets' started by William Lancee, Dec 3, 2018.

  1. William Lancee

    William Lancee Member Licensed User

    There are numerous implementations of Gamma and Beta functions in many computer languages
    This is my implementation for B4X as applied to statistics.

    http://en.wikipedia.org/wiki/Gamma_function
    http://www.stat.tamu.edu/~jnewton/604/chap3.pdf
    http://www.aip.de/groups/soe/local/numres/bookfpdf/f6-1.pdf

    Code:
    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.114)),1,4,4,False))
        
    Log("right-tail F(3,5) = 33.2025   p=" & NumberFormat2((1-cdfF(33.202535)), 1,4,4,False))
        
    Log("2-tailed Pearson R(14) = 0.742    p=" & NumberFormat2((1-cdfR(0.74214)), 1,4,4,False))
        
    Log("2-tailed t-Value(13) = 4.221   p=" & NumberFormat2((1-cdfT(4.22113)), 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.0Then
            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.0Then Return 0.0
        
    Dim result As Double
        
    If (x = 0.0 Or x = 1.0Then
            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.5053203294167724.01409824083091, _
            -
    1.2317395724501550.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
     
Loading...
  1. This site uses cookies to help personalise content, tailor your experience and to keep you logged in if you register.
    By continuing to use this site, you are consenting to our use of cookies.
    Dismiss Notice