Emulate BETADIST function

mark14 · Dec 22, 2010

Hi everyone,

I have another application that I'm developing in Oracle Application Express using Oracle 10g. I developed a 'prototype' of the Oracle application in Excel and used the BETADIST function in this prototype. However, Oracle doesn't have this function, and I basically need to re-write it.

I've read all the documentation on the BETADIST function and how to use it in Excel, so I have a pretty good handle on that. I understand that you provide your x, alpha, beta, and optionally, lower and upper bounds, and the function returns the cumulative beta probability density function.

For example, =BETADIST(0.2,1.57,2) gives me 0.1803.

The problem is, I'm no statistician, and have no idea what it is doing behind the scenes to get the 0.1803 value in order to recreate it. I've checked the Wikipedia pages, among others, and all of those formulas are greek to me.

Can anyone help break down this function for me and explain how that value is logically derived in my example above, using those numbers? That will greatly help me figure out how to rewrite it.

Thanks in advance!
Mark

shg · Dec 22, 2010

This is hacked together from bits and pieces at hand, and minimally tested:

Code:

Option Explicit
 
Const n             As Long = 200    ' increase for accuracy, decrease for speed
 
Public aa           As Double
Public bb           As Double
 
Function BetaDist1(x As Double, a As Double, b As Double)
    Dim d1          As Double
    Dim d2          As Double
    Dim n1          As Long
    Dim n2          As Long
 
    aa = a
    bb = b
    n1 = x * n
    n2 = n - n1
 
    d1 = SimpsonInt(0, x, n1)
    d2 = SimpsonInt(x, 1, n2)
    BetaDist1 = d1 / (d1 + d2)
End Function
 
