Mercurial > octave
changeset 1026:9fc405c8c06c
[project @ 1995-01-11 21:17:01 by jwe]
| author | jwe |
|---|---|
| date | Wed, 11 Jan 1995 21:17:01 +0000 |
| parents | f558749713f1 |
| children | 005998569449 |
| files | scripts/elfun/gcd.m scripts/elfun/lcm.m scripts/linear-algebra/null.m scripts/linear-algebra/orth.m scripts/signal/fftconv.m scripts/signal/fftfilt.m scripts/specfun/beta.m scripts/specfun/betai.m scripts/specfun/gammai.m scripts/statistics/kurtosis.m scripts/statistics/skewness.m |
| diffstat | 11 files changed, 281 insertions(+), 118 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/elfun/gcd.m Wed Jan 11 20:52:10 1995 +0000 +++ b/scripts/elfun/gcd.m Wed Jan 11 21:17:01 1995 +0000 @@ -1,13 +1,32 @@ +# Copyright (C) 1995 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, 675 Mass Ave, Cambridge, MA 02139, USA. + function [g, v] = gcd (a, ...) - + +# usage: gcd (a, ...) +# # [g [, v]] = gcd (a) returns the greatest common divisor g of the # entries of the integer vector a, and an integer vector v such that # g = v(1) * a(k) + ... + v(k) * a(k). # # [g [, v]] = gcd (a1, ..., ak) is the same with a = [a1, ..., ak]. -# Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Sep 16, 1994 -# Copyright Dept of Statistics and Probability Theory TU Wien +# Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Sep 16, 1994. if (nargin > 1) va_start; @@ -17,7 +36,7 @@ endif if (round (a) != a) - error ("gcd: all arguments must be integer"); + error ("gcd: all arguments must be integer"); endif g = abs (a(1));
--- a/scripts/elfun/lcm.m Wed Jan 11 20:52:10 1995 +0000 +++ b/scripts/elfun/lcm.m Wed Jan 11 21:17:01 1995 +0000 @@ -1,11 +1,30 @@ +# Copyright (C) 1995 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, 675 Mass Ave, Cambridge, MA 02139, USA. + function l = lcm (a, ...) +# usage: lcm (a, ...) +# # lcm (a) returns the least common multiple of the entries of the # integer vector a. # lcm (a1, ..., ak) is the same as lcm([a1, ..., ak]). -# Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Sep 16, 1994 -# Copyright Dept of Statistics and Probability Theory TU Wien +# Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Sep 16, 1994. if (nargin > 1) va_start; @@ -24,7 +43,7 @@ a = abs (a); l = a (1); for k = 1:(length (a) - 1) - l = l * a(k+1) / gcd(l, a(k+1)); + l = l * a(k+1) / gcd (l, a(k+1)); endfor endif
--- a/scripts/linear-algebra/null.m Wed Jan 11 20:52:10 1995 +0000 +++ b/scripts/linear-algebra/null.m Wed Jan 11 21:17:01 1995 +0000 @@ -1,3 +1,21 @@ +# Copyright (C) 1995 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, 675 Mass Ave, Cambridge, MA 02139, USA. + function retval = null (A, tol) # usage: null (A, tol) @@ -10,8 +28,7 @@ # max (size (A)) * sigma_max (A) * eps, where sigma_max (A) is the # maximal singular value of A. -# written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Dec 24, 1993 -# copyright Dept of Probability Theory and Statistics TU Wien +# Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Dec 24, 1993. [U, S, V] = svd (A); @@ -22,7 +39,7 @@ if (nargin == 1) tol = max (size (A)) * s (1) * eps; else (nargin != 2) - usage ("null(A [, tol])"); + usage ("null (A [, tol])"); endif rank = sum (s > tol);
--- a/scripts/linear-algebra/orth.m Wed Jan 11 20:52:10 1995 +0000 +++ b/scripts/linear-algebra/orth.m Wed Jan 11 21:17:01 1995 +0000 @@ -1,3 +1,21 @@ +# Copyright (C) 1995 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, 675 Mass Ave, Cambridge, MA 02139, USA. + function retval = orth (A, tol) # usage: orth (A, tol) @@ -10,8 +28,7 @@ # max (size (A)) * sigma_max (A) * eps, where sigma_max (A) is the # maximal singular value of A. -# written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Dec 24, 1993 -# copyright Dept of Probability Theory and Statistics TU Wien +# Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Dec 24, 1993. [U, S, V] = svd (A);
--- a/scripts/signal/fftconv.m Wed Jan 11 20:52:10 1995 +0000 +++ b/scripts/signal/fftconv.m Wed Jan 11 21:17:01 1995 +0000 @@ -1,6 +1,24 @@ +# Copyright (C) 1995 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, 675 Mass Ave, Cambridge, MA 02139, USA. + function c = fftconv (a, b, N) -# usage: fftconv (a, b [, N]) +# usage: fftconv (a, b [, N]) # # c = fftconv (a, b) returns the convolution of the vectors a and b, # a vector with length equal to length (a) + length (b) - 1. @@ -10,11 +28,10 @@ # The computation uses the FFT by calling fftfilt. If the optional # argument N is specified, an N-point FFT is used. -# Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Sep 3, 1994 -# Copyright Dept of Statistics and Probability Theory TU Wien +# Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Sep 3, 1994. if (nargin < 2 || nargin > 3) - usage (" fftconv (b, x [, N])"); + usage ("fftconv (b, x [, N])"); endif if (is_matrix (a) || is_matrix (b)) @@ -32,7 +49,7 @@ c = fftfilt (a, b); else if !(is_scalar (N)) - error ("fftconv: N has to be a scalar"); + error ("fftconv: N has to be a scalar"); endif c = fftfilt (a, b, N); endif
--- a/scripts/signal/fftfilt.m Wed Jan 11 20:52:10 1995 +0000 +++ b/scripts/signal/fftfilt.m Wed Jan 11 21:17:01 1995 +0000 @@ -1,3 +1,21 @@ +# Copyright (C) 1995 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, 675 Mass Ave, Cambridge, MA 02139, USA. + function y = fftfilt (b, x, N) # usage: fftfilt (b, x [, N]) @@ -7,7 +25,6 @@ # b using an N-point FFT. # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Sep 3, 1994 -# Copyright Dept of Statistics and Probability Theory TU Wien # Reference: Oppenheim & Schafer (1989). Discrete-time Signal # Processing (Chapter 8). Prentice-Hall. @@ -24,8 +41,8 @@ [r_x, c_x] = size (x); [r_b, c_b] = size (b); - if !( (min ([r_x c_x]) == 1) || (min ([r_b c_b]) == 1) ) - error ("fftfilt: both x and b should be vectors."); + if (! (min ([r_x c_x]) == 1 || min ([r_b c_b]) == 1)) + error ("fftfilt: both x and b should be vectors"); endif l_x = r_x * c_x; l_b = r_b * c_b; @@ -35,25 +52,25 @@ return; endif - x = reshape (x, 1, l_x); - b = reshape (b, 1, l_b); + x = reshape (x, 1, l_x); + b = reshape (b, 1, l_b); if (nargin == 2) - # use FFT with the smallest power of 2 which is >= length (x) + - # length (b) - 1 as number of points ... +# Use FFT with the smallest power of 2 which is >= length (x) + +# length (b) - 1 as number of points ... N = 2^(ceil (log (l_x + l_b - 1) / log(2))); y = ifft (fft (x, N) .* fft(b, N)); else - # use overlap-add method ... +# Use overlap-add method ... if !(is_scalar (N)) - error ("fftfilt: N has to be a scalar"); + error ("fftfilt: N has to be a scalar"); endif - N = 2^(ceil (log (max ([N l_b])) / log(2))); - L = N - l_b + 1; - B = fft (b, N); - R = ceil (l_x / L); - y = zeros (1, l_x); - for r=1:R; + N = 2^(ceil (log (max ([N l_b])) / log(2))); + L = N - l_b + 1; + B = fft (b, N); + R = ceil (l_x / L); + y = zeros (1, l_x); + for r = 1:R; lo = (r - 1) * L + 1; hi = min (r * L, l_x); tmp = ifft (fft (x(lo:hi), N) .* B); @@ -62,14 +79,15 @@ endfor endif - y = reshape (y(1:l_x), r_x, c_x); + y = reshape (y(1:l_x), r_x, c_x); - # final cleanups: if both x and b are real respectively integer, y - # should also be so - if !(any (imag (x)) || any (imag (b))) +# Final cleanups: if both x and b are real respectively integer, y +# should also be + + if (! (any (imag (x)) || any (imag (b)))) y = real (y); endif - if !(any (x - round (x)) || any (b - round (b))) + if (! (any (x - round (x)) || any (b - round (b)))) y = round (y); endif
--- a/scripts/specfun/beta.m Wed Jan 11 20:52:10 1995 +0000 +++ b/scripts/specfun/beta.m Wed Jan 11 21:17:01 1995 +0000 @@ -1,15 +1,35 @@ -function retval = beta(a, b) +# Copyright (C) 1995 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, 675 Mass Ave, Cambridge, MA 02139, USA. + +function retval = beta (a, b) -# usage: beta(a, b) +# usage: beta (a, b) # # Returns the beta function beta(a,b) = gamma(a) * gamma(b) / gamma(a+b) # of a and b. # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Jun 13, 1993 -# Updated by KH (Kurt.Hornik@ci.tuwien.ac.at) on Aug 13, 1994 -# Copyright Dept of Probability Theory and Statistics TU Wien - retval = exp(lgamma(a) + lgamma(b) - lgamma(a+b)); + if (nargin != 2) + usage ("beta (a, b)"); + endif + + retval = exp (lgamma (a) + lgamma (b) - lgamma (a+b)); endfunction
--- a/scripts/specfun/betai.m Wed Jan 11 20:52:10 1995 +0000 +++ b/scripts/specfun/betai.m Wed Jan 11 21:17:01 1995 +0000 @@ -1,6 +1,24 @@ +# Copyright (C) 1995 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, 675 Mass Ave, Cambridge, MA 02139, USA. + 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. @@ -8,7 +26,6 @@ # If x is a scalar, a and b must be of compatible dimensions. # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Aug 2, 1994. -# Copyright Dept of Probability Theory and Statistics TU Wien. # Computation is based on the series expansion # betai(a, b, x) @@ -24,78 +41,79 @@ error("betai: a and b must both be positive."); endif [nr, nc] = size(x); - if (min([nr nc]) == 0) + if (min ([nr, nc]) == 0) 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]."); + 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 !(is_scalar(a) && is_scalar(b)) - error("betai: if x is not a scalar, a and b must be scalars."); + if (! (is_scalar (a) && is_scalar (b))) + error ("betai: if x is not a scalar, a and b must be scalars"); endif - n = nr * nc; - x = reshape(x, 1, n); - y = zeros(1, n); - c = exp(lgamma(a+b) - lgamma(a) - lgamma(b)); + n = nr * nc; + x = reshape (x, 1, n); + y = zeros (1, n); + c = exp (lgamma (a+b) - lgamma (a) - lgamma (b)); - y(find(x == 1)) = ones(1, sum(x == 1)); + y (find (x == 1)) = ones (1, sum (x == 1)); - # now do the series computation. The error when truncating at term K - # is always less than 2^(-K), hence the following choice of K. - K = ceil(-log(eps) / log(2)); - k = (1:K)'; +# Now do the series computation. The error when truncating at term K +# is always less than 2^(-K), hence the following choice of K. - ind = find((x > 0) & (x <= 1/2)); - len = length(ind); + K = ceil (-log (eps) / log (2)); + k = (1:K)'; + + ind = find ((x > 0) & (x <= 1/2)); + len = length (ind); if (len > 0) tmp = cumprod((1 - b./k) * x(ind)) ./ ((a + k) * ones(1, len)); y(ind) = c * exp(a * log(x(ind))) .* (1/a + sum(tmp)); endif - ind = find((x > 1/2) & (x < 1)); + ind = find ((x > 1/2) & (x < 1)); len = length(ind); if (len > 0) - tmp = cumprod((1 - a./k) * (1 - x(ind))) ./ ((b + k) * ones(1, len)); - y(ind) = 1 - c * exp(b * log(1-x(ind))) .* (1/b + sum(tmp)); + tmp = cumprod ((1 - a./k) * (1 - x(ind))) ./ ((b + k) * ones(1, len)); + y(ind) = 1 - c * exp (b * log (1-x(ind))) .* (1/b + sum (tmp)); endif - y = reshape(y, nr, nc); + y = reshape (y, nr, nc); else - [ra, ca] = size(a); - [rb, cb] = size(b); - if !((ra == rb) && (ca == cb)) - error("betai: a and b must have the same size."); + [ra, ca] = size (a); + [rb, cb] = size (b); + if (! (ra == rb && ca == cb)) + error ("betai: a and b must have the same size"); endif - - n = ra * ca; - a = reshape(a, 1, n); - b = reshape(b, 1, n); - c = exp(lgamma(a+b) - lgamma(a) - lgamma(b)); + + n = ra * ca; + a = reshape (a, 1, n); + b = reshape (b, 1, n); + c = exp (lgamma (a+b) - lgamma (a) - lgamma (b)); if (x == 0) - y = zeros(1, n); + y = zeros (1, n); elseif (x == 1) - y = ones(1, n); + y = ones (1, n); else - K = ceil(-log(eps) / log(2)); - k = (1:K)' * ones(1, n); - h = ones(K, 1); - if ((x > 0) && (x <= 1/2)) - tmp = cumprod((1 - (h * b) ./ k) * x) ./ ((h * a) + k); - y = c .* exp(a * log(x)) .* (1 ./ a + sum(tmp)); + K = ceil (-log (eps) / log (2)); + k = (1:K)' * ones (1, n); + h = ones (K, 1); + if (x > 0 && x <= 1/2) + tmp = cumprod ((1 - (h * b) ./ k) * x) ./ ((h * a) + k); + y = c .* exp (a * log(x)) .* (1 ./ a + sum (tmp)); else - tmp = cumprod((1 - (h * a) ./ k) * (1-x)) ./ ((h * b) + k); - y = 1 - c .* exp(b * log(1-x)) .* (1 ./ b + sum(tmp)); + tmp = cumprod ((1 - (h * a) ./ k) * (1-x)) ./ ((h * b) + k); + y = 1 - c .* exp (b * log (1-x)) .* (1 ./ b + sum (tmp)); endif endif - y = reshape(y, ra, ca); + y = reshape (y, ra, ca); endif
--- a/scripts/specfun/gammai.m Wed Jan 11 20:52:10 1995 +0000 +++ b/scripts/specfun/gammai.m Wed Jan 11 21:17:01 1995 +0000 @@ -1,24 +1,44 @@ -function y = gammai(a, x) +# Copyright (C) 1995 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, 675 Mass Ave, Cambridge, MA 02139, USA. + +function y = gammai (a, x) -# usage: gammai(a, x) +# 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). +# +# 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. # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Aug 13, 1994 -# Copyright Dept of Probability Theory and Statistics TU Wien if (nargin != 2) - usage (" gammai(a, x)"); + usage ("gammai (a, x)"); endif - [r_a, c_a] = size(a); - [r_x, c_x] = size(x); + [r_a, c_a] = size (a); + [r_x, c_x] = size (x); e_a = r_a * c_a; e_x = r_x * c_x; @@ -26,38 +46,38 @@ # whenever a and x have the same size or a or x is scalar. # We do this by reducing the latter cases to the former. - if ((e_a == 0) || (e_x == 0)) - error("gammai: both a and x must be nonempty"); + if (e_a == 0 || e_x == 0) + error ("gammai: both a and x must be nonempty"); endif - if ((r_a == r_x) && (c_a == c_x)) + if (r_a == r_x && c_a == c_x) n = e_a; - a = reshape(a, 1, n); - x = reshape(x, 1, n); + a = reshape (a, 1, n); + x = reshape (x, 1, n); r_y = r_a; c_y = c_a; elseif (e_a == 1) n = e_x; - a = a * ones(1, n); - x = reshape(x, 1, n); + a = a * ones (1, n); + x = reshape (x, 1, n); r_y = r_x; c_y = c_x; elseif (e_x == 1) n = e_a; - a = reshape(a, 1, n); - x = x * ones(1, n); + a = reshape (a, 1, n); + x = x * ones (1, n); r_y = r_a; c_y = c_a; else - error("gammai: a and x must have the same size if neither is scalar"); + error ("gammai: a and x must have the same size if neither is scalar"); endif - # Now we can do sanity checking ... +# Now we can do sanity checking ... if (any (a <= 0) || any (a == Inf)) - error ("gammai: all entries of a must be positive anf finite"); + error ("gammai: all entries of a must be positive anf finite"); endif if (any (x < 0)) - error ("gammai: all entries of x must be nonnegative"); + error ("gammai: all entries of x must be nonnegative"); endif y = zeros(1, n); @@ -65,14 +85,14 @@ # For x < a + 1, use summation. The below choice of k should ensure # that the overall error is less than eps ... - S = find((x > 0) & (x < a + 1)); - s = length(S); + S = find ((x > 0) & (x < a + 1)); + s = length (S); if (s > 0) - k = ceil(- max([a(S), x(S)]) * log(eps)); + 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))); - y(S) = exp(-x(S) + a(S) .* log(x(S))) .* (1 + sum(A)) ./ gamma(a(S)+1); + y(S) = exp (-x(S) + a(S) .* log (x(S))) .* (1 + sum (A)) ./ gamma (a(S)+1); endif # For x >= a + 1, use the continued fraction. @@ -87,18 +107,18 @@ c_old = 0; c_new = v(1,:) ./ v(2,:); n = 1; - while (max(abs(c_old ./ c_new - 1)) > 10 * eps) + while (max (abs (c_old ./ c_new - 1)) > 10 * eps) c_old = c_new; - u = v + u .* (ones(2, 1) * (n - a(S))); - v = u .* (ones(2, 1) * x(S)) + n * v; + u = v + u .* (ones (2, 1) * (n - a(S))); + v = u .* (ones (2, 1) * x(S)) + n * v; c_new = v(1,:) ./ v(2,:); n = n + 1; endwhile - y(S) = 1 - exp(-x(S) + a(S) .* log(x(S))) .* c_new ./ gamma(a(S)); + y(S) = 1 - exp (-x(S) + a(S) .* log (x(S))) .* c_new ./ gamma (a(S)); endif - y(find(x == Inf)) = ones(1, sum(x == Inf)); + y (find (x == Inf)) = ones (1, sum (x == Inf)); - y = reshape(y, r_y, c_y); + y = reshape (y, r_y, c_y); endfunction \ No newline at end of file
--- a/scripts/statistics/kurtosis.m Wed Jan 11 20:52:10 1995 +0000 +++ b/scripts/statistics/kurtosis.m Wed Jan 11 21:17:01 1995 +0000 @@ -22,7 +22,7 @@ # # If x is a vector of length N, return the kurtosis # -# kurtosis(x) = N^(-1) std(x)^(-4) SUM_i (x(i)-mean(x))^4 - 3 +# kurtosis(x) = N^(-1) std(x)^(-4) SUM_i (x(i)-mean(x))^4 - 3 # # of x. # @@ -30,7 +30,6 @@ # column. # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Jul 29, 1994. -# Copyright Dept of Probability Theory and Statistics TU Wien, Austria. if (nargin != 1) usage ("kurtosis (x)");
--- a/scripts/statistics/skewness.m Wed Jan 11 20:52:10 1995 +0000 +++ b/scripts/statistics/skewness.m Wed Jan 11 21:17:01 1995 +0000 @@ -30,7 +30,6 @@ # column. # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Jul 29, 1994. -# Copyright Dept of Probability Theory and Statistics TU Wien, Austria. if (nargin != 1) usage ("skewness (x)");