home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / scripts / specfun / gammai.m < prev    next >
Text File  |  1996-09-28  |  4KB  |  124 lines
  1. # Copyright (C) 1995 John W. Eaton
  2. # 
  3. # This file is part of Octave.
  4. # 
  5. # Octave is free software; you can redistribute it and/or modify it
  6. # under the terms of the GNU General Public License as published by the
  7. # Free Software Foundation; either version 2, or (at your option) any
  8. # later version.
  9. # 
  10. # Octave is distributed in the hope that it will be useful, but WITHOUT
  11. # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12. # FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  13. # for more details.
  14. # 
  15. # You should have received a copy of the GNU General Public License
  16. # along with Octave; see the file COPYING.  If not, write to the Free
  17. # Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  
  19. function y = gammai (a, x)
  20.   
  21. # usage: gammai (a, x)
  22. #
  23. # Computes the incomplete gamma function
  24. #
  25. #   gammai (a, x) 
  26. #     = (integral from 0 to x of exp(-t) t^(a-1) dt) / gamma(a).
  27. #
  28. # If a is scalar, then gammai(a, x) is returned for each element of x
  29. # and vice versa.
  30. #
  31. # If neither a nor x is scalar, the sizes of a and x must agree, and
  32. # gammai is applied pointwise.
  33.   
  34. # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Aug 13, 1994
  35.  
  36.   if (nargin != 2)
  37.     usage ("gammai (a, x)");
  38.   endif
  39.   
  40.   [r_a, c_a] = size (a);
  41.   [r_x, c_x] = size (x);
  42.   e_a = r_a * c_a;
  43.   e_x = r_x * c_x;
  44.   
  45.   # The following code is rather ugly.  We want the function to work
  46.   # whenever a and x have the same size or a or x is scalar.  
  47.   # We do this by reducing the latter cases to the former.
  48.   
  49.   if (e_a == 0 || e_x == 0)
  50.     error ("gammai: both a and x must be nonempty");
  51.   endif
  52.   if (r_a == r_x && c_a == c_x)
  53.     n   = e_a;
  54.     a   = reshape (a, 1, n);
  55.     x   = reshape (x, 1, n);
  56.     r_y = r_a;
  57.     c_y = c_a;
  58.   elseif (e_a == 1)
  59.     n   = e_x;
  60.     a   = a * ones (1, n);
  61.     x   = reshape (x, 1, n);
  62.     r_y = r_x;
  63.     c_y = c_x;
  64.   elseif (e_x == 1)
  65.     n   = e_a;
  66.     a   = reshape (a, 1, n);
  67.     x   = x * ones (1, n);
  68.     r_y = r_a;
  69.     c_y = c_a;
  70.   else
  71.     error ("gammai: a and x must have the same size if neither is scalar"); 
  72.   endif
  73.  
  74. # Now we can do sanity checking ...
  75.   
  76.   if (any (a <= 0) || any (a == Inf))
  77.     error ("gammai: all entries of a must be positive anf finite");
  78.   endif
  79.   if (any (x < 0))
  80.     error ("gammai: all entries of x must be nonnegative");
  81.   endif
  82.   
  83.   y = zeros(1, n);
  84.  
  85. # For x < a + 1, use summation.  The below choice of k should ensure
  86. # that the overall error is less than eps ... 
  87.  
  88.   S = find ((x > 0) & (x < a + 1));
  89.   s = length (S);
  90.   if (s > 0)
  91.     k   = ceil (- max ([a(S), x(S)]) * log (eps));
  92.     K   = (1:k)';
  93.     M   = ones(k, 1);
  94.     A   = cumprod((M * x(S)) ./ (M * a(S) + K * ones(1, s)));
  95.     y(S) = exp (-x(S) + a(S) .* log (x(S))) .* (1 + sum (A)) ./ gamma (a(S)+1);
  96.   endif
  97.  
  98. # For x >= a + 1, use the continued fraction.
  99. # Note, however, that this converges MUCH slower than the series
  100. # expansion for small a and x not too large!
  101.  
  102.   S = find((x >= a + 1) & (x < Inf));
  103.   s = length(S);
  104.   if (s > 0)
  105.     u   = [zeros(1, s); ones(1, s)];
  106.     v   = [ones(1, s); x(S)];
  107.     c_old = 0;
  108.     c_new = v(1,:) ./ v(2,:);
  109.     n   = 1;
  110.     while (max (abs (c_old ./ c_new - 1)) > 10 * eps)
  111.       c_old = c_new;
  112.       u = v + u .* (ones (2, 1) * (n - a(S)));
  113.       v = u .* (ones (2, 1) * x(S)) + n * v;
  114.       c_new = v(1,:) ./ v(2,:);
  115.       n = n + 1;
  116.     endwhile
  117.     y(S) = 1 - exp (-x(S) + a(S) .* log (x(S))) .* c_new ./ gamma (a(S));
  118.   endif
  119.   
  120.   y (find (x == Inf)) = ones (1, sum (x == Inf));
  121.   
  122.   y = reshape (y, r_y, c_y);
  123.  
  124. endfunction