# B4J Code SnippetGamma and Beta functions for computing Cumulative Distribution Functions

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

B4X:
``````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``````

Replies
4
Views
700
Replies
4
Views
2K
Replies
11
Views
7K
Replies
37
Views
10K
Replies
109
Views
26K