|
715
|
1 function y = gammai(a, x) |
|
|
2 |
|
904
|
3 # usage: gammai(a, x) |
|
|
4 # |
|
|
5 # Computes the incomplete gamma function |
|
|
6 # gammai(a, x) |
|
|
7 # = (integral from 0 to x of exp(-t) t^(a-1) dt) / gamma(a). |
|
|
8 # If a is scalar, then gammai(a, x) is returned for each element of x |
|
|
9 # and vice versa. |
|
|
10 # If neither a nor x is scalar, the sizes of a and x must agree, and |
|
|
11 # gammai is applied pointwise. |
|
715
|
12 |
|
904
|
13 # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Aug 13, 1994 |
|
|
14 # Copyright Dept of Probability Theory and Statistics TU Wien |
|
715
|
15 |
|
|
16 if (nargin != 2) |
|
904
|
17 usage (" gammai(a, x)"); |
|
715
|
18 endif |
|
|
19 |
|
|
20 [r_a, c_a] = size(a); |
|
|
21 [r_x, c_x] = size(x); |
|
|
22 e_a = r_a * c_a; |
|
|
23 e_x = r_x * c_x; |
|
|
24 |
|
|
25 # The following code is rather ugly. We want the function to work |
|
|
26 # whenever a and x have the same size or a or x is scalar. |
|
|
27 # We do this by reducing the latter cases to the former. |
|
|
28 |
|
|
29 if ((e_a == 0) || (e_x == 0)) |
|
|
30 error("gammai: both a and x must be nonempty"); |
|
|
31 endif |
|
|
32 if ((r_a == r_x) && (c_a == c_x)) |
|
|
33 n = e_a; |
|
|
34 a = reshape(a, 1, n); |
|
|
35 x = reshape(x, 1, n); |
|
|
36 r_y = r_a; |
|
|
37 c_y = c_a; |
|
|
38 elseif (e_a == 1) |
|
|
39 n = e_x; |
|
|
40 a = a * ones(1, n); |
|
|
41 x = reshape(x, 1, n); |
|
|
42 r_y = r_x; |
|
|
43 c_y = c_x; |
|
|
44 elseif (e_x == 1) |
|
|
45 n = e_a; |
|
|
46 a = reshape(a, 1, n); |
|
|
47 x = x * ones(1, n); |
|
|
48 r_y = r_a; |
|
|
49 c_y = c_a; |
|
|
50 else |
|
|
51 error("gammai: a and x must have the same size if neither is scalar"); |
|
|
52 endif |
|
|
53 |
|
|
54 # Now we can do sanity checking ... |
|
|
55 |
|
904
|
56 if (any (a <= 0) || any (a == Inf)) |
|
|
57 error ("gammai: all entries of a must be positive anf finite"); |
|
715
|
58 endif |
|
904
|
59 if (any (x < 0)) |
|
|
60 error ("gammai: all entries of x must be nonnegative"); |
|
715
|
61 endif |
|
|
62 |
|
|
63 y = zeros(1, n); |
|
|
64 |
|
904
|
65 # For x < a + 1, use summation. The below choice of k should ensure |
|
|
66 # that the overall error is less than eps ... |
|
|
67 |
|
715
|
68 S = find((x > 0) & (x < a + 1)); |
|
|
69 s = length(S); |
|
|
70 if (s > 0) |
|
|
71 k = ceil(- max([a(S), x(S)]) * log(eps)); |
|
|
72 K = (1:k)'; |
|
|
73 M = ones(k, 1); |
|
|
74 A = cumprod((M * x(S)) ./ (M * a(S) + K * ones(1, s))); |
|
|
75 y(S) = exp(-x(S) + a(S) .* log(x(S))) .* (1 + sum(A)) ./ gamma(a(S)+1); |
|
|
76 endif |
|
|
77 |
|
904
|
78 # For x >= a + 1, use the continued fraction. |
|
|
79 # Note, however, that this converges MUCH slower than the series |
|
|
80 # expansion for small a and x not too large! |
|
|
81 |
|
715
|
82 S = find((x >= a + 1) & (x < Inf)); |
|
|
83 s = length(S); |
|
|
84 if (s > 0) |
|
|
85 u = [zeros(1, s); ones(1, s)]; |
|
|
86 v = [ones(1, s); x(S)]; |
|
|
87 c_old = 0; |
|
|
88 c_new = v(1,:) ./ v(2,:); |
|
|
89 n = 1; |
|
|
90 while (max(abs(c_old ./ c_new - 1)) > 10 * eps) |
|
|
91 c_old = c_new; |
|
|
92 u = v + u .* (ones(2, 1) * (n - a(S))); |
|
|
93 v = u .* (ones(2, 1) * x(S)) + n * v; |
|
|
94 c_new = v(1,:) ./ v(2,:); |
|
|
95 n = n + 1; |
|
|
96 endwhile |
|
|
97 y(S) = 1 - exp(-x(S) + a(S) .* log(x(S))) .* c_new ./ gamma(a(S)); |
|
|
98 endif |
|
|
99 |
|
|
100 y(find(x == Inf)) = ones(1, sum(x == Inf)); |
|
|
101 |
|
|
102 y = reshape(y, r_y, c_y); |
|
|
103 |
|
|
104 endfunction |
