betaRegularized function

double betaRegularized(double a, double b, double x)

Returns the regularized lower incomplete beta function I_x(a,b) = 1/Beta(a,b) * int(t^(a-1)*(1-t)^(b-1),t=0..x) for real a > 0, b > 0, 1 >= x >= 0.

Implementation

double betaRegularized(double a, double b, double x) {
  if (a < 0.0) {
    throw ArgumentError.value(a, 'a', messages.argumentNotNegative);
  }

  if (b < 0.0) {
    throw ArgumentError.value(b, 'b', messages.argumentNotNegative);
  }

  if (x < 0.0 || x > 1.0) {
    throw ArgumentError.value(x, 'x', messages.argumentInIntervalXYInclusive);
  }

  var bt = (x == 0.0 || x == 1.0)
      ? 0.0
      : exp(gammaLn(a + b) -
          gammaLn(a) -
          gammaLn(b) +
          (a * log(x)) +
          (b * log(1.0 - x)));

  var symmetryTransformation = x >= (a + 1.0) / (a + b + 2.0);

  /* Continued fraction representation */
  var eps = doublePrecision;
  var fpmin = increment(0.0) / eps;

  if (symmetryTransformation) {
    x = 1.0 - x;
    var swap = a;
    a = b;
    b = swap;
  }

  var qab = a + b;
  var qap = a + 1.0;
  var qam = a - 1.0;
  var c = 1.0;
  var d = 1.0 - (qab * x / qap);

  if (d.abs() < fpmin) {
    d = fpmin;
  }

  d = 1.0 / d;
  var h = d;

  for (int m = 1, m2 = 2; m <= 50000; m++, m2 += 2) {
    var aa = m * (b - m) * x / ((qam + m2) * (a + m2));
    d = 1.0 + (aa * d);

    if (d.abs() < fpmin) {
      d = fpmin;
    }

    c = 1.0 + (aa / c);
    if (c.abs() < fpmin) {
      c = fpmin;
    }

    d = 1.0 / d;
    h *= d * c;
    aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
    d = 1.0 + (aa * d);

    if (d.abs() < fpmin) {
      d = fpmin;
    }

    c = 1.0 + (aa / c);

    if (c.abs() < fpmin) {
      c = fpmin;
    }

    d = 1.0 / d;
    var del = d * c;
    h *= del;

    if ((del - 1.0).abs() <= eps) {
      return symmetryTransformation ? 1.0 - (bt * h / a) : bt * h / a;
    }
  }

  return symmetryTransformation ? 1.0 - (bt * h / a) : bt * h / a;
}