gammapInv function

double gammapInv(num p, num a)

Implementation

double gammapInv(num p, num a) {
  const epsilon = 1.0e-8;
  final a1 = a - 1.0;
  var gln = gammaLn(a);
  var x = 0.0;
  var afac = 0.0;
  var lna1 = 0.0;

  if (p >= 1.0) {
    return max(100, a + 100 * sqrt(a));
  } else if (p <= 0) {
    return 0.0;
  } else if (a > 1.0) {
    lna1 = log(a1);
    afac = exp(a1 * (lna1 - 1.0) - gln);
    final pp = p < 0.5 ? p : 1.0 - p;
    final t = sqrt(-2 * log(pp));
    x = (2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t;
    if (p < 0.5) {
      x = -x;
    }
    x = max(1.0e-3,
        a * pow(1.0 - 1.0 / (9.0 * a) - x / (3.0 * sqrt(a)), 3).toDouble());
  } else {
    final t = 1.0 - a * (0.253 + a * 0.12);
    if (p < t) {
      x = pow(p / t, 1.0 / a).toDouble();
    } else {
      x = 1.0 - log(1.0 - (p - t) / (1.0 - t));
    }
  }
  for (var j = 0; j < 12; j++) {
    if (x <= 0.0) {
      return 0.0;
    }
    final err = lowRegGamma(a, x) - p;
    var t = a > 1.0
        ? afac * exp(-(x - a1) + a1 * (log(x) - lna1))
        : exp(-x + a1 * log(x) - gln);
    final u = err / t;
    x -= (t = u / (1.0 - 0.5 * min(1.0, u * ((a - 1.0) / x - 1.0))));
    if (x <= 0.0) {
      x = 0.5 * (x + t);
    }
    if (t.abs() < epsilon * x) {
      break;
    }
  }
  return x;
}