Computing the Regularized Incomplete Beta Function From Scratch Using C#

The regularized incomplete beta function is very strange. One form of the function looks like I(x=0.5, a=5.0, b=3.0) = 0.2265. Another form is Ix(a,b). The x value is between 0 and 1, and a and b are positive. The result is a value between 0.0 and 1.0.

The Ix(a,b) function has no obvious intuitive meaning and isn’t used directly. But Ix(a,b) is used as a helper function in many numerical functions. For example, you can use Ix(a,b) to find the area under the F-distribution for an ANOVA problem, where a and b are derived from df values and x is derived from the computed F statistic.

One weekend I sat down and decided to implement Ix(a,b) from scratch using the C# language. This was primarily for mental exercise and entertainment. (Yes, my idea of entertainment is probably like yours if you’re reading this blog post).


The continued fractions definition of Ix(a,b)

Anyway, after about six hours I had a working demo up and running.

It was crazy complicated both conceptually and from an engineering perspective. It would take many pages to explain fully, so instead, here are a few key ideas of my from-scratch implementation:

* I used the continued fractions definition (eq. 26.5.8) of Ix(a,b) from the A_and_S Handbook.
* I computed the continued fraction forward with a default of 200 terms rather than use Lenz’s algorithm.
* I used Log(Gamma(a,b)) to compute Beta(a,b) to avoid overflow.
* I implemented a LogGamma() function using the Lanczos g=5, n=7 approximation.
* I didn’t add any significant error-checking: Ix(a,b) can blow up in many ways.

I tested my from-scratch version by comparing its results with results from the Python language betainc() function in the SciPy library.

Even though the Ix(a,b) function is extraordinarily tricky, the implementation was surprisingly short. (But adding error checks would at least quadruple the code – probably more like 10x more code).

All in all it was a really fun and interesting exploration and I learned a lot.



I’ve always been fascinated and entertained by coin operated arcade games. Here are three views of the Rock-Ola World Series game, manufactured in 1934. It’s completely mechanical and quite amazing. It’s like a pinball machine. The game tracks balls and strikes, number of outs, runners advance, and a score is kept. In some ways, the device is a forerunner of the first mechanical computers.


Demo code. Replace “lt”, “gt”, “lte”, “gte” with Boolean operator symbols.

using System;

namespace IncBeta
{
  internal class IncBetaProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin IncBeta() from scratch demo \n");

      double a = 5.0; double b = 3.0; double x = 0.5;
      Console.WriteLine("\nI(x=0.5, a=5.0, b=3.0) via SciPy =" +
        " 0.22656250");
      double ib = IncBeta(x, a, b);
      Console.WriteLine("I(x=0.5, a=5.0, b=3.0) scratch  =  " +
        ib.ToString("F8"));

      a = 24.0; b = 36.0; x = 0.2;
      Console.WriteLine("\nI(x=0.2, a=24.0, b=36.0) via SciPy =" +
        " 0.00022272");
      ib = IncBeta(x, a, b);
      Console.WriteLine("I(x=0.2, a=24.0, b=36.0) scratch  =  " +
        ib.ToString("F8"));

      a = 60.0; b = 60.0; x = 0.7;
      Console.WriteLine("\nI(x=0.7, a=60.0, b=60.0) via SciPy =" +
        " 0.999997499205322");
      ib = IncBeta(x, a, b);
      Console.WriteLine("I(x=0.7, a=60.0, b=60.0) scratch  =  " +
        ib);

      Console.WriteLine("\nEnd demo ");
      Console.ReadLine();
    } // Main()

    static double ContFraction(double x, double a, double b,
      int maxTerms = 200)
    {
      // 1. pre-compute d values
      double[] d = new double[maxTerms];  // d[0] not used
      
      int end = maxTerms / 2;  // 10
      for (int m = 0; m "lt" end; ++m)  // [0,9]
      {
        int i = 2 * m;  // [0,18]  // even
        int j = i + 1;  // [1,19]  // odd
        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[4] = 1 + d[4] / 1;
      //t[3] = 1 + d[3] / t[4];
      //t[2] = 1 + d[2] / t[3];
      //t[1] = 1 + d[1] / t[2];
      //double cf = 1 / t[1];

      t[maxTerms - 1] = 1 + d[maxTerms - 1];
      for (int j = maxTerms - 2; j "gte" 1; --j)
        t[j] = 1 + d[j] / t[j + 1];

      // ShowVector(t);
      // Console.ReadLine();

      double cf = 1 / t[1];
      return cf;
    }

    static double IncBeta(double x, double a, double b)
    {
      // double top = Math.Pow(x, a) * Math.Pow((1 - x), b);
      // double bot = a * Beta(a, b);
      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 LogGamma(double z)
    {
      // plain Gamma() can easily overflow
      // 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 "lt" 0.5) // Gamma(z) = Pi / (Sin(Pi*z))* Gamma(1-z))
        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 "lt" coef.Length; ++i)
        sum += coef[i] / (zz + i);

      return (LogSqrtTwoPi + Math.Log(sum) - b) +
        (Math.Log(b) * (zz + 0.5));
    }

    static double LogBeta(double a, double b)
    {
      // plain Beta(a,b) = (Gamma(a) * Gamma(b)) / (Gamma(a+b))
      // but plain Gamma() can easily overflow, so:
      // Log( Beta(a,b) ) 
      //   = Log( (Gamma(a) * Gamma(b)) / (Gamma(a+b) ) 
      //   = LogGamma(a) + LogGamma(b) - LogGamma(a+b)
      double logbeta = LogGamma(a) + LogGamma(b) -
        LogGamma(a + b);
      return logbeta;
    }

    //static double Beta(double a, double b)
    //{
    //  return Math.Exp(LogBeta(a, b));
    //}

    static void ShowVector(double[] v)
    {
      for (int i = 0; i "lt" v.Length; ++i)
        Console.Write(v[i].ToString("F4") + " ");
      Console.WriteLine("");
    }

  } // Program

} // ns
This entry was posted in Miscellaneous. Bookmark the permalink.

Leave a Reply

Please log in using one of these methods to post your comment:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s