(double a, double x)
| 806 | } |
| 807 | |
| 808 | private static double lnLowIncGamma(double a, double x) |
| 809 | { |
| 810 | |
| 811 | /* |
| 812 | * We compute the log of the lower incomplete gamma function by taking the log of |
| 813 | * |
| 814 | * oo |
| 815 | * ===== |
| 816 | * -x a \ Gamma(a) n |
| 817 | * y(a, x) = e x > ---------------- x |
| 818 | * / Gamma(a + 1 + n) |
| 819 | * ===== |
| 820 | * n = 0 |
| 821 | * |
| 822 | * Which becomes |
| 823 | * |
| 824 | * |
| 825 | * / oo \ |
| 826 | * |===== | |
| 827 | * |\ Gamma(a) n| |
| 828 | * log(y(a, x)) = -x + log(x) a + log| > ---------------- x | |
| 829 | * |/ Gamma(a + 1 + n) | |
| 830 | * |===== | |
| 831 | * \n = 0 / |
| 832 | * |
| 833 | * To reduce over flow of the gammas and exponentation in the summation, we compute the sum as |
| 834 | * |
| 835 | * oo |
| 836 | * ===== |
| 837 | * \ LogGamma(a) - LogGamma(a + 1 + n) + ln(x) n |
| 838 | * > e |
| 839 | * / |
| 840 | * ===== |
| 841 | * n = 0 |
| 842 | * |
| 843 | * Testin with x=0.5 to 100 (in increments of 0.5) |
| 844 | * a=0.15, max relative error is ~ 6e-13, with x < 30 having a relative error smaller than 5e-14 |
| 845 | * a=0.5, maximum relative and absolute eror is ~1.5e-12 and ~1.5e-12 repsectivly, the error starts getting above e-15 when x = 25, and x=50 the error is up to 2.2e-13 |
| 846 | * a=10, maximum relative and absolute eror is ~4e-14 and ~4e-13 repsectivly, the error starts getting abovee-15 when x = 49. For x near zero (up to x ~ 2.5) the error is higher before droping, ~10^-14 |
| 847 | * a=25, maximum relative error is ~9.99e-15. From x~ 0 to 14, the error is worse, then droping to ~e-16 . |
| 848 | * a=50, accuracy starting to degrad badly. From x~ 0 to 18 the error goes from 1.3 to 1e-7, the erro grows at an exponential rate as x -> 0. As x-> Infinity the error gets back down to ~5e-16 |
| 849 | */ |
| 850 | |
| 851 | //Sumation first |
| 852 | |
| 853 | double lnGa = lnGamma(a); |
| 854 | /** |
| 855 | * This value will be updated by the property Gamma(z+1) = Gamma(z) * z, which - when taken the log of, is <br> |
| 856 | * LnGamma(z+1) = LnGamma(z) + ln(z) |
| 857 | */ |
| 858 | double lnGan = lnGa+log(a); |
| 859 | double n = 0; |
| 860 | /** |
| 861 | * this is the n * ln(x) term. it will be updated by adding the log of x at each step |
| 862 | */ |
| 863 | double lnXN = 0; |
| 864 | double lnX = log(x); |
| 865 |
no test coverage detected