changeset 1586:54350c98055c

[project @ 1995-10-31 05:22:53 by jwe]
author jwe
date Tue, 31 Oct 1995 05:24:40 +0000
parents 100413a7e8a2
children dd087a402811
files scripts/specfun/betai.m scripts/specfun/gammai.m
diffstat 2 files changed, 16 insertions(+), 14 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/specfun/betai.m	Tue Oct 31 05:14:17 1995 +0000
+++ b/scripts/specfun/betai.m	Tue Oct 31 05:24:40 1995 +0000
@@ -16,9 +16,9 @@
 # along with Octave; see the file COPYING.  If not, write to the Free
 # Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 
-function y = betai(a, b, x)
+function y = betai (a, b, x)
   
-# usage: betai(a, b, x)
+# usage: betai (a, b, x)
 #
 # Returns the incomplete beta function
 #   betai (a, b, x) = BETA(a,b)^(-1) INT_0^x t^(a-1) (1-t)^(b-1) dt.
@@ -37,18 +37,18 @@
     usage (" betai (a, b, x)");
   endif
 
-  if !((a > 0) && (b > 0))
-    error("betai:  a and b must both be positive.");
+  if (! (a > 0 && b > 0))
+    error ("betai: a and b must both be positive");
   endif
-  [nr, nc] = size(x);
+  [nr, nc] = size (x);
   if (min ([nr, nc]) == 0)
-    error ("betai:  x must not be empty.");
+    error ("betai: x must not be empty.");
   endif
   if (any (x < 0) || any (x > 1))
     error ("betai: all entries of x must be in [0,1].");
   endif
 
-  if ((nr > 1) || (nc > 1))
+  if (nr > 1 || nc > 1)
     
     if (! (is_scalar (a) && is_scalar (b)))
       error ("betai: if x is not a scalar, a and b must be scalars");
--- a/scripts/specfun/gammai.m	Tue Oct 31 05:14:17 1995 +0000
+++ b/scripts/specfun/gammai.m	Tue Oct 31 05:24:40 1995 +0000
@@ -80,7 +80,7 @@
     error ("gammai: all entries of x must be nonnegative");
   endif
   
-  y = zeros(1, n);
+  y = zeros (1, n);
 
 # For x < a + 1, use summation.  The below choice of k should ensure
 # that the overall error is less than eps ... 
@@ -90,8 +90,8 @@
   if (s > 0)
     k   = ceil (- max ([a(S), x(S)]) * log (eps));
     K   = (1:k)';
-    M   = ones(k, 1);
-    A   = cumprod((M * x(S)) ./ (M * a(S) + K * ones(1, s)));
+    M   = ones (k, 1);
+    A   = cumprod ((M * x(S)) ./ (M * a(S) + K * ones(1, s)));
     y(S) = exp (-x(S) + a(S) .* log (x(S))) .* (1 + sum (A)) ./ gamma (a(S)+1);
   endif
 
@@ -99,11 +99,13 @@
 # Note, however, that this converges MUCH slower than the series
 # expansion for small a and x not too large!
 
-  S = find((x >= a + 1) & (x < Inf));
-  s = length(S);
+  S = find ((x >= a + 1) & (x < Inf));
+  s = length (S);
   if (s > 0)
-    u   = [zeros(1, s); ones(1, s)];
-    v   = [ones(1, s); x(S)];
+    t1 = zeros (1, s);
+    t2 = ones (1, s);
+    u   = [t1; t2];
+    v   = [t2; x(S)];
     c_old = 0;
     c_new = v(1,:) ./ v(2,:);
     n   = 1;