protected static double gammaQuantileNewtonRefinement(final double logpt,final double k,final double theta,final int maxit,double x){
  final double EPS_N=1e-15;
  if (x <= 0) {
    x=Double.MIN_NORMAL;
  }
  double logpc=logcdf(x,k,theta);
  if (x == Double.MIN_NORMAL && logpc > logpt * (1. + 1e-7)) {
    return 0.;
  }
  if (logpc == Double.NEGATIVE_INFINITY) {
    return 0.;
  }
  for (int i=0; i < maxit; i++) {
    final double logpe=logpc - logpt;
    if (Math.abs(logpe) < Math.abs(EPS_N * logpt)) {
      break;
    }
    final double g=logpdf(x,k,theta);
    if (g == Double.NEGATIVE_INFINITY) {
      break;
    }
    final double newx=x - logpe * Math.exp(logpc - g);
    logpc=logcdf(newx,k,theta);
    if (Math.abs(logpc - logpt) > Math.abs(logpe) || (i > 0 && Math.abs(logpc - logpt) == Math.abs(logpe))) {
      break;
    }
    x=newx;
  }
  return x;
}

protected static double gammaQuantileNewtonRefinement ( final double logpt , final double k , final double theta , final int maxit , double x ) { final double EPS_N = 1e-15 ; if ( x <= 0 ) { x = Double . MIN_NORMAL ; } double logpc = logcdf ( x , k , theta ) ; if ( x == Double . MIN_NORMAL && logpc > logpt * ( 1. + 1e-7 ) ) { return 0. ; } if ( logpc == Double . NEGATIVE_INFINITY ) { return 0. ; } for ( int i = 0 ; i < maxit ; i ++ ) { final double logpe = logpc - logpt ; if ( Math . abs ( logpe ) < Math . abs ( EPS_N * logpt ) ) { break ; } final double g = logpdf ( x , k , theta ) ; if ( g == Double . NEGATIVE_INFINITY ) { break ; } final double newx = x - logpe * Math . exp ( logpc - g ) ; logpc = logcdf ( newx , k , theta ) ; if ( Math . abs ( logpc - logpt ) > Math . abs ( logpe ) || ( i > 0 && Math . abs ( logpc - logpt ) == Math . abs ( logpe ) ) ) { break ; } x = newx ; } return x ; }
refinement of chisquared probit using newton iterations . a trick used by gnu r to improve precision .
