@@ -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
@@ -945,7 +934,8 @@ function __gamma_inc_inv(a::Float64, minpq::Float64, pcase::Bool)
945
934
haseta = false
946
935
947
936
logp = pcase ? log (minpq) : log1p (- minpq)
948
- logr = (1.0 / a)* (logp + logabsgamma (a + 1.0 )[1 ])
937
+ loggamma1pa = a <= 1.0 ? loggamma1p (a) : loggamma (a + 1.0 )
938
+ logr = (logp + loggamma1pa) / a
949
939
if logr < log (0.2 * (1 + a)) # small value of p
950
940
x0 = gamma_inc_inv_psmall (a, logr)
951
941
elseif ! pcase && minpq < min (0.02 , exp (- 1.5 * a)/ gamma (a)) && a < 10 # small q
@@ -955,7 +945,7 @@ function __gamma_inc_inv(a::Float64, minpq::Float64, pcase::Bool)
955
945
elseif abs (a - 1.0 ) < 1.0e-4
956
946
x0 = pcase ? - log1p (- minpq) : - log (minpq)
957
947
elseif a < 1.0 # small value of a
958
- x0 = gamma_inc_inv_asmall (a, logp )
948
+ x0 = exp (logr )
959
949
else # large a
960
950
haseta = true
961
951
x0, fp = gamma_inc_inv_alarge (a, minpq, pcase)
0 commit comments