|
2537
|
1 ## Copyright (C) 1995, 1996 Kurt Hornik |
|
|
2 ## |
|
|
3 ## This program is free software; you can redistribute it and/or modify |
|
|
4 ## it under the terms of the GNU General Public License as published by |
|
2313
|
5 ## the Free Software Foundation; either version 2, or (at your option) |
|
|
6 ## any later version. |
|
2537
|
7 ## |
|
|
8 ## This program is distributed in the hope that it will be useful, but |
|
2313
|
9 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
2537
|
11 ## General Public License for more details. |
|
|
12 ## |
|
2313
|
13 ## You should have received a copy of the GNU General Public License |
|
2537
|
14 ## along with this file. If not, write to the Free Software Foundation, |
|
|
15 ## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
|
1026
|
16 |
|
2311
|
17 ## usage: gammai (a, x) |
|
|
18 ## |
|
|
19 ## Computes the incomplete gamma function |
|
|
20 ## |
|
2537
|
21 ## gammai(a, x) |
|
|
22 ## = (integral from 0 to x of exp(-t) t^(a-1) dt) / gamma(a). |
|
2311
|
23 ## |
|
|
24 ## If a is scalar, then gammai(a, x) is returned for each element of x |
|
|
25 ## and vice versa. |
|
|
26 ## |
|
|
27 ## If neither a nor x is scalar, the sizes of a and x must agree, and |
|
|
28 ## gammai is applied pointwise. |
|
|
29 |
|
2312
|
30 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> |
|
|
31 ## Created: 13 August 1994 |
|
|
32 ## Adapted-By: jwe |
|
|
33 |
|
1026
|
34 function y = gammai (a, x) |
|
2537
|
35 |
|
715
|
36 if (nargin != 2) |
|
1026
|
37 usage ("gammai (a, x)"); |
|
715
|
38 endif |
|
2537
|
39 |
|
|
40 [retval, a, x] = common_size (a, x); |
|
|
41 if (retval > 0) |
|
|
42 error ("gammai: a and x must be of common size or scalar"); |
|
|
43 endif |
|
|
44 |
|
|
45 [r, c] = size (x); |
|
|
46 s = r * c; |
|
|
47 x = reshape (x, 1, s); |
|
|
48 a = reshape (a, 1, s); |
|
|
49 y = zeros (1, s); |
|
2325
|
50 |
|
2537
|
51 k = find (!(a > 0) | isnan (x)); |
|
|
52 if any (k) |
|
|
53 y(k) = NaN * ones (1, length (k)); |
|
715
|
54 endif |
|
2537
|
55 |
|
|
56 k = find ((x == Inf) & (a > 0)); |
|
|
57 if any (k) |
|
|
58 y(k) = ones (1, length (k)); |
|
715
|
59 endif |
|
2537
|
60 |
|
|
61 ## For x < a + 1, use summation. The below choice of L should ensure |
|
|
62 ## that the overall error is less than eps ... |
|
|
63 k = find((x > 0) & (x < a + 1)); |
|
|
64 if any (k) |
|
|
65 L = ceil (- max ([a(k), x(k)]) * log (eps)); |
|
|
66 A = cumprod ((ones (L, 1) * x(k)) ... |
|
|
67 ./ (ones (L, 1) * a(k) + (1 : L)' * ones (1, length (k)))); |
|
|
68 y(k) = exp (-x(k) + a(k) .* log (x(k))) ... |
|
|
69 .* (1 + sum (A)) ./ gamma (a(k) + 1); |
|
715
|
70 endif |
|
|
71 |
|
2303
|
72 ## For x >= a + 1, use the continued fraction. |
|
|
73 ## Note, however, that this converges MUCH slower than the series |
|
|
74 ## expansion for small a and x not too large! |
|
2537
|
75 k = find ((x >= a + 1) & (x < Inf) & (a > 0)); |
|
|
76 if any (k) |
|
|
77 len = length (k); |
|
|
78 u = [zeros (1, len); ones (1, len)]; |
|
|
79 v = [ones (1, len); x(k)]; |
|
715
|
80 c_old = 0; |
|
2537
|
81 c_new = v(1, :) ./ v(2, :); |
|
715
|
82 n = 1; |
|
1026
|
83 while (max (abs (c_old ./ c_new - 1)) > 10 * eps) |
|
715
|
84 c_old = c_new; |
|
2537
|
85 u = v + u .* (ones (2, 1) * (n - a(k))); |
|
|
86 v = u .* (ones (2, 1) * x(k)) + n * v; |
|
|
87 c_new = v(1, :) ./ v(2, :); |
|
715
|
88 n = n + 1; |
|
|
89 endwhile |
|
2537
|
90 y(k) = 1 - exp (-x(k) + a(k) .* log (x(k))) .* c_new ... |
|
|
91 ./ gamma (a(k)); |
|
715
|
92 endif |
|
2537
|
93 |
|
|
94 y = reshape (y, r, c); |
|
2325
|
95 |
|
2554
|
96 endfunction |
