protected static double regularizedIncBetaQuadrature(double alpha,double beta,double x){
  final double alphapbeta=alpha + beta;
  final double a1=alpha - 1.0;
  final double b1=beta - 1.0;
  final double mu=alpha / alphapbeta;
  final double lnmu=Math.log(mu);
  final double lnmuc=Math.log1p(-mu);
  double t=Math.sqrt(alpha * beta / (alphapbeta * alphapbeta * (alphapbeta + 1.0)));
  final double xu;
  if (x > alpha / alphapbeta) {
    if (x >= 1.0) {
      return 1.0;
    }
    xu=Math.min(1.0,Math.max(mu + 10.0 * t,x + 5.0 * t));
  }
 else {
    if (x <= 0.0) {
      return 0.0;
    }
    xu=Math.max(0.0,Math.min(mu - 10.0 * t,x - 5.0 * t));
  }
  double sum=0.0;
  for (int i=0; i < GAUSSLEGENDRE_Y.length; i++) {
    t=x + (xu - x) * GAUSSLEGENDRE_Y[i];
    sum+=GAUSSLEGENDRE_W[i] * Math.exp(a1 * (Math.log(t) - lnmu) + b1 * (Math.log1p(-t) - lnmuc));
  }
  double ans=sum * (xu - x) * Math.exp(a1 * lnmu - GammaDistribution.logGamma(alpha) + b1 * lnmuc - GammaDistribution.logGamma(beta) + GammaDistribution.logGamma(alphapbeta));
  return ans > 0 ? 1.0 - ans : -ans;
}

protected static double regularizedIncBetaQuadrature ( double alpha , double beta , double x ) { final double alphapbeta = alpha + beta ; final double a1 = alpha - 1.0 ; final double b1 = beta - 1.0 ; final double mu = alpha / alphapbeta ; final double lnmu = Math . log ( mu ) ; final double lnmuc = Math . log1p ( - mu ) ; double t = Math . sqrt ( alpha * beta / ( alphapbeta * alphapbeta * ( alphapbeta + 1.0 ) ) ) ; final double xu ; if ( x > alpha / alphapbeta ) { if ( x >= 1.0 ) { return 1.0 ; } xu = Math . min ( 1.0 , Math . max ( mu + 10.0 * t , x + 5.0 * t ) ) ; } else { if ( x <= 0.0 ) { return 0.0 ; } xu = Math . max ( 0.0 , Math . min ( mu - 10.0 * t , x - 5.0 * t ) ) ; } double sum = 0.0 ; for ( int i = 0 ; i < GAUSSLEGENDRE_Y . length ; i ++ ) { t = x + ( xu - x ) * GAUSSLEGENDRE_Y [ i ] ; sum += GAUSSLEGENDRE_W [ i ] * Math . exp ( a1 * ( Math . log ( t ) - lnmu ) + b1 * ( Math . log1p ( - t ) - lnmuc ) ) ; } double ans = sum * ( xu - x ) * Math . exp ( a1 * lnmu - GammaDistribution . logGamma ( alpha ) + b1 * lnmuc - GammaDistribution . logGamma ( beta ) + GammaDistribution . logGamma ( alphapbeta ) ) ; return ans > 0 ? 1.0 - ans : - ans ; }
returns the regularized incomplete beta function i_x ( a , b ) by quadrature , based on the book " numerical recipes " .
