[jira] [Commented] (NUMBERS-181) Update the RegularizedBeta function using the Boost implementation

["Alex Herbert (Jira)" <ji...@apache.org>]

    [ https://issues.apache.org/jira/browse/NUMBERS-181?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=17479015#comment-17479015 ] 

Alex Herbert commented on NUMBERS-181:
--------------------------------------

The ported Boost beta function has been tested using the reference data from the Boost test suite and supplemented with extra data to hit full coverage and add spot checks. The new implementation has been compared to the current implementations on the same data. This was performed to identify the maximum and root mean squared ULP error on each test dataset.
h2. Regularized Beta

The current gamma package implements the regularized beta function using the classic continued fraction representation (see [http://functions.wolfram.com/06.21.10.0001.01]). This requires a leading power term divided by the beta function. Since the current gamma package has no beta function this is currently implemented using the log beta function. It could also be implemented using the log gamma implementation (this is done in commons statistics for the regularized beta derivative).
{code:java}
beta(a, b) = exp(log(gamma(a)) + log(gamma(b)) - log(gamma(a+b)))
           = exp(log(beta(a, b)){code}
The power term is implemented in a separate method in the new Boost version. Thus it is possible to test the continued fraction version with a new updated computation for the leading power terms.
||Dataset||Current|| ||New|| ||New continued fraction|| ||
|| ||max||RMS||max||RMS||max||RMS||
|Small Integer|113.282|16.4048|7.03929|1.10146|243.542|13.6968|
|Small|55.1872|8.42909|12.9271|3.01663|40.8128|4.63655|
|Medium|496.337|64.7578|187.541|24.6896|183.541|23.8939|
|Large|472128|17719.2|2155.54|93.9803|136300|7023.06|

The current version is outperformed by the new version. The new implementation has much better performance on the small and large dataset. There is a lot of side stepping using recursive properties of the beta and gamma functions to move the function into a region that can be computed with high accuracy. It also uses other representations than the classic continued fraction to compute the result depending on the parameters.

Note that the continued fraction version using the new Boost beta functions to compute the power terms performs well on the medium dataset. However this does not handle the other datasets which are special (integer) or have more diverse values for the parameters a and b.
h2. Beta vs LogBeta

The current package has no beta function implementation. I tested the new version against exp(log(beta(a, b))
||Dataset||exp(log(beta))|| ||beta|| ||
|| ||max||RMS||max||RMS||
|Small|103|25.191|3.64656|1.54721|
|Medium|263.239|42.5923|193.292|32.9232|
|Divergent|25.1787|3.94835|15.7892|3.25813|

The new beta function is more accurate than using exp with the log beta. The difference is greatest on the small dataset (a and b < 1e-3) where the expected result of the beta function is always above 1e6. Here the use of exp will convert a small log value to a large value and any errors are increased.
h2. Regularized Beta Derivative

The is no implementation for this in the current gamma package. This function is the probability density function of the beta distribution. A version of the function is implemented using log gamma in the commons statistics. I copied this implementation and create a new version using log beta.
||Dataset||LogGamma|| ||LogBeta|| ||beta derivative|| ||
|| ||max||RMS||max||RMS||max||RMS||
|Small Integer|199.6|48.1154|164.272|28.8628|54.155|12.0911|
|Small|68.0044|9.51327|68.0044|7.67129|19.5118|3.87678|
|Medium|1443.5|286.979|459.872|77.3431|137.08|29.9185|
|Large|8093700|236188|26893.3|1088.98|3673.03|236.86|

The current version in commons statistics should have been implemented using the LogBeta function. It may originate from Commons Math before a LogBeta class existed. However the new beta derivative method is much better than the LogBeta version.
h2. API Additions

The boost function can replace the implementation of RegularizedBeta.value(x, a, b). It has versions that are regularized and the full (non-regularized) incomplete beta; each has a complement function.

There are new functions for:
 * Beta
 * Incomplete beta
 * Complement of the incomplete beta
 * Complement of the regularized incomplete beta
 * Derivative of the regularized incomplete beta

For reference the current beta functions and relevant equivalent gamma functions are implemented in the API as:
{code:java}
LogBeta.value(a, b)
RegularizedBeta.value(x, a, b)

IncompleteGamma.Lower.value(a, x)
IncompleteGamma.Upper.value(a, x) 
RegularizedGamma.P.value(a, x)
RegularizedGamma.P.derivative(a, x)
RegularizedGamma.Q.value(a, x)
RegularizedGamma.Q.derivative(a, x){code}
Unlike the regularized incomplete gamma functions which are commonly known as P and Q there is no common term for the complement of the regularized incomplete beta. It is often just defined using the regularized incomplete beta I{~}x{~}(a, b):
{code:java}
ibetac(x, a, b) = 1 - I[x](a, b) = I[1 - x](b, a){code}
Note that the regularized incomplete beta is typically denoted using a subscript x as I{~}x{~}(a, b) or without as I(x; a, b). Hence the use of x as the first argument is following convention in the current API (and contrasts with the trailing x in the gamma functions).

Boost defines the functions as:
{noformat}
beta(a, b)
// full incomplete uses the 'beta' name
beta(a, b, x)
betac(a, b, x)
// regularized incomplete uses 'ibeta'
ibeta(a, b, x)
ibetac(a, b, x)
ibeta_derivative(a, b, x){noformat}
Matlab implements this function using betainc(x, a, b) with an optional 'lower'/'upper' argument for the complement:
{noformat}
betainc(x, a, b) == betainc(x, a, b, 'lower')
betainc(x, a, b, 'upper'){noformat}
Python's SciPy has betainc but no complement:
{noformat}
scipy.special.betainc(a, b, x){noformat}
The following are possible APIs:
{code:java}
Beta.value(a, b)


IncompleteBeta.Lower.value(x, a, b)
IncompleteBeta.Upper.value(x, a, b)
RegularizedBeta.Lower.value(x, a, b)
RegularizedBeta.Lower.derivative(x, a, b)
RegularizedBeta.Upper.value(x, a, b)
RegularizedBeta.Upper.derivative(x, a, b)

OR

// No inner classes

IncompleteBeta.value(x, a, b) 
IncompleteBeta.complement(x, a, b)
RegularizedBeta.value(x, a, b)       *** current API ***
RegularizedBeta.complement(x, a, b)
RegularizedBeta.derivative(x, a, b)

OR

// Separate class for the complement.
// But where should the derivative be placed?

RegularizedBetac.value(x, a, b)


OR

// Add the incomplete functions to the Beta class

Beta.value(x, a, b)
Beta.complement(x, a, b){code}
I would favour the option with no inner classes. This would add a complement and derivative method to RegularizedBeta and introduce a new IncompleteBeta class and a Beta class. It avoids using the terms upper and lower and uses value and complement as the method names.

 

> Update the RegularizedBeta function using the Boost implementation
> ------------------------------------------------------------------
>
>                 Key: NUMBERS-181
>                 URL: https://issues.apache.org/jira/browse/NUMBERS-181
>             Project: Commons Numbers
>          Issue Type: Improvement
>          Components: gamma
>    Affects Versions: 1.1
>            Reporter: Alex Herbert
>            Priority: Major
>             Fix For: 1.1
>
>
> The gamma functions were updated to use the Boost C++ implementation in NUMBERS-174. This was done for the methods where accuracy was improved over the existing implementation.
> Investigate porting the Boost C++ implementation of the beta function. This may improve accuracy over the existing implementation.
> The implementation may solve the accuracy of the RegularizedBeta when a and b are largeish (NUMBERS-52) and will solve the computation of the complement of the RegularizedBeta function (NUMBERS-169).



--
This message was sent by Atlassian Jira
(v8.20.1#820001)