The .NET library (formerly called .NET Core) does not have any built-in functions for classic statistical evaluations. But it is possible to implement such functionality from scratch. This post shows C# execution of three of the most important classic data functions: log-beta, log-gamma, and regulated incomplete beta.
The regulated insufficient beta function is usually formed as Ix(a,b) or I(x;a,b). It is used in a variety of ways, including statistical assessment of variability (ANOVA). The function Ix(a,b) calls the log beta function. The log beta function calls the log gamma function.
The first part of the result shows the value of I(x;a,b) for x=0.5, a=4, 0, b = 7.0, so 0.828125. The regulated incomplete beta function has no obvious analysis. It’s best to think of this as an abstract mathematical function that validates an x-value between 0.0 and 1.0 and a and b values that are positive. The return value is a number between 0.0 and 1.0, so this function is said to be “regulated”.
The second part of the result shows the LogBet value (3.0, 5.0) which is -4.653960. Also this function has no visible analysis and is also the best idea as it is simply a function that accepts 2 favorable values. There is a simple beta(a,b) function, but the return value from beta() can be large, so the default way is to use LogBeta(), which is essentially a log of the beta() function.
The third part of the result shows the LogGamma(4.0) value, which is 6.0000. The regular Gamma() function is a generalization of the Factorial() function. If n is an integer, then Gamma(n) = Factorial(n-1). Nonetheless, Gamma(z) is defined for the preferred real Z. Because the return value of Gamma(z) can be astronomically large, this is the primary method of using LogGamma(z).
This short article assumes you have moderate to much better C# programming skills, but don’t assume you know anything about beta, gamma, or regulated under-beta features. The full test code is presented in this post and is also available in the attached download. All regular bug checking has been duly eliminated to keep the main ideas as clear as possible.
To create the test program, I used Aesthetic Workshop 2022 (free community edition) with .NET 5 (to avoid the annoying new .NET 6 console theme). Nevertheless, the test version has no significant dependencies, so that all reasonably current VS and .NET versions can be used.
LogGamma() function
Notification that Gamma(n) = Factorial(n-1), for example Gamma(4.0) = Factorial(3) = 3 * 2 * 1 = 6.0.
LogGamma(z) cannot be calculated exactly when z is not an integer, so there are many estimation formulas. The demo program uses the so-called Lanczos formula g = 5, n = 7. The submission has been received.
static double LogGamma(double z)
{
// Lanczos approximation g=5, n=7
double[] coef = new double[7] { 1.000000000190015,
76.18009172947146, -86.50532032941677,
24.01409824083091, -1.231739572450155,
0.1208650973866179e-2, -0.5395239384953e-5 };
double LogSqrtTwoPi = 0.91893853320467274178;
if (z < 0.5)
return Math.Log(Math.PI / Math.Sin(Math.PI * z)) -
LogGamma(1.0 - z);
double zz = z - 1.0;
double b = zz + 5.5; // g + 0.5
double sum = coef[0];
for (int i = 1; i < coef.Length; ++i)
sum += coef[i] / (zz + i);
return (LogSqrtTwoPi + Math.Log(sum) - b) +
(Math.Log(b) * (zz + 0.5));
}
The demo program calls LogGamma() with the following declarations:
double z = 4.0;
double lg = LogGamma(z);
Console.WriteLine("LogGamma(4.0) = ");
Console.WriteLine(lg.ToString("F6"));
double g = Math.Exp(lg);
Console.WriteLine("Gamma(4.0) = ");
Console.WriteLine(g.ToString("F6"));
Computing Gamma(4.0) by calling Math.Exp() in LogGamma(4.0) is done simply to show that LogGamma() is just a log of a gamma function.
The process does not perform error monitoring. In a non-demo situation, you probably want to check the value of the input specification to make sure it’s positive.
LogBeta() function
And (x;a,b) the regulated incomplete beta function calls LogBeta(), which then calls LogGamma(). The test application of the LogBeta() function is:
static double LogBeta(double a, double b)
{
double logbeta = LogGamma(a) + LogGamma(b) -
LogGamma(a + b);
return logbeta;
}
The test program calls the LogBeta() function as follows:
double lb = LogBeta(3.0, 5.0);
Console.WriteLine("LogBeta(3.0, 5.0) = ");
Console.WriteLine(lb.ToString("F6"));
Mathematically, Beta(a,b) = (Gamma(a) * Gamma(b)) / Gamma(a + b). Taking log() on both sides of the equation and using the facts that log(x/y) = log(x) – log(y) and also log(x * y) = log(x) + log(y ) get the result in LogBeta() function application.
Regulated incomplete beta feature.
The regulated insufficient beta function I(x;a,b) is extremely difficult to perform. The demo program uses the so-called Estimated Amount of Amounts Produced (CF).
The approximation I(x;a,b) CF is the mathematical position of two fundamental terms. The term on the left includes beta (a, b). The term on the right has an infinite variety of terms d[i], where the value of each term depends on whether [i] is even or odd.
static double RegIncBeta(double x, double a, double b)
{
// pick the form of RegIncompleteBeta() that converges best
if (x < (a + 1.0) / (a + b + 2.0))
return RegIncompleteBeta(x, a, b);
else
return 1.0 - RegIncompleteBeta(1 - x, b, a);
}
static double RegIncompleteBeta(double x, double a, double b)
{
// regularized incomplete beta
// calls LogBeta() and helper ContFraction()
double cf = ContFraction(x, a, b);
double logTop = (a * Math.Log(x)) +
(b * Math.Log(1 - x));
double logBot = Math.Log(a) + LogBeta(a, b);
double logLeft = logTop - logBot;
double left = Math.Exp(logLeft);
return left * cf; // should be [0.0, 1.0]
}
static double ContFraction(double x, double a, double b,
int maxTerms = 100)
{
// 1. pre-compute 100 d values
double[] d = new double[maxTerms]; // d[0] not used
int end = maxTerms / 2; // 50
for (int m = 0; m < end; ++m) // [0,49]
{
int i = 2 * m; // even di
int j = i + 1; // odd di
d[i] = (m * (b - m) * x) / ((a + 2 * m - 1) *
(a + 2 * m));
d[j] = -1 * ((a + m) * (a + b + m) * x) / ((a + 2 * m) *
(a + 2 * m + 1));
}
// 2. work backwards
double[] t = new double[maxTerms]; // t[0] not used
t[maxTerms - 1] = 1 + d[maxTerms - 1];
for (int j = maxTerms - 2; j >= 1; --j)
t[j] = 1 + d[j] / t[j + 1];
// ShowVector(t);
// Console.ReadLine();
double cf = 1 / t[1];
return cf;
}
There is a top-level RegIncBeta() function. Calls the RegIncompleteBeta() function, which does all the work. The wrapper technique is based on the math of home ownership that I(x;a,b) = 1 – I(x-1;b,a) and calls the version that converges faster. The RegIncompleteBeta() function calls the helper function ContFraction(), which evaluates the term CF. To avoid beta(a,b) calculations that can easily overflow, the RegIncompleteBeta() function evaluates the LogBeta(a,b) and then undoes the logging operation by calling Math.Exp().
The demo program calls the function I(x;a,b) with these statements:
double x = 0.5;
double a = 4.0;
double b = 7.0;
double rib = RegIncBeta(x, a, b);
Console.WriteLine("Ix(a,b) for (x=0.5, a=4.0, b=7.0) = ");
Console.WriteLine(rib.ToString("F6"));
The helper function ContFraction() initially calculates 100 d[i] values and then calculates the term CF from the last d[i] value to the original d[i] value. This technique can fail on some combinations of x, an, and also b with underflow or numeric overflow, so in a non-demo scenario you probably intend to wrap your code in a try-catch block.
Unlike using 100 precomputed terms d[i] for the helper function ContFraction(), you can look at the value of the function RegIncBeta() starting with maxTerms = 2 and also by bringing maxTerms to RegIncBeta() the function converges to a safe value. For example with x = 0.5, a = 4.0, b = 7.0, RegIncBeta() handles with maxTerms = 10:.
maxTerms RegIncBeta()
=======================
2 0.812500
4 0.828613
6 0.828119
8 0.828125
10 0.828125
Another approach is to use the so-called Lentz formula to estimate the value of the CF component. Instead of precomputing an arbitrary multitude of terms d[i] and then working in reverse order to the last value of d[i], the Lentz formula uses some brilliant mathematical methods to proceed by checking on each iteration that the CF value was merged with constant value. However, the Lentz formula is usually very brittle and also does not converge for some problems.
Completing
There are many external code libraries that contain applications for complex statistical functions such as beta, log-beta, gamma, log-gamma, and adjusted incomplete beta. For example, the scipy Python language collection includes these functions. You can use such external collections in many scenarios. However, applying statistics functions from scratch with C# is used when you want to avoid external dependency for technological or copyright reasons.