changeset 3123:e3fc19fa9e69

[project @ 1997-12-10 06:15:37 by jwe]
author jwe
date Wed, 10 Dec 1997 06:15:37 +0000
parents c2d111b3f1bf
children 38684be52a3e
files scripts/specfun/betai.m scripts/specfun/betainc.m scripts/specfun/gammai.m scripts/specfun/gammainc.m
diffstat 4 files changed, 0 insertions(+), 251 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/specfun/betai.m	Mon Dec 01 06:52:23 1997 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,93 +0,0 @@
-## Copyright (C) 1995, 1996  Kurt Hornik
-## 
-## This program is free software; you can redistribute it and/or modify
-## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2, or (at your option)
-## any later version.
-## 
-## This program is distributed in the hope that it will be useful, but
-## WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-## General Public License for more details. 
-## 
-## You should have received a copy of the GNU General Public License
-## along with this file.  If not, write to the Free Software Foundation,
-## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
-
-## 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.
-
-## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
-## Created: 2 August 1994
-## Adapted-By: jwe
-
-## Computation is based on the series expansion
-##   betai(a, b, x) 
-##     = \frac{x^a}{B(a, b)}
-##         \sum_{l=0}^\infty \frac{(1-b)\cdots(l-b)}{a+l} \frac{x^l}{l!}
-## for x <= 1/2.  For x > 1/2, betai(a, b, x) = 1 - betai(b, a, 1-x).
-
-function y = betai (a, b, x)
-  
-  if (nargin != 3)
-    usage ("betai (a, b, x)");
-  endif
-
-  [retval, x, a, b] = common_size (x, a, b);
-  if (retval > 0)
-    error ("betai:  a, b and x must be of common size or scalar");
-  endif
-  
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  a = reshape (a, 1, s);
-  b = reshape (b, 1, s);
-  y = zeros (1, s);
-
-  k = find ((x < 0) | (x > 1) | !(a > 0) | !(b > 0) | isnan (x));
-  if any (k)
-    y(k) = NaN * ones (1, length (k));
-  endif
-
-  k = find ((x == 1) & (a > 0) & (b > 0));
-  if any (k)
-    y(k) = ones (1, length (k));
-  endif
-
-  ## Now do the series computations.  
-  ## The error when truncating at term L is always less than 2^(-L),
-  ## hence the following choice of L. 
-
-  L = ceil (-log (eps) / log (2));
-  h = ones (L, 1);
-  
-  k = find ((x >= 0) & (x <= 1/2) & (a > 0) & (b > 0));
-  if any (k)
-    l   = (1 : L)' * ones (1, length (k));
-    tmp = cumprod ((1 - (h * b(k)) ./ l) .* (h * x(k))) ...
-	./ ((h * a(k)) + l);
-    y(k) = exp (a(k) .* log (x(k))) .* (1 ./ a(k) + sum (tmp)) ...
-	./ beta (a(k), b(k));
-  endif
-  
-  k = find ((x > 1/2) & (x < 1) & (a > 0) & (b > 0));
-  if any (k)
-    l   = (1 : L)' * ones (1, length (k));
-    tmp = cumprod ((1 - (h * a(k)) ./ l) .* (h * (1 - x(k)))) ...
-	./ (h * b(k) + l);
-    y(k) = 1 - exp (b(k) .* log (1 - x(k))) ...
-	.* (1 ./ b(k) + sum (tmp)) ./ beta (a(k), b(k));
-  endif
-
-  y = reshape (y, r, c);
-  
-endfunction
-
-
-
-
-
--- a/scripts/specfun/betainc.m	Mon Dec 01 06:52:23 1997 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,30 +0,0 @@
-## Copyright (C) 1996, 1997 John W. Eaton
-##
-## This file is part of Octave.
-##
-## Octave is free software; you can redistribute it and/or modify it
-## under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2, or (at your option)
-## any later version.
-##
-## Octave is distributed in the hope that it will be useful, but
-## WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-## General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## 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.
-
-## Usage: betainc (x, a, b)
-##
-## See also: betai
-
-## Author: jwe
-
-function y = betainc (x, a, b)
-
-  y = betai (a, b, x);
-
-endfunction
--- a/scripts/specfun/gammai.m	Mon Dec 01 06:52:23 1997 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,98 +0,0 @@
-## Copyright (C) 1995, 1996  Kurt Hornik
-## 
-## This program is free software; you can redistribute it and/or modify
-## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2, or (at your option)
-## any later version.
-## 
-## This program is distributed in the hope that it will be useful, but
-## WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-## General Public License for more details. 
-## 
-## You should have received a copy of the GNU General Public License
-## along with this file.  If not, write to the Free Software Foundation,
-## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
-
-## usage: gammai (a, x)
-##
-## Computes the incomplete gamma function
-##
-##    gammai(a, x) 
-##      = (integral from 0 to x of exp(-t) t^(a-1) dt) / gamma(a).
-##
-## If a is scalar, then gammai(a, x) is returned for each element of x
-## and vice versa.
-##
-## If neither a nor x is scalar, the sizes of a and x must agree, and
-## gammai is applied pointwise.
-
-## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
-## Created: 13 August 1994
-## Adapted-By: jwe
-
-function y = gammai (a, x)
-  
-  if (nargin != 2)
-    usage ("gammai (a, x)");
-  endif
-  
-  [retval, a, x] = common_size (a, x);
-  if (retval > 0)
-    error ("gammai:  a and x must be of common size or scalar");
-  endif
-  
-  [r, c] = size (x);
-  s = r * c;
-  x = reshape (x, 1, s);
-  a = reshape (a, 1, s);
-  y = zeros (1, s);
-
-  k = find (!(a > 0) | isnan (x));
-  if any (k)
-    y(k) = NaN * ones (1, length (k));
-  endif
-  
-  k = find ((x == Inf) & (a > 0));
-  if any (k)
-    y(k) = ones (1, length (k));
-  endif
-  
-  ## For x < a + 1, use summation.  The below choice of L should ensure
-  ## that the overall error is less than eps ... 
-  k = find((x > 0) & (x < a + 1));
-  if any (k)
-    L = ceil (- max ([a(k), x(k)]) * log (eps));
-    A = cumprod ((ones (L, 1) * x(k)) ...
-	./ (ones (L, 1) * a(k) + (1 : L)' * ones (1, length (k))));
-    y(k) = exp (-x(k) + a(k) .* log (x(k))) ...
-	.* (1 + sum (A)) ./ gamma (a(k) + 1);
-  endif
-
-  ## For x >= a + 1, use the continued fraction.
-  ## Note, however, that this converges MUCH slower than the series
-  ## expansion for small a and x not too large!
-  k = find ((x >= a + 1) & (x < Inf) & (a > 0));
-  if any (k)
-    len = length (k);
-    t0  = zeros (1, len);
-    t1  = ones (1, len);
-    u   = [t0; t1];
-    v   = [t1; x(k)];
-    c_old = 0;
-    c_new = v(1, :) ./ v(2, :);
-    n   = 1;
-    while (max (abs (c_old ./ c_new - 1)) > 10 * eps)
-      c_old = c_new;
-      u = v + u .* (ones (2, 1) * (n - a(k)));
-      v = u .* (ones (2, 1) * x(k)) + n * v;
-      c_new = v(1, :) ./ v(2, :);
-      n = n + 1;
-    endwhile
-    y(k) = 1 - exp (-x(k) + a(k) .* log (x(k))) .* c_new ...
-	./ gamma (a(k));
-  endif
-  
-  y = reshape (y, r, c);
-
-endfunction
--- a/scripts/specfun/gammainc.m	Mon Dec 01 06:52:23 1997 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,30 +0,0 @@
-## Copyright (C) 1996, 1997 John W. Eaton
-##
-## This file is part of Octave.
-##
-## Octave is free software; you can redistribute it and/or modify it
-## under the terms of the GNU General Public License as published by
-## the Free Software Foundation; either version 2, or (at your option)
-## any later version.
-##
-## Octave is distributed in the hope that it will be useful, but
-## WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-## General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## 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.
-
-## Usage: gammainc (x, a)
-##
-## See also: gammai
-
-## Author: jwe
-
-function y = gammainc (x, a)
-
-  y = gammai (a, x);
-
-endfunction