public static double incompleteBeta()

in samoa-api/src/main/java/org/apache/samoa/moa/core/Statistics.java [812:901]
73 lines of code
22 McCabe index (conditional complexity)

  public static double incompleteBeta(double aa, double bb, double xx) {

    double a, b, t, x, xc, w, y;
    boolean flag;

    if (aa <= 0.0 || bb <= 0.0)
      throw new ArithmeticException("ibeta: Domain error!");

    if ((xx <= 0.0) || (xx >= 1.0)) {
      if (xx == 0.0)
        return 0.0;
      if (xx == 1.0)
        return 1.0;
      throw new ArithmeticException("ibeta: Domain error!");
    }

    flag = false;
    if ((bb * xx) <= 1.0 && xx <= 0.95) {
      t = powerSeries(aa, bb, xx);
      return t;
    }

    w = 1.0 - xx;

    /* Reverse a and b if x is greater than the mean. */
    if (xx > (aa / (aa + bb))) {
      flag = true;
      a = bb;
      b = aa;
      xc = xx;
      x = w;
    } else {
      a = aa;
      b = bb;
      xc = w;
      x = xx;
    }

    if (flag && (b * x) <= 1.0 && x <= 0.95) {
      t = powerSeries(a, b, x);
      if (t <= MACHEP)
        t = 1.0 - MACHEP;
      else
        t = 1.0 - t;
      return t;
    }

    /* Choose expansion for better convergence. */
    y = x * (a + b - 2.0) - (a - 1.0);
    if (y < 0.0)
      w = incompleteBetaFraction1(a, b, x);
    else
      w = incompleteBetaFraction2(a, b, x) / xc;

    /*
     * Multiply w by the factor a b _ _ _ x (1-x) | (a+b) / ( a | (a) | (b) ) .
     */

    y = a * Math.log(x);
    t = b * Math.log(xc);
    if ((a + b) < MAXGAM && Math.abs(y) < MAXLOG && Math.abs(t) < MAXLOG) {
      t = Math.pow(xc, b);
      t *= Math.pow(x, a);
      t /= a;
      t *= w;
      t *= gamma(a + b) / (gamma(a) * gamma(b));
      if (flag) {
        if (t <= MACHEP)
          t = 1.0 - MACHEP;
        else
          t = 1.0 - t;
      }
      return t;
    }
    /* Resort to logarithms. */
    y += t + lnGamma(a + b) - lnGamma(a) - lnGamma(b);
    y += Math.log(w / a);
    if (y < MINLOG)
      t = 0.0;
    else
      t = Math.exp(y);

    if (flag) {
      if (t <= MACHEP)
        t = 1.0 - MACHEP;
      else
        t = 1.0 - t;
    }
    return t;
  }