Commit 3f0221c3 authored by Raymond Toy's avatar Raymond Toy
Browse files

Fix #9: Implement normalized incomplete gamma functions

Adds two new functions:

- incomplete-gamma-upper-normalized
- incomplete-gamma-lower-normalized

These are the corresponding incomplete gamma functions but normalized
by gamma(a).
parent c1c007eb
......@@ -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[0] = 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!
Please register or to comment