Skip to content
logo-mike
Menu
  • Home
  • ASP NET Core
  • ASP NET 5
  • ASP NET 7
  • Javascript
  • C#
  • Other
Menu
Beta Function

The LogBeta and LogGamma Functions Using C#

Posted on August 10, 2022

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.

READ  Learn How to Use Structures in C#

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.

Recent Posts

  • Microsoft 365 Basic – What’s in the new Office suite? Is it worth buying?
  • Office 365 vs Office 2021 – on-premises or cloud – which should you choose?
  • Microsoft Office 365 – List of releases, most important changes for the company
  • The office will become Microsoft 365
  • How to Check If the Docker Daemon or a Container Is Running
  • .NET 7 WebAssembly Plans: Mixed-Mode AOT, Multi-Threading, Web Crypto
  • Opdex Launches Stratis AMM DEX to Support GameFi Ecosystem
  • Minimal Bitcoin Miner in C#
  • Privacy Policy
  • Terms and Conditions
  • Cookies Policy
  • Contact Us