ibetaInv function

double ibetaInv(double a, double b, double p)

Computes inverse of incomplete beta function https://malishoaib.wordpress.com/2014/05/30/inverse-of-incomplete-beta-function-computational-statisticians-wet-dream/

Implementation

double ibetaInv(double a, double b, double p) {
  final double a1 = a - 1.0;
  final double b1 = b - 1.0;
  const double error = 1e-8;

  if (p <= 0.0) return 0.0;
  if (p >= 1.0) return 1.0;

  double x;

  if (a >= 1.0 && b >= 1.0) {
    double pp;

    if (p < 0.5)
      pp = p;
    else
      pp = 1.0 - p;

    double t = math.sqrt(-2 * math.log(pp));
    x = (2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t;
    if (p < 0.5) x = -x;
    double al = ((x * x) - 3.0) / 6.0;
    double h = 2.0 / (1.0 / (2.0 * a - 1.0) + 1.0 / (2.0 * b - 1.0));
    double w = (x * math.sqrt(al + h) / h) -
        (1.0 / (2.0 * b - 1) - 1.0 / (2.0 * a - 1.0)) *
            (al + 5.0 / 6.0 - 2.0 / (3.0 * h));
    x = a / (a + b * math.exp(2.0 * w));
  } else {
    double lna = math.log(a / (a + b));
    double lnb = math.log(b / (a + b));
    double t = math.exp(a * lna) / a;
    double u = math.exp(b * lnb) / b;
    double w = t + u;
    if (p < t / w)
      x = math.pow(a * w * p, 1.0 / a).toDouble();
    else
      x = 1.0 - math.pow(b * w * (1.0 - p), 1.0 / b);
  }

  double afac = -lgamma(a).lgamma - lgamma(b).lgamma + lgamma(a + b).lgamma;
  double j = 0.0;
  for (int i = 0; i < 10; i++) {
    if (x == 1.0) return x;
    double err = ibetaReg(a, b, x) - p;
    double t = math.exp(a1 * math.log(x) + b1 * math.log(1.0 - x) + afac);
    double u = err / t;
    t = u / (1.0 - 0.5 * math.min(1.0, u * (a1 / x - b1 / (1.0 - x))));
    x -= t;
    if (x <= 0.0) x = 0.5 * (x + t);
    if (x >= 1.0) x = 0.5 * (x + t + 1.0);
    if ((t.abs() < error * x) && (j > 0)) break;
  }

  return x;
}