Mercurial > hg > octave-avbm
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 +++ b/scripts/specfun/betai.m @@ -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 +++ b/scripts/specfun/gammai.m @@ -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;