@@ -724,17 +724,6 @@ function gamma_inc_inv_qsmall(a::Float64, q::Float64)
724
724
return x0
725
725
end
726
726
727
- """
728
- gamma_inc_inv_asmall(a, logp)
729
-
730
- Compute `x0` - initial approximation when `a` is small.
731
- Here the solution `x` of ``P(a,x)=p=\\ exp(logp)`` satisfies ``x_{l} < x < x_{u}``
732
- where ``x_{l} = (p\\ Gamma(a+1))^{1/a}`` and ``x_{u} = -\\ log{(1 - p\\ Gamma(a+1))}``, and is used as starting value for Newton iteration.
733
- """
734
- function gamma_inc_inv_asmall (a:: Float64 , logp:: Float64 )
735
- return exp ((logp + loggamma1p (a)) / a)
736
- end
737
-
738
727
"""
739
728
gamma_inc_inv_alarge(a, minpq, pcase)
740
729
@@ -941,7 +930,8 @@ function __gamma_inc_inv(a::Float64, minpq::Float64, pcase::Bool)
941
930
haseta = false
942
931
943
932
logp = pcase ? log (minpq) : log1p (- minpq)
944
- logr = (1.0 / a)* (logp + logabsgamma (a + 1.0 )[1 ])
933
+ loggamma1pa = a <= 1.0 ? loggamma1p (a) : loggamma (a + 1.0 )
934
+ logr = (logp + loggamma1pa) / a
945
935
if logr < log (0.2 * (1 + a)) # small value of p
946
936
x0 = gamma_inc_inv_psmall (a, logr)
947
937
elseif ! pcase && minpq < min (0.02 , exp (- 1.5 * a)/ gamma (a)) && a < 10 # small q
@@ -951,7 +941,7 @@ function __gamma_inc_inv(a::Float64, minpq::Float64, pcase::Bool)
951
941
elseif abs (a - 1.0 ) < 1.0e-4
952
942
x0 = pcase ? - log1p (- minpq) : - log (minpq)
953
943
elseif a < 1.0 # small value of a
954
- x0 = gamma_inc_inv_asmall (a, logp )
944
+ x0 = exp (logr )
955
945
else # large a
956
946
haseta = true
957
947
x0, fp = gamma_inc_inv_alarge (a, minpq, pcase)
0 commit comments