Merge branch 'issue-9-normalized-incomplete-gamma' into 'master'

Fix #9: Implement normalized incomplete gamma functions

Closes #9

See merge request !13
parents c1c007eb 3f0221c3
Pipeline #2358 passed with stage
in 2 minutes and 34 seconds
 ... ... @@ -10,6 +10,7 @@ before_script: test: artifacts: when: always paths: - test.log script: ... ...
 ... ... @@ -130,8 +130,34 @@ supports these: - log-gamma - gamma - psi - incomplete-gamma-lower - incomplete-gamma-lower a z &key normalized-p If normalized-p is NIL, return the [value](https://dlmf.nist.gov/8.2.E1): math \gamma(a, z) = \int_0^z t^{a-1} e^{-t} dt  Otherwise return the normalized [value](https://dlmf.nist.gov/8.2.E4): math P(a, z) = \frac{1}{\Gamma(a)} \int_0^z t^{a-1} e^{-t} dt  - incomplete-gamma-upper If normalized-p is NIL, return the [value](https://dlmf.nist.gov/8.2.E2): math \Gamma(a, z) = \int_z^{\infty} t^{a-1} e^{-t} dt  Otherwise return the normalized [value]https://dlmf.nist.gov/8.2.E4): math Q(a, z) = \frac{1}{\Gamma(a)} \int_z^{\infty} t^{a-1} e^{-t} dt  - incomplete-gamma-lower-normalized a z Returns the [value](https://dlmf.nist.gov/8.2.E4): math P(a, z) = \frac{1}{\Gamma(a)} \int_0^z t^{a-1} e^{-t} dt  - incomplete-gamma-upper-normalized a z Returns the [value](https://dlmf.nist.gov/8.2.E4): math Q(a, z) = \frac{1}{\Gamma(a)} \int_z^{\infty} t^{a-1} e^{-t} dt  - fresnel-s - fresnel-c - sin-integral ... ...
 ... ... @@ -474,54 +474,172 @@ return (* sum (expt z a))))) ;; Tail of the incomplete gamma function. (defun incomplete-gamma-upper (a z) (defun incomplete-gamma-upper (a z &key (normalized-p nil)) "Upper tail of the incomplete gamma function defined by: integrate(t^(a-1)*exp(-t), t, z, inf)" (with-floating-point-contagion (a z) (if (and (realp a) (<= a 0)) ;; incomplete_gamma_tail(v, z) = z^v*exp_integral_e(1-a,z) (* (expt z a) (exp-integral-e (- 1 a) z)) (if (and (zerop (imagpart a)) (zerop (imagpart z))) ;; For real values, we split the result to compute either the ;; tail directly from the continued fraction or from gamma(a) ;; - incomplete-gamma. The continued fraction doesn't ;; converge on the negative real axis, so we can't use that ;; there. And accuracy appears to be better if z is "small". ;; We take this to mean |z| < |a-1|. Note that |a-1| is the ;; peak of the integrand. (if (and (> (abs z) (abs (- a 1))) (not (minusp (realpart z)))) (cf-incomplete-gamma-tail a z) (- (gamma a) (cf-incomplete-gamma a z))) ;; If the argument is close enough to the negative real axis, ;; the continued fraction for the tail is not very accurate. ;; Use the incomplete gamma function to evaluate in this ;; region. (Arbitrarily selected the region to be a sector. ;; But what is the correct size of this sector?) (if (<= (abs (phase z)) 3.1) (cf-incomplete-gamma-tail a z) (- (gamma a) (cf-incomplete-gamma a z))))))) (defun incomplete-gamma-lower (a z) integrate(t^(a-1)*exp(-t), t, z, inf) If NORMALIZED-P is non-NIL, return the normalized value: integrate(t^(a-1)*exp(-t), t, z, inf) / gamma(a)" (if normalized-p (incomplete-gamma-upper-normalized a z) (with-floating-point-contagion (a z) (if (and (realp a) (<= a 0)) ;; incomplete_gamma_tail(v, z) = z^v*exp_integral_e(1-a,z) (* (expt z a) (exp-integral-e (- 1 a) z)) (if (and (zerop (imagpart a)) (zerop (imagpart z))) ;; For real values, we split the result to compute either the ;; tail directly from the continued fraction or from gamma(a) ;; - incomplete-gamma. The continued fraction doesn't ;; converge on the negative real axis, so we can't use that ;; there. And accuracy appears to be better if z is "small". ;; We take this to mean |z| < |a-1|. Note that |a-1| is the ;; peak of the integrand. (if (and (> (abs z) (abs (- a 1))) (not (minusp (realpart z)))) (cf-incomplete-gamma-tail a z) (- (gamma a) (cf-incomplete-gamma a z))) ;; If the argument is close enough to the negative real axis, ;; the continued fraction for the tail is not very accurate. ;; Use the incomplete gamma function to evaluate in this ;; region. (Arbitrarily selected the region to be a sector. ;; But what is the correct size of this sector?) (if (<= (abs (phase z)) 3.1) (cf-incomplete-gamma-tail a z) (- (gamma a) (cf-incomplete-gamma a z)))))))) (defun incomplete-gamma-lower (a z &key (normalized-p nil)) "Incomplete gamma function defined by: integrate(t^(a-1)*exp(-t), t, 0, z)" integrate(t^(a-1)*exp(-t), t, 0, z) If NORMALIZED-P is non-NIL, return the normalized value: integrate(t^(a-1)*exp(-t), t, 0, z) / gamma(a)" (if normalized-p (incomplete-gamma-lower-normalized a z) (with-floating-point-contagion (a z) (if (and (< (abs a) 1) (< (abs z) 1)) (s-incomplete-gamma a z) (if (and (realp a) (realp z)) (if (< z (- a 1)) (cf-incomplete-gamma a z) (- (gamma a) (cf-incomplete-gamma-tail a z))) ;; The continued fraction doesn't converge very fast if a ;; and z are small. In this case, use the series ;; expansion instead, which converges quite rapidly. (if (< (abs z) (abs a)) (cf-incomplete-gamma a z) (- (gamma a) (cf-incomplete-gamma-tail a z)))))))) ;; The continued fraction for the normalized incomplete gamma upper ;; function, Q(a, z), valid for all z except for the negative real ;; axis: ;; ;; Q(a, z) = z^a*exp(-z)/gamma(a) * 1/cf(a,z) ;; ;; where cf is the continued fraction: ;; ;; b[n] = 2*n+1-a+z ;; a[n] = n*(a-n) ;; ;; See http://functions.wolfram.com/06.08.10.0003.01 (defun cf-incomplete-gamma-upper-normalized (a z) "Continued fraction for the normalized incomplete gamma upper function" (with-floating-point-contagion (a z) (let* ((az (- z a)) (cf (lentz #'(lambda (n) (+ n n 1 az)) #'(lambda (n) (* n (- a n))))) (scale-term (handler-case (/ (* (expt z a) (exp (- z))) (gamma a)) (arithmetic-error () ;; z^a*exp(-z)/gamma(a) can overflow. Use logs to compute ;; the result: exp(a*log(z) - z - log gamma(a)) (exp (- (* a (log z)) z (log-gamma a))))))) ;; If the continued fraction returns a negative value when both ;; a and z are positive real numbers, the continued fraction ;; didn't actually converge. (a = 200, z = 100 is one such ;; case). In that case, we can use ;; incomplete-gamma-lower-normalized to compute the final ;; value. (cond ((and (realp a) (realp z) (plusp a) (plusp z) (realp cf) (minusp cf)) ;; Use the fact that P(a,z) + Q(a,z) = 1. (- 1 (incomplete-gamma-lower-normalized a z))) (t (/ scale-term cf)))))) (defun incomplete-gamma-upper-normalized (a z) "Normalized incomplete gamma upper function defined by: Q(a,z) = 1/gamma(a)*integrate(t^(a-1)*exp(-t), t, z, inf). This is the incomplete-gamma-upper function divided by gamma(a)" ;; Using maxima, ;; plot2d(gamma_incomplete_regularized(a,a-0.375),[z,100,1000]) ;; shows that 0.5018 >= Q(a,z) >= 0.5002. As z increases, Q(a,z) ;; decreases. Thus, when z <= a-0.375, use 1 - P(a,z) (if (and (realp a) (realp z) (<= z (- a 0.375))) (- 1 (incomplete-gamma-lower-normalized a z)) (cf-incomplete-gamma-upper-normalized a z))) ;; The continued fraction for the normalized incomplete gamma lower ;; function, P(a, z). First note that ;; ;; Q(a, z) = 1 - z^a*exp(-z)/gamma(a) * 1/cf(a,z) ;; ;; where cf is the continued fraction: ;; ;; b = a ;; b[n] = a + z + n ;; a[n] = -(a + n - 1)*z ;; ;; See http://functions.wolfram.com/06.08.10.0007.01 ;; ;; Since P(a, z) = 1 - Q(a, z), we have a continued fraction for P(a, z) (defun cf-incomplete-gamma-lower-normalized (a z) (with-floating-point-contagion (a z) (if (and (< (abs a) 1) (< (abs z) 1)) (s-incomplete-gamma a z) (if (and (realp a) (realp z)) (if (< z (- a 1)) (cf-incomplete-gamma a z) (- (gamma a) (cf-incomplete-gamma-tail a z))) ;; The continued fraction doesn't converge very fast if a ;; and z are small. In this case, use the series ;; expansion instead, which converges quite rapidly. (if (< (abs z) (abs a)) (cf-incomplete-gamma a z) (- (gamma a) (cf-incomplete-gamma-tail a z))))))) (let* ((-z (- z)) (a+z (+ a z)) (cf (lentz #'(lambda (n) (if (zerop n) a (+ a+z n))) #'(lambda (n) (* -z (+ a n -1))))) (scale-term (handler-case (/ (* (expt z a) (exp -z)) (gamma a)) (arithmetic-error () ;; z^a*exp(-z)/gamma(a) can overflow. Use logs to compute ;; the result: exp(a*log(z) - z - log gamma(a)) (exp (- (* a (log z)) z (log-gamma a))))))) (/ scale-term cf)))) (defun incomplete-gamma-lower-normalized (a z) "Normalized incomplete gamma lower function defined by: P(a,z) = 1/gamma(a)*integrate(t^(a-1)*exp(-t), t, 0, z) This is the incomplete-gamma-lower function divided by gamma(a)" ;; Assume that if a > z, incomplete-gamma-lower is greater than 1/2, ;; so that we can achieve greater accuracy using P(a,z) = 1 - ;; Q(a,z). (if (and (realp a) (realp z) (> z (- a 0.375))) (- 1 (incomplete-gamma-upper-normalized a z)) (cf-incomplete-gamma-lower-normalized a z))) (defun erf (z) "Error function: ... ...
 ... ... @@ -1778,3 +1778,66 @@ (check-accuracy 198 (erfc #q-4) #q1.9999999845827420997199811478403265131159514278547464108088316570950057869589732) nil) (rt:deftest incomplete-gamma-upper-normalized.d.a2 ;; incomplete-gamma-upper-normalized(2, z) = (1 + z)*exp(-z) (loop for k from 1 below 10 for z = (float k 1d0) for y = (incomplete-gamma-upper-normalized 2 z) for true = (* (+ z 1) (exp (- z))) for result = (check-accuracy 52.33 y true) when result append (list (list (list k z) result))) nil) (rt:deftest incomplete-gamma-upper-normalized.q.a2 ;; incomplete-gamma-upper-normalized(2, z) = (1 + z)*exp(-z) (loop for k from 1 below 10 for z = (float k #q1q0) for y = (incomplete-gamma-upper-normalized 2 z) for true = (* (+ z 1) (exp (- z))) for result = (check-accuracy 203.70 y true) when result append (list (list (list k z) result))) nil) (rt:deftest incomplete-gamma-lower-normalized.d.a2 ;; incomplete-gamma-lower-normalized(2, z) = 1-(1 + z)*exp(-z) (loop for k from 1 below 10 for z = (float k 1d0) for y = (incomplete-gamma-lower-normalized 2 z) for true = (- 1 (* (+ z 1) (exp (- z)))) for result = (check-accuracy 49.82 y true) when result append (list (list (list k z) result))) nil) (rt:deftest incomplete-gamma-lower-normalized.q.a2 ;; incomplete-gamma-lower-normalized(2, z) = 1-(1 + z)*exp(-z) (loop for k from 1 below 10 for z = (float k #q1q0) for y = (incomplete-gamma-lower-normalized 2 z) for true = (- 1 (* (+ z 1) (exp (- z)))) for result = (check-accuracy 203.83 y true) when result append (list (list (list k z) result))) nil) ;; incomplete-gamma-upper-normalized(200, 100) = ;; ;; exp(-100*36451035284924577938246208798747009164319474757880246359883694555113407009453436064573518999387789077985197279221655719227002367495061633272603038249747260895707250896595889294145309676586627989388740458641362406969609459453916777341749316070359589697827702813520519796940239276744754778199440304584107317957027129587503199/1356006206645357299077422810994072904566969809700681604285727988319939931024001696953196916719184549697395496290863162742676361760549235149195411231740418104602504325580502523311497039304043141691060121240640609954226541318710631103275528465092597490136227936213123455950399178299 ;; ;; according to maxima. (rt:deftest incomplete-gamma-upper-normalized.d.a200 (let ((true (* (exp -100d0) 36451035284924577938246208798747009164319474757880246359883694555113407009453436064573518999387789077985197279221655719227002367495061633272603038249747260895707250896595889294145309676586627989388740458641362406969609459453916777341749316070359589697827702813520519796940239276744754778199440304584107317957027129587503199/1356006206645357299077422810994072904566969809700681604285727988319939931024001696953196916719184549697395496290863162742676361760549235149195411231740418104602504325580502523311497039304043141691060121240640609954226541318710631103275528465092597490136227936213123455950399178299)) (result (incomplete-gamma-upper-normalized 200d0 100d0))) (check-accuracy 59.89 result true)) nil) (rt:deftest incomplete-gamma-upper-normalized.q.a200 (let ((true (* (exp #q-100q0) 36451035284924577938246208798747009164319474757880246359883694555113407009453436064573518999387789077985197279221655719227002367495061633272603038249747260895707250896595889294145309676586627989388740458641362406969609459453916777341749316070359589697827702813520519796940239276744754778199440304584107317957027129587503199/1356006206645357299077422810994072904566969809700681604285727988319939931024001696953196916719184549697395496290863162742676361760549235149195411231740418104602504325580502523311497039304043141691060121240640609954226541318710631103275528465092597490136227936213123455950399178299)) (result (incomplete-gamma-upper-normalized #q200 #q100))) (check-accuracy 211.99 result true)) nil)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!