protected static double gammaQuantileNewtonRefinement ( final double logpt , final double k , final double theta , final int maxit , double x ) { final double EPS_N = 1e-15 ; // Precision threshold // 0 is not possible, try MIN_NORMAL instead if ( x <= 0 ) { x = Double . MIN_NORMAL ; } // Current estimation 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. ; } // Refine by newton iterations for ( int i = 0 ; i < maxit ; i ++ ) { // Error of current approximation final double logpe = logpc - logpt ; if ( Math . abs ( logpe ) < Math . abs ( EPS_N * logpt ) ) { break ; } // Step size is controlled by PDF: final double g = logpdf ( x , k , theta ) ; if ( g == Double . NEGATIVE_INFINITY ) { break ; } final double newx = x - logpe * FastMath . exp ( logpc - g ) ; // New estimate: logpc = logcdf ( newx , k , theta ) ; if ( Math . abs ( logpc - logpt ) > Math . abs ( logpe ) || ( i > 0 && Math . abs ( logpc - logpt ) == Math . abs ( logpe ) ) ) { // no further improvement break ; } x = newx ; } return x ; }
Refinement of ChiSquared probit using Newton iterations .