Function SimpsonInt(ti As Double, tf As Double, ByVal n As Long) As Double
    ' shg 2006
 
    ' Returns the integral of Func (below) from ti to tf _
      using Composite Simpson's Rule over n intervals
    Dim i           As Double  ' index
    Dim dH          As Double  ' step size
    Dim dOdd        As Double  ' sum of Func(i), i = 1, 3, 5, 7, ... n-1, i.e., n/2 values
    Dim dEvn        As Double  ' sum of Func(i), i =   2, 4, 6,  ... n-2  i.e., n/2 - 1 values
    ' 1 + (n/2) + (n/2 - 1) + 1 = n+1 function evaluations
 
    If n < 1 Then Exit Function
 
    If n And 1 Then n = n + 1    ' n must be even
    dH = (tf - ti) / n
 
    For i = 1 To n - 1 Step 2
        dOdd = dOdd + Func(ti + i * dH)
    Next i
 
    For i = 2 To n - 2 Step 2
        dEvn = dEvn + Func(ti + i * dH)
    Next i
 
    SimpsonInt = (Func(ti) + 4# * dOdd + 2# * dEvn + Func(tf)) * dH / 3#    ' weighted sum
End Function
 
Function Func(t As Double) As Double
    Func = t ^ (aa - 1) * (1 - t) ^ (bb - 1)
End Function

It gives reasonable accuracy, which could be improved at the expense of speed by increasing the constant n:

Code:

      -------A-------- ---B--- ---C--- ---D--- ---E--- ---F---
  1                  x     0.1     0.3     0.5     0.7     0.9
  2                  a    1.57    1.57    1.57    1.57    1.57
  3                  b       2       2       2       2       2
  4                                                           
  5    BETADIST(x,a,b) 0.06495 0.31702 0.60120 0.84027 0.98060
  6   BetaDist1(x,a,b) 0.06489 0.31698 0.60118 0.84026 0.98060

mark14 · Dec 22, 2010

Thanks shg4421! I'll see if I digest that logic and implement into my app.

shg · Dec 23, 2010

Here's a more complete emulation of BETADIST.

Code:

Option Explicit
 
Const nInt          As Long = 500    ' increase for accuracy, decrease for speed
 
'======================= B e t a D i s t  E x a m p l e  =======================
' shg 2010-1223
 
' Beta distribution parameters
Public mAlpha       As Double
Public mBeta        As Double
 
Private Function Func(x As Double) As Double
    ' The evaluated function for BetaDist
    Func = x ^ (mAlpha - 1) * (1 - x) ^ (mBeta - 1)
End Function
 
' The UDF/VBA function itself
Function BetaDist1(ByVal x As Double, _
                   alpha As Double, beta As Double, _
                   Optional A As Double = 0#, Optional B As Double = 1#) As Variant
 
    Dim n1          As Long     ' intervals in left-side integration
    Dim n2          As Long     ' intervals in right-side integration
    Dim d1          As Double   ' result of left-side integration
    Dim d2          As Double   ' result of right-side integration
 
    If alpha <= 0# Or beta <= 0# Or x < A Or x > B Or A = B Then
        BetaDist1 = CVErr(xlErrNum)
 
    Else
        x = (x - A) / (B - A)   ' scale to interval
 
        mAlpha = alpha
        mBeta = beta
 
        n1 = x * nInt           ' divide iterations between left and right sides
        n2 = nInt - n1
 
        d1 = SimpsonInt(0#, x, n1)  ' get left-side integral
        d2 = SimpsonInt(x, 1#, n2)  ' get right-side integral
        BetaDist1 = d1 / (d1 + d2)  ' compute ratio of left side to total
    End If
End Function
 
'============= G e n e r a l   S i m p s o n   I n t e g r a t o r =============
Private Function SimpsonInt(xi As Double, xf As Double, ByVal n As Long) As Double
    ' shg 2006
 
    ' Returns the integral of Func (specific to problem) over the interval xi to xf
    ' using the Composite Simpson's Rule over n intervals
    Dim i           As Double  ' index
    Dim dH          As Double  ' step size
    Dim dOdd        As Double  ' sum of Func(i), i = 1, 3, 5, 7, ..., n-1 (n/2 values)
    Dim dEvn        As Double  ' sum of Func(i), i =   2, 4, 6, ..., n-2  (n/2 - 1 values)
    ' That's a total of  1 + (n/2) + (n/2 - 1) + 1 = n + 1 function evaluations
 
    If n < 1 Then Exit Function
 
    If n And 1 Then n = n + 1    ' n must be even
    dH = (xf - xi) / n
 
    For i = 1 To n - 2 Step 2
        dOdd = dOdd + Func(xi + i * dH)
        dEvn = dEvn + Func(xi + (i + 1) * dH)
    Next i
    dOdd = dOdd + Func(xi + i * dH)
 
    SimpsonInt = (Func(xi) + 4 * dOdd + 2 * dEvn + Func(xf)) * dH / 3    ' weighted sum
End Function

mark14 · Dec 27, 2010

Thanks again shg4421. I'm not sure what the # symbols are for. Can you help explain what they are doing, like in the following lines...

Code:

If alpha <= 0# Or beta <= 0#

Code:

d1 = SimpsonInt(0#, x, n1)

Thanks!

shg · Dec 27, 2010

Type declaration character for a Double.

mark14 · Dec 28, 2010

shg4421, I've been able to convert your VBA code into Oracle PL/SQL, and after some testing and tweaking, it works like a charm!

Thanks so much for your help with this question. There's no telling how long it would've taken me to figure this out if not for you.

Thanks again!
Mark

shg · Dec 28, 2010

Good job, mark, you're welcome.

mark14 · Feb 20, 2011

Hello shg4421 (and others),

I've run into a problem with this BetaDist1 emulation of BETADIST, where the function gives me the #VALUE! error for some particular parameter combinations.

Here's an example of one that causes the error:
BETADIST(0.3, 401, 1001, 0, 1). These parameters are valid and the BETADIST function returns 0.876099. However, when using the same parameters in BetaDist1, it errors out.

I stepped through the code, and it appears the main cause of the error is a result of attempting to divide by zero in the line below when both d1 and d2 end up equalling zero, as is the case with this example.

Code:

d1 = SimpsonInt(0#, x, n1)  ' get left-side integral
d2 = SimpsonInt(x, 1#, n2)  ' get right-side integral
BetaDist1 = d1 / (d1 + d2)  ' compute ratio of left side to total

However, that's not really the root cause. I assume that the root cause is that, for some reason, d1 and d2 are equalling zero when they shouldn't. Any idea how we can fix this so that BetaDist1 gives me the same answer as BETADIST in cases such as these? If BETADIST can do it, I really need BetaDist1 to do it.

Thanks for your help!
Mark

shg · Feb 23, 2011

The problem is with A, B, and the function that computed the PDF for integration:

Code:

Func = x ^ (mAlpha - 1) * (1 - x) ^ (mBeta - 1)

With alpha = 400 and beta = 1000, you're raising numbers < 1 to huge powers that go beyond (normalized) floating point numbers ability to represent. For example, in Excel,

0.49^1000 ~ 9E-302

0.5^1000 = 0

There's obviously a way to rearrange the calculation to circumvent this (the native function manages), but I'm not seeing it at the moment, mark.

Emulate BETADIST function

mark14

New Member

Excel Facts

shg

MrExcel MVP

mark14

New Member

shg

MrExcel MVP

mark14

New Member

shg

MrExcel MVP

mark14

New Member

shg

MrExcel MVP

mark14

New Member

shg

MrExcel MVP

Similar threads

Forum statistics

Share this page

Emulate BETADIST function

New Member

Excel Facts

MrExcel MVP

New Member

MrExcel MVP

New Member

MrExcel MVP

New Member

MrExcel MVP

New Member

MrExcel MVP

Similar threads

Forum statistics

Share this page

We've detected that you are using an adblocker.

Which adblocker are you using?

Disable AdBlock

Disable AdBlock Plus

Disable uBlock Origin

Disable uBlock