Mercurial > hg > octave-terminal
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
deleted file mode 100644 --- a/scripts/specfun/betai.m +++ /dev/null @@ -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 - - - - -
deleted file mode 100644 --- a/scripts/specfun/betainc.m +++ /dev/null @@ -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
deleted file mode 100644 --- a/scripts/specfun/gammai.m +++ /dev/null @@ -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
deleted file mode 100644 --- a/scripts/specfun/gammainc.m +++ /dev/null @@ -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