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 / betai.m < prev    next >
Text File  |  1996-09-28  |  4KB  |  126 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 = betai(a, b, x)
  20.   
  21. # usage: betai(a, b, x)
  22. #
  23. # Returns the incomplete beta function
  24. #   betai (a, b, x) = BETA(a,b)^(-1) INT_0^x t^(a-1) (1-t)^(b-1) dt.
  25. # If x has more than one component, both a and b must be scalars.
  26. # If x is a scalar, a and b must be of compatible dimensions.
  27.   
  28. # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Aug 2, 1994.
  29.  
  30. # Computation is based on the series expansion
  31. #   betai(a, b, x) 
  32. #     = \frac{1}{B(a, b)} x^a 
  33. #         \sum_{k=0}^\infty \frac{(1-b)\cdots(k-b)}{a+k} \frac{x^k}{k!}
  34. # for x <= 1/2.  For x > 1/2, betai(a, b, x) = 1 - betai(b, a, 1-x).
  35.  
  36.   if (nargin <> 3)
  37.     usage (" betai (a, b, x)");
  38.   endif
  39.  
  40.   if !((a > 0) && (b > 0))
  41.     error("betai:  a and b must both be positive.");
  42.   endif
  43.   [nr, nc] = size(x);
  44.   if (min ([nr, nc]) == 0)
  45.     error ("betai:  x must not be empty.");
  46.   endif
  47.   if (any (x < 0) || any (x > 1))
  48.     error ("betai: all entries of x must be in [0,1].");
  49.   endif
  50.  
  51.   if ((nr > 1) || (nc > 1))
  52.     
  53.     if (! (is_scalar (a) && is_scalar (b)))
  54.       error ("betai: if x is not a scalar, a and b must be scalars");
  55.     endif
  56.  
  57.     n = nr * nc;
  58.     x = reshape (x, 1, n);
  59.     y = zeros (1, n);
  60.     c = exp (lgamma (a+b) - lgamma (a) - lgamma (b));
  61.  
  62.     y (find (x == 1)) = ones (1, sum (x == 1));
  63.  
  64. # Now do the series computation.  The error when truncating at term K
  65. # is always less than 2^(-K), hence the following choice of K.
  66.  
  67.     K = ceil (-log (eps) / log (2));
  68.     k = (1:K)';
  69.  
  70.     ind = find ((x > 0) & (x <= 1/2));
  71.     len = length (ind);
  72.     if (len > 0)
  73.       tmp    = cumprod((1 - b./k) * x(ind)) ./ ((a + k) * ones(1, len));
  74.       y(ind) = c * exp(a * log(x(ind))) .* (1/a + sum(tmp));
  75.     endif
  76.     
  77.     ind = find ((x > 1/2) & (x < 1));
  78.     len = length(ind);
  79.     if (len > 0)
  80.       tmp    = cumprod ((1 - a./k) * (1 - x(ind))) ./ ((b + k) * ones(1, len));
  81.       y(ind) = 1 - c * exp (b * log (1-x(ind))) .* (1/b + sum (tmp));
  82.     endif
  83.   
  84.     y = reshape (y, nr, nc);
  85.     
  86.   else
  87.     
  88.     [ra, ca] = size (a);
  89.     [rb, cb] = size (b);
  90.     if (! (ra == rb && ca == cb))
  91.       error ("betai: a and b must have the same size");
  92.     endif
  93.  
  94.     n = ra * ca;
  95.     a = reshape (a, 1, n);
  96.     b = reshape (b, 1, n);
  97.     c = exp (lgamma (a+b) - lgamma (a) - lgamma (b));
  98.     
  99.     if (x == 0)
  100.       y   = zeros (1, n);
  101.     elseif (x == 1)
  102.       y   = ones (1, n);
  103.     else
  104.       K = ceil (-log (eps) / log (2));
  105.       k = (1:K)' * ones (1, n);
  106.       h = ones (K, 1);
  107.       if (x > 0 && x <= 1/2)
  108.     tmp = cumprod ((1 - (h * b) ./ k) * x) ./ ((h * a) + k);
  109.     y   = c .* exp (a * log(x)) .* (1 ./ a + sum (tmp));
  110.       else
  111.     tmp = cumprod ((1 - (h * a) ./ k) * (1-x)) ./ ((h * b) + k);
  112.     y   = 1 - c .* exp (b * log (1-x)) .* (1 ./ b + sum (tmp));
  113.       endif
  114.     endif
  115.   
  116.     y = reshape (y, ra, ca);
  117.     
  118.   endif
  119.  
  120. endfunction
  121.  
  122.  
  123.  
  124.  
  125.  
  126.