C#

Function RMSBetaDist(x, alpha, beta, rate)
     'Main function where the cumulative Beta distribution is assembled

datasize = Range("Datsize")
Dim bt() As Double
Dim RMSBeta() As Double
Dim i As Integer

ReDim bt(datasize)
ReDim RMSBeta(datasize)

For i = 1 To datasize
If x(i, 1) < 0 Or x(i, 1) > 1 Then
    'Error = MsgBox("There is a problem with the x value of the beta function", vbOKOnly)
    'GoTo Abend
    x(i, 1) = 1
End If
If x(i, 1) = 0 Then
    bt(i) = 0
    RMSBeta(i) = 0
    GoTo Nextarray
ElseIf x(i, 1) = 1 Then
    bt(i) = 0
    RMSBeta(i) = 1
    GoTo Nextarray
Else
    bt(i) = Exp(RMSGammaln(alpha(i, 1) + beta(i, 1)) - RMSGammaln(alpha(i, 1)) - RMSGammaln(beta(i, 1)) + alpha(i, 1) * Log(x(i, 1)) + beta(i, 1) * Log(1 - x(i, 1)))
End If
If x(i, 1) < (alpha(i, 1) + 1) / (alpha(i, 1) + beta(i, 1) + 2) Then
    RMSBeta(i) = bt(i) * RMSBetaCF(x(i, 1), alpha(i, 1), beta(i, 1)) / alpha(i, 1)
Else                           'symmetry relation I(a,b)=1-I(b,a) where I() is the incomplete beta function
    RMSBeta(i) = 1 - bt(i) * RMSBetaCF(1 - x(i, 1), beta(i, 1), alpha(i, 1)) / beta(i, 1)
End If

Nextarray:
RMSBetaDist = RMSBetaDist + rate(i) * (1 - RMSBeta(i))
Next i

Abend:
End Function

Function RMSGammaln(xg)   'Lanczos gamma approximation
     'Used in the Beta distribution approximation

Dim y, temp, serial, PI, small As Double
Dim coeff(6) As Double
Dim j As Integer

PI = 3.14159265358979
coeff(0) = 76.18009173
coeff(1) = -86.50532033
coeff(2) = 24.01409822
coeff(3) = -1.231739516
coeff(4) = 0.00120858003
coeff(5) = -0.00000536382
small = 0
If xg < 1 Then
    small = xg
    y = xg - 1
Else
    y = xg - 1
End If
temp = y + 5.5
temp = temp - (y + 0.5) * Log(temp)
serial = 1

For j = 0 To 5
    y = y + 1
    serial = serial + coeff(j) / y
Next j
If small <> 0 Then   'formula for values of x < 1
    RMSGammaln = -temp + Log(Sqr(2 * PI) * serial)
Else
    RMSGammaln = -temp + Log(Sqr(2 * PI) * serial)
End If

End Function

Function RMSBetaCF(x2, alpha2, beta2)
    'Used to approximate the integral portion of the cumulative Beta distribution

Dim errorf As Double
Dim qap, qam, qab, em, temp, d As Double
Dim bz, bm, bp, bpp As Double
Dim az, am, ap, app, aold As Double
Dim m, MaxIterations As Integer

errorf = 0.0000001
MaxIterations = 1000
bm = 1
az = 1
am = 1
qab = alpha2 + beta2
qap = alpha2 + 1
qam = alpha2 - 1
bz = 1 - qab * x2 / qap

For m = 1 To MaxIterations
    em = m
    temp = em + em
    d = em * (beta2 - em) * x2 / ((qam + temp) * (alpha2 + temp))
    ap = az + d * am
    bp = bz + d * bm
    d = -(alpha2 + em) * (qab + em) * x2 / ((qap + temp) * (alpha2 + temp))
    app = ap + d * az
    bpp = bp + d * bz
    aold = az
    am = ap / bpp
    bm = bp / bpp
    az = app / bpp
    bz = 1
    If Abs(az - aold) < errorf * Abs(az) Then
        GoTo FinalStep
    End If
Next m
ErrorStep:
    'Error = MsgBox("Either alpha or beta is too big or the maximum iterations are too small for convergence", vbOKOnly, "Problem in RMSBetaCF function")
    RMSBetaCF = 1
    GoTo Abend
FinalStep:
    RMSBetaCF = az
Abend:
    'Msg = MsgBox(m & " iterations ", vbOKOnly)
End Function