Finds the value x such that #gammaP(double,double) P(a,x) = p . @param a any real value @param p and value in the range [0, 1] @return the inverse
(double p, double a)
| 703 | * @return the inverse |
| 704 | */ |
| 705 | public static double invGammaP(double p, double a) |
| 706 | { |
| 707 | //Method from Numerical Recipies 3rd edition p 263 |
| 708 | if(p < 0 || p > 1) |
| 709 | throw new ArithmeticException("Probability p must be in the range [0,1], "+ p + "is not valid"); |
| 710 | //First an estimate must be obtained |
| 711 | double am1 = a-1; |
| 712 | double lnGamA = lnGamma(a); |
| 713 | double x;//the to be returned |
| 714 | double afac = 1;//ONLY used when a>1 |
| 715 | double lna1 =1;//ONLY used when a>1 |
| 716 | if(a > 1) |
| 717 | { |
| 718 | lna1 = log(am1); |
| 719 | afac = exp(am1*(lna1-1)-lnGamA); |
| 720 | double pp = (p < 0.5) ? p : 1-p; |
| 721 | double t = sqrt(-2*log(pp)); |
| 722 | //Now our inital estimate |
| 723 | x = (2.30753+t*0.27061)/(1.+t*(0.99229+t*0.04481)) - t; |
| 724 | if(p < 0.5) |
| 725 | x = -x; |
| 726 | x = max(1e-3, a * pow(1.0 - (1.0/(9.*a)) - x/(3.*sqrt(a)) , 3) );//if the estimate is too small, increase it |
| 727 | |
| 728 | } |
| 729 | else |
| 730 | { |
| 731 | double t = 1.0 - a*(0.253+a*0.12); |
| 732 | if(p < t) |
| 733 | x = pow(p/t, 1.0/a); |
| 734 | else |
| 735 | x = 1-log(1-(p-t)/(1.0-t)); |
| 736 | } |
| 737 | |
| 738 | //Estimate obtained, now refinement |
| 739 | for(int j = 0; j < 12; j++) |
| 740 | { |
| 741 | if (x <= 0)//x is very small, return 0 b/c rounding errors and loss of precision will make accuracy imposible |
| 742 | return 0; |
| 743 | double err = gammaP(a, x) - p; |
| 744 | double t; |
| 745 | if(a > 1) |
| 746 | t = afac*exp(-(x-am1)+am1*(log(x)-lna1)); |
| 747 | else |
| 748 | t = exp(-x+am1*log(x)-lnGamA); |
| 749 | |
| 750 | double u = err/t; |
| 751 | t = u / (1-0.5*min(1 , u*(am1/x - 1)) );//Halley's method |
| 752 | x -= t; |
| 753 | if(x <= 0)//bais x from going negative (means we are getting a bad value) |
| 754 | x = (x+t)/2; |
| 755 | if(abs(t) < 1e-8*x) |
| 756 | break;//the error is the (1e-8)^2, if we reach this point we have already converged |
| 757 | } |
| 758 | |
| 759 | return x; |
| 760 | } |
| 761 | |
| 762 | /** |