Mercurial > octave-nkf
changeset 4859:265d566cc770
[project @ 2004-04-08 23:52:45 by jwe]
line wrap: on
line diff
--- a/scripts/statistics/distributions/discrete_cdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/discrete_cdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -33,7 +33,7 @@ usage ("discrete_cdf (x, v, p)"); endif - [r, c] = size (x); + sz = size (x); if (! isvector (v)) error ("discrete_cdf: v must be a vector"); @@ -43,16 +43,16 @@ error ("discrete_cdf: p must be a nonzero, nonnegative vector"); endif - n = r * c; + n = numel (x); m = length (v); x = reshape (x, n, 1); v = reshape (v, 1, m); p = reshape (p / sum (p), m, 1); - cdf = zeros (n, 1); + cdf = zeros (sz); k = find (isnan (x)); if (any (k)) - cdf (k) = NaN * ones (length (k), 1); + cdf (k) = NaN; endif k = find (!isnan (x)); if (any (k)) @@ -60,6 +60,4 @@ cdf (k) = ((x(k) * ones (1, m)) >= (ones (n, 1) * v)) * p; endif - cdf = reshape (cdf, r, c); - endfunction
--- a/scripts/statistics/distributions/discrete_inv.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/discrete_inv.m Thu Apr 08 23:52:45 2004 +0000 @@ -33,7 +33,7 @@ usage ("discrete_inv (x, v, p)"); endif - [r, c] = size (x); + sz = size (x); if (! isvector (v)) error ("discrete_inv: v must be a vector"); @@ -43,7 +43,7 @@ error ("discrete_inv: p must be a nonzero, nonnegative vector"); endif - n = r * c; + n = numel (x); x = reshape (x, 1, n); m = length (v); v = sort (v); @@ -52,12 +52,12 @@ ## Allow storage allocated for P to be reclaimed. p = []; - inv = NaN * ones (n, 1); + inv = NaN * ones (sz); if (any (k = find (x == 0))) - inv(k) = -Inf * ones (1, length (k)); + inv(k) = -Inf; endif if (any (k = find (x == 1))) - inv(k) = v(m) * ones (1, length (k)); + inv(k) = v(m) * ones (size (k)); endif if (any (k = find ((x > 0) & (x < 1)))) @@ -68,15 +68,13 @@ ## ## Vectorized code is: ## - ## inv(k) = v(sum ((ones (m, 1) * x(k)) > (s * ones (1, n))) + 1); + ## inv(k) = v(sum ((ones (m, 1) * x(k)) > (s * ones (1, n))) + 1); for q = 1:n - inv(q) = v(sum (x(q) > s) + 1); + inv(k(q)) = v(sum (x(k(q)) > s) + 1); endfor endif - inv = reshape (inv, r, c); - endfunction
--- a/scripts/statistics/distributions/discrete_pdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/discrete_pdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -33,7 +33,7 @@ usage ("discrete_pdf (x, v, p)"); endif - [r, c] = size (x); + sz = size (x); if (! isvector (v)) error ("discrete_pdf: v must be a vector"); @@ -43,16 +43,16 @@ error ("discrete_pdf: p must be a nonzero, nonnegative vector"); endif - n = r * c; + n = numel (x); m = length (v); x = reshape (x, n, 1); v = reshape (v, 1, m); p = reshape (p / sum (p), m, 1); - pdf = zeros (n, 1); + pdf = zeros (sz); k = find (isnan (x)); if (any (k)) - pdf (k) = NaN * ones (length (k), 1); + pdf (k) = NaN; endif k = find (!isnan (x)); if (any (k)) @@ -60,6 +60,4 @@ pdf (k) = ((x(k) * ones (1, m)) == (ones (n, 1) * v)) * p; endif - pdf = reshape (pdf, r, c); - endfunction
--- a/scripts/statistics/distributions/discrete_rnd.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/discrete_rnd.m Thu Apr 08 23:52:45 2004 +0000 @@ -19,24 +19,51 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} discrete_rnd (@var{n}, @var{v}, @var{p}) +## @deftypefnx {Function File} {} discrete_rnd (@var{v}, @var{p}, @var{r}, @var{c}) +## @deftypefnx {Function File} {} discrete_rnd (@var{v}, @var{p}, @var{sz}) ## Generate a row vector containing a random sample of size @var{n} from ## the univariate distribution which assumes the values in @var{v} with -## probabilities @var{p}. +## probabilities @var{p}. @var{n} must be a scalar. ## -## Currently, @var{n} must be a scalar. +## If @var{r} and @var{c} are given create a matrix with @var{r} rows and +## @var{c} columns. Or if @var{sz} is a vector, create a matrix of size +## @var{sz}. ## @end deftypefn ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> ## Description: Random deviates from a discrete distribution -function rnd = discrete_rnd (n, v, p) +function rnd = discrete_rnd (v, p, r, c) + + if (nargin == 4) + if (! (isscalar (r) && (r > 0) && (r == round (r)))) + error ("discrete_rnd: r must be a positive integer"); + endif + if (! (isscalar (c) && (c > 0) && (c == round (c)))) + error ("discrete_rnd: c must be a positive integer"); + endif + sz = [r, c]; + elseif (nargin == 3) + ## A potential problem happens here if all args are scalar, as + ## we can't distiguish between the command syntax. Thankfully this + ## case doesn't make much sense. So we assume the first syntax + ## if the first arg is scalar - if (nargin != 3) - usage ("discrete_rnd (n, v, p)"); - endif - - if (! isscalar (n)) - error ("discrete_rnd: n must be a scalar"); + if (isscalar (v)) + sz = [1, floor(v)]; + v = p; + p = r; + else + if (isscalar (r) && (r > 0)) + sz = [r, r]; + elseif (isvector(r) && all (r > 0)) + sz = r(:)'; + else + error ("discrete_rnd: r must be a postive integer or vector"); + endif + endif + else + usage ("discrete_rnd (n, v, p) | discrete_rnd (v, p, r, c)"); endif if (! isvector (v)) @@ -47,10 +74,24 @@ error ("discrete_rnd: p must be a nonzero, nonnegative vector"); endif + n = prod (sz); + m = length (v); u = rand (1, n); - m = length (p); s = reshape (cumsum (p / sum (p)), m, 1); + ## The following loop is a space/time tradeoff in favor of space, + ## since the dataset may be large. + ## + ## Vectorized code is: + ## rnd = v (1 + sum ((s * ones (1, n)) <= ((ones (m, 1) * u)))); + rnd = reshape (rnd, sz); + ## + ## Non-vectorized code is: + ## + ## rnd = zeros (sz); + ## for q=1:n + ## rnd (q) = v (sum (s <= u (q)) + 1); + ## endfor endfunction
--- a/scripts/statistics/distributions/empirical_rnd.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/empirical_rnd.m Thu Apr 08 23:52:45 2004 +0000 @@ -19,19 +19,35 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} empirical_rnd (@var{n}, @var{data}) +## @deftypefnx {Function File} {} empirical_rnd (@var{data}, @var{r}, @var{c}) +## @deftypefnx {Function File} {} empirical_rnd (@var{data}, @var{sz}) ## Generate a bootstrap sample of size @var{n} from the empirical ## distribution obtained from the univariate sample @var{data}. +## +## If @var{r} and @var{c} are given create a matrix with @var{r} rows and +## @var{c} columns. Or if @var{sz} is a vector, create a matrix of size +## @var{sz}. ## @end deftypefn ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> ## Description: Bootstrap samples from the empirical distribution -function rnd = empirical_rnd (n, data) +function rnd = empirical_rnd (data, r, c) + + if (nargin == 2) + if (isscalar(data)) + c = data; + data = r; + r = 1; + endif + elseif (nargin != 3) + usage ("empirical_rnd (n, data) | empirical_rnd (data, r, c)"); + endif if (! isvector (data)) error ("empirical_rnd: data must be a vector"); endif - rnd = discrete_rnd (n, data, ones (size (data)) / length (data)); + rnd = discrete_rnd (data, ones (size (data)) / length (data), r, c); endfunction
--- a/scripts/statistics/distributions/exponential_cdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/exponential_cdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -35,32 +35,40 @@ usage ("exponential_cdf (x, lambda)"); endif - [retval, x, l] = common_size (x, l); - if (retval > 0) - error ("exponential_cdf: x and lambda must be of common size or scalar"); + if (!isscalar (x) && !isscalar(l)) + [retval, x, l] = common_size (x, l); + if (retval > 0) + error ("exponential_cdf: x and lambda must be of common size or scalar"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - l = reshape (l, 1, s); - cdf = zeros (1, s); + if (isscalar (x)) + sz = size (l); + else + sz = size (x); + endif + + cdf = zeros (sz); k = find (isnan (x) | !(l > 0)); if (any (k)) - cdf(k) = NaN * ones (1, length (k)); + cdf(k) = NaN; endif k = find ((x == Inf) & (l > 0)); if (any (k)) - cdf(k) = ones (1, length (k)); + cdf(k) = 1; endif k = find ((x > 0) & (x < Inf) & (l > 0)); if (any (k)) - cdf (k) = 1 - exp (- l(k) .* x(k)); + if isscalar (l) + cdf (k) = 1 - exp (- l .* x(k)); + elseif isscalar (x) + cdf (k) = 1 - exp (- l(k) .* x); + else + cdf (k) = 1 - exp (- l(k) .* x(k)); + endif endif - cdf = reshape (cdf, r, c); - endfunction
--- a/scripts/statistics/distributions/exponential_inv.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/exponential_inv.m Thu Apr 08 23:52:45 2004 +0000 @@ -33,32 +33,40 @@ usage ("exponential_inv (x, lambda)"); endif - [retval, x, l] = common_size (x, l); - if (retval > 0) - error ("exponential_inv: x and lambda must be of common size or scalar"); + if (!isscalar (x) && !isscalar(l)) + [retval, x, l] = common_size (x, l); + if (retval > 0) + error ("exponential_inv: x and lambda must be of common size or scalar"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - l = reshape (l, 1, s); - inv = zeros (1, s); + if (isscalar (x)) + sz = size (l); + else + sz = size (x); + endif + + inv = zeros (sz); k = find (!(l > 0) | (x < 0) | (x > 1) | isnan (x)); if (any (k)) - inv(k) = NaN * ones (1, length (k)); + inv(k) = NaN; endif k = find ((x == 1) & (l > 0)); if (any (k)) - inv(k) = Inf * ones (1, length (k)); + inv(k) = Inf; endif k = find ((x > 0) & (x < 1) & (l > 0)); if (any (k)) - inv(k) = - log (1 - x(k)) ./ l(k); + if isscalar (l) + inv(k) = - log (1 - x(k)) ./ l; + elseif isscalar (x) + inv(k) = - log (1 - x) ./ l(k); + else + inv(k) = - log (1 - x(k)) ./ l(k); + endif endif - inv = reshape (inv, r, c); - endfunction
--- a/scripts/statistics/distributions/exponential_pdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/exponential_pdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -32,27 +32,34 @@ usage ("exponential_pdf (x, lambda)"); endif - [retval, x, l] = common_size (x, l); - if (retval > 0) - error ("exponential_pdf: x and lambda must be of common size or scalar"); + if (!isscalar (x) && !isscalar(l)) + [retval, x, l] = common_size (x, l); + if (retval > 0) + error ("exponential_pdf: x and lambda must be of common size or scalar"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - l = reshape (l, 1, s); - pdf = zeros (1, s); + if (isscalar (x)) + sz = size (l); + else + sz = size (x); + endif + pdf = zeros (sz); k = find (!(l > 0) | isnan (x)); if (any (k)) - pdf(k) = NaN * ones (1, length (k)); + pdf(k) = NaN; endif k = find ((x > 0) & (x < Inf) & (l > 0)); if (any (k)) - pdf(k) = l(k) .* exp (- l(k) .* x(k)); + if isscalar (l) + pdf(k) = l .* exp (- l .* x(k)); + elseif isscalar (x) + pdf(k) = l(k) .* exp (- l(k) .* x); + else + pdf(k) = l(k) .* exp (- l(k) .* x(k)); + endif endif - pdf = reshape (pdf, r, c); - endfunction
--- a/scripts/statistics/distributions/exponential_rnd.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/exponential_rnd.m Thu Apr 08 23:52:45 2004 +0000 @@ -19,9 +19,11 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} exponential_rnd (@var{lambda}, @var{r}, @var{c}) +## @deftypefn {Function File} {} exponential_rnd (@var{lambda}, @var{sz}) ## Return an @var{r} by @var{c} matrix of random samples from the ## exponential distribution with parameter @var{lambda}, which must be a -## scalar or of size @var{r} by @var{c}. +## scalar or of size @var{r} by @var{c}. Or if @var{sz} is a vector, +## create a matrix of size @var{sz}. ## ## If @var{r} and @var{c} are omitted, the size of the result matrix is ## the size of @var{lambda}. @@ -39,29 +41,48 @@ if (! (isscalar (c) && (c > 0) && (c == round (c)))) error ("exponential_rnd: c must be a positive integer"); endif - [retval, l] = common_size (l, zeros (r, c)); - if (retval > 0) - error ("exponential_rnd: lambda must be scalar or of size %d by %d", - r, c); + sz = [r, c]; + + if (any (size (l) != 1) && + ((length (size (nl)) != length (sz)) || any (size (l) != sz))) + error ("exponential_rnd: lambda must be scalar or of size [r, c]"); endif - elseif (nargin != 1) + elseif (nargin == 2) + if (isscalar (r) && (r > 0)) + sz = [r, r]; + elseif (isvector(r) && all (r > 0)) + sz = r(:)'; + else + error ("exponential_rnd: r must be a postive integer or vector"); + endif + + if (any (size (l) != 1) && + ((length (size (l)) != length (sz)) || any (size (l) != sz))) + error ("exponential_rnd: lambda must be scalar or of size sz"); + endif + elseif (nargin == 1) + sz = size (l); + else usage ("exponential_rnd (lambda, r, c)"); endif - [r, c] = size (l); - s = r * c; - l = reshape (l, 1, s); - rnd = zeros (1, s); - k = find (!(l > 0) | !(l < Inf)); - if (any (k)) - rnd(k) = NaN * ones (1, length (k)); - endif - k = find ((l > 0) & (l < Inf)); - if (any (k)) - rnd(k) = - log (1 - rand (1, length (k))) ./ l(k); + if (isscalar (l)) + if ((l > 0) && (l < Inf)) + rnd = - log (1 - rand (sz)) ./ l; + else + rnd = NaN * ones (sz); + endif + else + rnd = zeros (sz); + k = find (!(l > 0) | !(l < Inf)); + if (any (k)) + rnd(k) = NaN; + endif + k = find ((l > 0) & (l < Inf)); + if (any (k)) + rnd(k) = - log (1 - rand (size (k))) ./ l(k); + endif endif - rnd = reshape (rnd, r, c); - endfunction
--- a/scripts/statistics/distributions/f_cdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/f_cdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -33,33 +33,34 @@ usage ("f_cdf (x, m, n)"); endif - [retval, x, m, n] = common_size (x, m, n); - if (retval > 0) - error ("f_cdf: x, m and n must be of common size or scalar"); + if (!isscalar (m) || !isscalar (n)) + [retval, x, m, n] = common_size (x, m, n); + if (retval > 0) + error ("f_cdf: x, m and n must be of common size or scalar"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - m = reshape (m, 1, s); - n = reshape (n, 1, s); - cdf = zeros (1, s); + sz = size (x); + cdf = zeros (sz); k = find (!(m > 0) | !(n > 0) | isnan (x)); if (any (k)) - cdf(k) = NaN * ones (1, length (k)); + cdf(k) = NaN; endif k = find ((x == Inf) & (m > 0) & (n > 0)); if (any (k)) - cdf(k) = ones (1, length (k)); + cdf(k) = 1; endif k = find ((x > 0) & (x < Inf) & (m > 0) & (n > 0)); if (any (k)) - cdf(k) = 1 - betainc (1 ./ (1 + m(k) .* x(k) ./ n(k)), n(k) / 2, m(k) / 2); + if (isscalar (m) && isscalar (n)) + cdf(k) = 1 - betainc (1 ./ (1 + m .* x(k) ./ n), n / 2, m / 2); + else + cdf(k) = 1 - betainc (1 ./ (1 + m(k) .* x(k) ./ n(k)), n(k) / 2, + m(k) / 2); + endif endif - cdf = reshape (cdf, r, c); - endfunction
--- a/scripts/statistics/distributions/f_inv.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/f_inv.m Thu Apr 08 23:52:45 2004 +0000 @@ -33,34 +33,34 @@ usage ("f_inv (x, m, n)"); endif - [retval, x, m, n] = common_size (x, m, n); - if (retval > 0) - error ("f_inv: x, m and n must be of common size or scalar"); + if (!isscalar (m) || !isscalar (n)) + [retval, x, m, n] = common_size (x, m, n); + if (retval > 0) + error ("f_inv: x, m and n must be of common size or scalar"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - m = reshape (m, 1, s); - n = reshape (n, 1, s); - inv = zeros (1, s); + sz = size (x); + inv = zeros (sz); k = find ((x < 0) | (x > 1) | isnan (x) | !(m > 0) | !(n > 0)); if (any (k)) - inv(k) = NaN * ones (1, length (k)); + inv(k) = NaN; endif k = find ((x == 1) & (m > 0) & (n > 0)); if (any (k)) - inv(k) = Inf * ones (1, length (k)); + inv(k) = Inf; endif k = find ((x > 0) & (x < 1) & (m > 0) & (n > 0)); if (any (k)) - inv(k) = ((1 ./ beta_inv (1 - x(k), n(k) / 2, m(k) / 2) - 1) - .* n(k) ./ m(k)); + if (isscalar (m) && isscalar (n)) + inv(k) = ((1 ./ beta_inv (1 - x(k), n / 2, m / 2) - 1) .* n ./ m); + else + inv(k) = ((1 ./ beta_inv (1 - x(k), n(k) / 2, m(k) / 2) - 1) + .* n(k) ./ m(k)); + endif endif - inv = reshape (inv, r, c); - -endfunction \ No newline at end of file +endfunction
--- a/scripts/statistics/distributions/f_pdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/f_pdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -33,31 +33,34 @@ usage ("f_pdf (x, m, n)"); endif - [retval, x, m, n] = common_size (x, m, n); - if (retval > 0) - error ("f_pdf: x, m and n must be of common size or scalar"); + if (!isscalar (m) || !isscalar (n)) + [retval, x, m, n] = common_size (x, m, n); + if (retval > 0) + error ("f_pdf: x, m and n must be of common size or scalar"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - m = reshape (m, 1, s); - n = reshape (n, 1, s); - pdf = zeros (1, s); + sz = size (x); + pdf = zeros (sz); k = find (isnan (x) | !(m > 0) | !(n > 0)); if (any (k)) - pdf(k) = NaN * ones (1, length (k)); + pdf(k) = NaN; endif k = find ((x > 0) & (x < Inf) & (m > 0) & (n > 0)); if (any (k)) - tmp = m(k) .* x(k) ./ n(k); - pdf(k) = (exp ((m(k) / 2 - 1) .* log (tmp) - - ((m(k) + n(k)) / 2) .* log (1 + tmp)) - .* (m(k) ./ n(k)) ./ beta (m(k) / 2, n(k) / 2)); + if (isscalar (m) && isscalar (n)) + tmp = m / n * x(k); + pdf(k) = (exp ((m / 2 - 1) .* log (tmp) + - ((m + n) / 2) .* log (1 + tmp)) + .* (m / n) ./ beta (m / 2, n / 2)); + else + tmp = m(k) .* x(k) ./ n(k); + pdf(k) = (exp ((m(k) / 2 - 1) .* log (tmp) + - ((m(k) + n(k)) / 2) .* log (1 + tmp)) + .* (m(k) ./ n(k)) ./ beta (m(k) / 2, n(k) / 2)); + endif endif - pdf = reshape (pdf, r, c); - endfunction
--- a/scripts/statistics/distributions/f_rnd.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/f_rnd.m Thu Apr 08 23:52:45 2004 +0000 @@ -19,9 +19,12 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} f_rnd (@var{m}, @var{n}, @var{r}, @var{c}) +## @deftypefnx {Function File} {} f_rnd (@var{m}, @var{n}, @var{sz}) ## Return an @var{r} by @var{c} matrix of random samples from the F ## distribution with @var{m} and @var{n} degrees of freedom. Both ## @var{m} and @var{n} must be scalar or of size @var{r} by @var{c}. +## If @var{sz} is a vector the random samples are in a matrix of +## size @var{sz}. ## ## If @var{r} and @var{c} are omitted, the size of the result matrix is ## the common size of @var{m} and @var{n}. @@ -32,6 +35,16 @@ function rnd = f_rnd (m, n, r, c) + if (nargin > 1) + if (!isscalar(m) || !isscalar(n)) + [retval, m, n] = common_size (m, n); + if (retval > 0) + error ("f_rnd: m and n must be of common size or scalar"); + endif + endif + endif + + if (nargin == 4) if (! (isscalar (r) && (r > 0) && (r == round (r)))) error ("f_rnd: r must be a positive integer"); @@ -39,37 +52,52 @@ if (! (isscalar (c) && (c > 0) && (c == round (c)))) error ("f_rnd: c must be a positive integer"); endif - [retval, m, n] = common_size (m, n, zeros (r, c)); - if (retval > 0) - error ("f_rnd: m and n must be scalar or of size %d by %d", r, c); + sz = [r, c]; + + if (any (size (m) != 1) && + ((length (size (m)) != length (sz)) || any (size (m) != sz))) + error ("f_rnd: m and n must be scalar or of size [r,c]"); + endif + elseif (nargin == 3) + if (isscalar (r) && (r > 0)) + sz = [r, r]; + elseif (isvector(r) && all (r > 0)) + sz = r(:)'; + else + error ("f_rnd: r must be a postive integer or vector"); + endif + + if (any (size (m) != 1) && + ((length (size (m)) != length (sz)) || any (size (m) != sz))) + error ("f_rnd: m and n must be scalar or of size sz"); endif elseif (nargin == 2) - [retval, m, n] = common_size (m, n); - if (retval > 0) - error ("f_rnd: m and n must be of common size or scalar"); - endif + sz = size(a); else usage ("f_rnd (m, n, r, c)"); endif - [r, c] = size (m); - s = r * c; - m = reshape (m, 1, s); - n = reshape (n, 1, s); - rnd = zeros (1, s); + + if (isscalar (m) && isscalar (n)) + if ((m > 0) && (m < Inf) && (n > 0) && (n < Inf)) + rnd = f_inv (rand (sz), m, n); + else + rnd = NaN * ones (sz); + endif + else + rnd = zeros (sz); - k = find (!(m > 0) | !(m < Inf) | - !(n > 0) | !(n < Inf)); - if (any (k)) - rnd(k) = NaN * ones (1, length (k)); + k = find (!(m > 0) | !(m < Inf) | + !(n > 0) | !(n < Inf)); + if (any (k)) + rnd(k) = NaN; + endif + + k = find ((m > 0) & (m < Inf) & + (n > 0) & (n < Inf)); + if (any (k)) + rnd(k) = f_inv (rand (size (k)), m(k), n(k)); + endif endif - k = find ((m > 0) & (m < Inf) & - (n > 0) & (n < Inf)); - if (any (k)) - rnd(k) = f_inv (rand (1, length (k)), m(k), n(k)); - endif - - rnd = reshape (rnd, r, c); - -endfunction \ No newline at end of file +endfunction
--- a/scripts/statistics/distributions/geometric_cdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/geometric_cdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -32,32 +32,34 @@ usage ("geometric_cdf (x, p)"); endif - [retval, x, p] = common_size (x, p); - if (retval > 0) - error ("geometric_cdf: x and p must be of common size or scalar"); + if (!isscalar (x) && !isscalar (p)) + [retval, x, p] = common_size (x, p); + if (retval > 0) + error ("geometric_cdf: x and p must be of common size or scalar"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - p = reshape (p, 1, s); - cdf = zeros (1, s); + cdf = zeros (size (x)); k = find (isnan (x) | !(p >= 0) | !(p <= 1)); if (any (k)) - cdf(k) = NaN * ones (1, length (k)); + cdf(k) = NaN; endif k = find ((x == Inf) & (p >= 0) & (p <= 1)); if (any (k)) - cdf(k) = ones (1, length (k)); + cdf(k) = 1; endif k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (p > 0) & (p <= 1)); if (any (k)) - cdf(k) = 1 - ((1 - p(k)) .^ (x(k) + 1)); + if (isscalar (x)) + cdf(k) = 1 - ((1 - p(k)) .^ (x + 1)); + elseif (isscalar (p)) + cdf(k) = 1 - ((1 - p) .^ (x(k) + 1)); + else + cdf(k) = 1 - ((1 - p(k)) .^ (x(k) + 1)); + endif endif - cdf = reshape (cdf, r, c); - endfunction
--- a/scripts/statistics/distributions/geometric_inv.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/geometric_inv.m Thu Apr 08 23:52:45 2004 +0000 @@ -32,33 +32,34 @@ usage ("geometric_inv (x, p)"); endif - [retval, x, p] = common_size (x, p); - if (retval > 0) - error ("geometric_inv: x and p must be of common size or scalar"); + if (!isscalar (x) && !isscalar (p)) + [retval, x, p] = common_size (x, p); + if (retval > 0) + error ("geometric_inv: x and p must be of common size or scalar"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - p = reshape (p, 1, s); - inv = zeros (1, s); + inv = zeros (size (x)); k = find (!(x >= 0) | !(x <= 1) | !(p >= 0) | !(p <= 1)); if (any (k)) - inv(k) = NaN * ones (1, length (k)); + inv(k) = NaN; endif k = find ((x == 1) & (p >= 0) & (p <= 1)); if (any (k)) - inv(k) = Inf * ones (1, length (k)); + inv(k) = Inf; endif k = find ((x > 0) & (x < 1) & (p > 0) & (p <= 1)); if (any (k)) - inv(k) = max (ceil (log (1 - x(k)) ./ log (1 - p(k))) - 1, - zeros (1, length (k))); + if (isscalar (x)) + inv(k) = max (ceil (log (1 - x) ./ log (1 - p(k))) - 1, 0); + elseif (isscalar (p)) + inv(k) = max (ceil (log (1 - x(k)) / log (1 - p)) - 1, 0); + else + inv(k) = max (ceil (log (1 - x(k)) ./ log (1 - p(k))) - 1, 0); + endif endif - inv = reshape (inv, r, c); - endfunction
--- a/scripts/statistics/distributions/geometric_pdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/geometric_pdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -32,33 +32,35 @@ usage ("geometric_pdf (x, p)"); endif - [retval, x, p] = common_size (x, p); - if (retval > 0) - error ("geometric_pdf: x and p must be of common size or scalar"); + if (!isscalar (x) && !isscalar (p)) + [retval, x, p] = common_size (x, p); + if (retval > 0) + error ("geometric_pdf: x and p must be of common size or scalar"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - p = reshape (p, 1, s); - cdf = zeros (1, s); + pdf = zeros (size (x)); k = find (isnan (x) | !(p >= 0) | !(p <= 1)); if (any (k)) - pdf(k) = NaN * ones (1, length (k)); + pdf(k) = NaN; endif ## Just for the fun of it ... k = find ((x == Inf) & (p == 0)); if (any (k)) - pdf(k) = ones (1, length (k)); + pdf(k) = 1; endif k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (p > 0) & (p <= 1)); if (any (k)) - pdf(k) = p(k) .* ((1 - p(k)) .^ x(k)); + if (isscalar (x)) + pdf(k) = p(k) .* ((1 - p(k)) .^ x); + elseif (isscalar (p)) + pdf(k) = p .* ((1 - p) .^ x(k)); + else + pdf(k) = p(k) .* ((1 - p(k)) .^ x(k)); + endif endif - pdf = reshape (pdf, r, c); - endfunction
--- a/scripts/statistics/distributions/geometric_rnd.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/geometric_rnd.m Thu Apr 08 23:52:45 2004 +0000 @@ -19,12 +19,14 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} geometric_rnd (@var{p}, @var{r}, @var{c}) +## @deftypefnx {Function File} {} geometric_rnd (@var{p}, @var{sz}) ## Return an @var{r} by @var{c} matrix of random samples from the ## geometric distribution with parameter @var{p}, which must be a scalar ## or of size @var{r} by @var{c}. ## -## If @var{r} and @var{c} are omitted, the size of the result matrix is -## the size of @var{p}. +## If @var{r} and @var{c} are given create a matrix with @var{r} rows and +## @var{c} columns. Or if @var{sz} is a vector, create a matrix of size +## @var{sz}. ## @end deftypefn ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> @@ -39,34 +41,59 @@ if (! (isscalar (c) && (c > 0) && (c == round (c)))) error ("geometric_rnd: c must be a positive integer"); endif - [retval, p] = common_size (p, zeros (r, c)); - if (retval > 0) - error ("geometric_rnd: p must be scalar or of size %d by %d", r, c); + sz = [r, c]; + + if (any (size (p) != 1) && ((length (size (p)) != length (sz)) || + any (size (p) != sz))) + error ("geometric_rnd: p must be scalar or of size [r, c]"); endif + elseif (nargin == 2) + if (isscalar (r) && (r > 0)) + sz = [r, r]; + elseif (isvector(r) && all (r > 0)) + sz = r(:)'; + else + error ("geometric_rnd: r must be a postive integer or vector"); + endif + + if (any (size (p) != 1) && ((length (size (p)) != length (sz)) || + any (size (p) != sz))) + error ("geometric_rnd: n must be scalar or of size sz"); + endif + elseif (nargin == 1) + sz = size(n); elseif (nargin != 1) usage ("geometric_rnd (p, r, c)"); endif - [r, c] = size (p); - s = r * c; - p = reshape (p, 1, s); - rnd = zeros (1, s); + + if (isscalar (p)) + if (!(p >= 0) || !(p <= 1)) + rnd = NaN * ones (sz); + elseif (p == 0) + rnd = Inf * ones (sz); + elseif ((p > 0) & (p < 1)); + rnd = floor (log (rand (sz)) / log (1 - p)); + else + rnd = zeros (sz); + endif + else + rnd = zeros (sz); - k = find (!(p >= 0) | !(p <= 1)); - if (any (k)) - rnd(k) = NaN * ones (1, length (k)); + k = find (!(p >= 0) | !(p <= 1)); + if (any (k)) + rnd(k) = NaN * ones (1, length (k)); + endif + + k = find (p == 0); + if (any (k)) + rnd(k) = Inf * ones (1, length (k)); + endif + + k = find ((p > 0) & (p < 1)); + if (any (k)) + rnd(k) = floor (log (rand (size (k))) ./ log (1 - p(k))); + endif endif - k = find (p == 0); - if (any (k)) - rnd(k) = Inf * ones (1, length (k)); - endif - - k = find ((p > 0) & (p < 1)); - if (any (k)) - rnd(k) = floor (log (rand (1, length (k))) ./ log (1 - p(k))); - endif - - rnd = reshape (rnd, r, c); - endfunction
--- a/scripts/statistics/distributions/hypergeometric_cdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/hypergeometric_cdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -36,7 +36,11 @@ function cdf = hypergeometric_cdf (x, m, t, n) if (nargin != 4) - usage ("hypergeometrix_cdf (x, m, t, n)"); + usage ("hypergeometric_cdf (x, m, t, n)"); + endif + + if (!isscalar (m) || !isscalar (t) || !isscalar (n)) + error ("hypergeometric_cdf: m, t and n must all be positive integers"); endif if ((m < 0) | (t < 0) | (n <= 0) | (m != round (m)) |
--- a/scripts/statistics/distributions/hypergeometric_inv.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/hypergeometric_inv.m Thu Apr 08 23:52:45 2004 +0000 @@ -36,6 +36,10 @@ usage ("hypergeometric_inv (x, m, t, n)"); endif + if (!isscalar (m) || !isscalar (t) || !isscalar (n)) + error ("hypergeometrix_inv: m, t and n must all be positive integers"); + endif + if ((m < 0) | (t < 0) | (n <= 0) | (m != round (m)) | (t != round (t)) | (n != round (n)) | (m > t) | (n > t)) inv = NaN * ones (size (x))
--- a/scripts/statistics/distributions/hypergeometric_pdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/hypergeometric_pdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -37,18 +37,15 @@ usage ("hypergeometric_pdf (x, m, t, n)"); endif - [retval, x, m, t, n] = common_size (x, m, t, n); - if (retval > 0) - error ("hypergeometric_pdf: x, m, t, and n must be of common size or scalar"); + if (!isscalar (m) || !isscalar (t) || !isscalar (n)) + [retval, x, m, t, n] = common_size (x, m, t, n); + if (retval > 0) + error ("hypergeometric_pdf: x, m, t, and n must be of common size or scalar"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - m = reshape (m, 1, s); - t = reshape (t, 1, s); - n = reshape (n, 1, s); - pdf = zeros * ones (1, s); + pdf = zeros (size (x)); + ## everything in i1 gives NaN i1 = ((m < 0) | (t < 0) | (n <= 0) | (m != round (m)) | (t != round (t)) | (n != round (n)) | (m > t) | (n > t)); @@ -56,12 +53,21 @@ i2 = ((x != round (x)) | (x < 0) | (x > m) | (n < x) | (n-x > t-m)); k = find (i1); if (any (k)) - pdf (k) = NaN * ones (size (k)); + if (isscalar (m) && isscalar (t) && isscalar (n)) + pdf = NaN * ones ( size (x)); + else + pdf (k) = NaN; + endif endif k = find (!i1 & !i2); if (any (k)) - pdf (k) = (bincoeff (m(k), x(k)) .* bincoeff (t(k)-m(k), n(k)-x(k)) - ./ bincoeff (t(k), n(k))); + if (isscalar (m) && isscalar (t) && isscalar (n)) + pdf (k) = (bincoeff (m, x(k)) .* bincoeff (t-m, n-x(k)) + / bincoeff (t, n)); + else + pdf (k) = (bincoeff (m(k), x(k)) .* bincoeff (t(k)-m(k), n(k)-x(k)) + ./ bincoeff (t(k), n(k))); + endif endif endfunction
--- a/scripts/statistics/distributions/hypergeometric_rnd.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/hypergeometric_rnd.m Thu Apr 08 23:52:45 2004 +0000 @@ -19,25 +19,61 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} hypergeometric_rnd (@var{n_size}, @var{m}, @var{t}, @var{n}) -## Generate a row vector containing a random sample of size @var{n_size} from -## the hypergeometric distribution with parameters @var{m}, @var{t}, and -## @var{n}. +## @deftypefnx {Function File} {} hypergeometric_rnd (@var{@var{m}, @var{t}, @var{n}, @var{r}, @var{c}) +## @deftypefnx {Function File} {} hypergeometric_rnd (@var{@var{m}, @var{t}, @var{n}, @var{sz}) +## Generate a row vector containing a random sample of size @var{n_size} +## from the hypergeometric distribution with parameters @var{m}, @var{t}, +## and @var{n}. +## +## If @var{r} and @var{c} are given create a matrix with @var{r} rows and +## @var{c} columns. Or if @var{sz} is a vector, create a matrix of size +## @var{sz}. ## ## The parameters @var{m}, @var{t}, and @var{n} must positive integers ## with @var{m} and @var{n} not greater than @var{t}. ## @end deftypefn -function rnd = hypergeometric_rnd (N, m, t, n) +## function rnd = hypergeometric_rnd (N, m, t, n) +function rnd = hypergeometric_rnd (m, t, n, r, c) - if (nargin != 4) - usage ("hypergeometric_rnd (N, m, t, n)"); + if (nargin == 5) + if (! (isscalar (r) && (r > 0) && (r == round (r)))) + error ("hypergeometric_rnd: r must be a positive integer"); + endif + if (! (isscalar (c) && (c > 0) && (c == round (c)))) + error ("hypergeometric_rnd: c must be a positive integer"); + endif + sz = [r, c]; + elseif (nargin == 4) + ## A potential problem happens here if all args are scalar, as + ## we can distiguish between the command syntax. This is quite + ## ambigous! I assume that if the last arg is a vector then + ## then third form is assumed. This means that you can't define + ## and r-by-r matrix with a single scalar! + + if (isscalar (r)) + sz = [1, floor(m)]; + m = t; + t = n; + n = r; + elseif (isvector(r) && all (r > 0)) + sz = r(:)'; + else + error ("hypergeometric_rnd: r must be a vector of positive integers"); + endif + else + usage ("hypergeometric_rnd (N, m, t, n) | hypergeometric_rnd (m, t, n, r, c)"); + endif + + if (!isscalar (m) || !isscalar (t) || !isscalar (n)) + error ("hypergeometric_cdf: m, t and n must all be positive integers"); endif if ((m < 0) | (t < 0) | (n <= 0) | (m != round (m)) | (t != round (t)) | (n != round (n)) | (m > t) | (n > t)) - rnd = NaN * ones (1, N) + rnd = NaN * ones (sz) else - rnd = discrete_rnd (N, 0 : n, hypergeometric_pdf (0 : n, m, t, n)); + rnd = discrete_rnd (0 : n, hypergeometric_pdf (0 : n, m, t, n), sz); endif endfunction
--- a/scripts/statistics/distributions/kolmogorov_smirnov_cdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/kolmogorov_smirnov_cdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -57,17 +57,20 @@ endif endif - [nr, nc] = size (x); - if (min (nr, nc) == 0) + n = numel (x); + if (n == 0) error ("kolmogorov_smirnov_cdf: x must not be empty"); endif - n = nr * nc; - x = reshape (x, 1, n); - cdf = zeros (1, n); + cdf = zeros (size (x)); + ind = find (x > 0); if (length (ind) > 0) - y = x(ind); + if (size(ind,2) < size(ind,1)) + y = x(ind.'); + else + y = x(ind); + endif K = ceil (sqrt (- log (tol) / 2) / min (y)); k = (1:K)'; A = exp (- 2 * k.^2 * y.^2); @@ -76,6 +79,4 @@ cdf(ind) = 1 + 2 * sum (A); endif - cdf = reshape (cdf, nr, nc); - endfunction
--- a/scripts/statistics/distributions/laplace_cdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/laplace_cdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -32,19 +32,16 @@ usage ("laplace_cdf (x)"); endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - cdf = zeros (1, s); + cdf = zeros (size (x)); k = find (isnan (x)); if (any (k)) - cdf(k) = NaN * ones (1, length (k)); + cdf(k) = NaN; endif k = find (x == Inf); if (any (k)) - cdf(k) = ones (1, length (k)); + cdf(k) = 1; endif k = find ((x > -Inf) & (x < Inf)); @@ -52,6 +49,4 @@ cdf(k) = (1 + sign (x(k)) .* (1 - exp (- abs (x(k))))) / 2; endif - cdf = reshape (cdf, r, c); - -endfunction \ No newline at end of file +endfunction
--- a/scripts/statistics/distributions/laplace_inv.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/laplace_inv.m Thu Apr 08 23:52:45 2004 +0000 @@ -32,19 +32,16 @@ usage ("laplace_inv (x)"); endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - inv = (-Inf) * ones (1, s); + inv = (-Inf) * ones (size (x)); k = find (isnan (x) | (x < 0) | (x > 1)); if (any (k)) - inv(k) = NaN * ones (1, length (k)); + inv(k) = NaN; endif k = find (x == 1); if (any (k)) - inv(k) = Inf * ones (1, length (k)); + inv(k) = Inf; endif k = find ((x > 0) & (x < 1)); @@ -53,6 +50,4 @@ - (x(k) > 1/2) .* log (2 * (1 - x(k)))); endif - inv = reshape (inv, r, c); - -endfunction \ No newline at end of file +endfunction
--- a/scripts/statistics/distributions/laplace_pdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/laplace_pdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -32,14 +32,11 @@ usage ("laplace_pdf (x)"); endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - pdf = zeros (1, s); + pdf = zeros (size (x)); k = find (isnan (x)); if (any (k)) - pdf(k) = NaN * ones (1, length (k)); + pdf(k) = NaN; endif k = find ((x > -Inf) & (x < Inf)); @@ -47,6 +44,4 @@ pdf(k) = exp (- abs (x(k))) / 2; endif - pdf = reshape (pdf, r, c); - -endfunction \ No newline at end of file +endfunction
--- a/scripts/statistics/distributions/laplace_rnd.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/laplace_rnd.m Thu Apr 08 23:52:45 2004 +0000 @@ -19,8 +19,10 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} laplace_rnd (@var{r}, @var{c}) +## @deftypefnx {Function File} {} laplace_rnd (@var{sz}); ## Return an @var{r} by @var{c} matrix of random numbers from the -## Laplace distribution. +## Laplace distribution. Or is @var{sz} is a vector, create a matrix of +## @var{sz}. ## @end deftypefn ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> @@ -28,18 +30,27 @@ function rnd = laplace_rnd (r, c) - if (nargin != 2) + if (nargin == 2) + if (! (isscalar (r) && (r > 0) && (r == round (r)))) + error ("laplace_rnd: r must be a positive integer"); + endif + if (! (isscalar (c) && (c > 0) && (c == round (c)))) + error ("laplace_rnd: c must be a positive integer"); + endif + sz = [r, c]; + elseif (nargin == 1) + if (isscalar (r) && (r > 0)) + sz = [r, r]; + elseif (isvector(r) && all (r > 0)) + sz = r(:)'; + else + error ("laplace_rnd: r must be a postive integer or vector"); + endif + else usage ("laplace_rnd (r, c)"); endif - if (! (isscalar (r) && (r > 0) && (r == round (r)))) - error ("laplace_rnd: r must be a positive integer"); - endif - if (! (isscalar (c) && (c > 0) && (c == round (c)))) - error ("laplace_rnd: c must be a positive integer"); - endif - - tmp = rand (r, c); + tmp = rand (sz); rnd = ((tmp < 1/2) .* log (2 * tmp) - (tmp > 1/2) .* log (2 * (1 - tmp)));
--- a/scripts/statistics/distributions/logistic_inv.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/logistic_inv.m Thu Apr 08 23:52:45 2004 +0000 @@ -32,24 +32,21 @@ usage ("logistic_inv (x)"); endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - inv = zeros (1, s); + inv = zeros (size (x)); k = find ((x < 0) | (x > 1) | isnan (x)); if (any (k)) - inv(k) = NaN * ones (1, length (k)); + inv(k) = NaN; endif k = find (x == 0); if (any (k)) - inv(k) = (-Inf) * ones (1, length (k)); + inv(k) = -Inf; endif k = find (x == 1); if (any (k)) - inv(k) = Inf * ones (1, length (k)); + inv(k) = Inf; endif k = find ((x > 0) & (x < 1)); @@ -57,6 +54,4 @@ inv (k) = - log (1 ./ x(k) - 1); endif - inv = reshape (inv, r, c); - endfunction
--- a/scripts/statistics/distributions/logistic_rnd.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/logistic_rnd.m Thu Apr 08 23:52:45 2004 +0000 @@ -19,8 +19,10 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} logistic_rnd (@var{r}, @var{c}) +## @deftypefn {Function File} {} logistic_rnd (@var{sz}) ## Return an @var{r} by @var{c} matrix of random numbers from the -## logistic distribution. +## logistic distribution. Or is @var{sz} is a vector, create a matrix of +## @var{sz}. ## @end deftypefn ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> @@ -28,17 +30,27 @@ function rnd = logistic_rnd (r, c) - if (nargin != 2) + + if (nargin == 2) + if (! (isscalar (r) && (r > 0) && (r == round (r)))) + error ("logistic_rnd: r must be a positive integer"); + endif + if (! (isscalar (c) && (c > 0) && (c == round (c)))) + error ("logistic_rnd: c must be a positive integer"); + endif + sz = [r, c]; + elseif (nargin == 1) + if (isscalar (r) && (r > 0)) + sz = [r, r]; + elseif (isvector(r) && all (r > 0)) + sz = r(:)'; + else + error ("logistic_rnd: r must be a postive integer or vector"); + endif + else usage ("logistic_rnd (r, c)"); endif - if (! (isscalar (r) && (r > 0) && (r == round (r)))) - error ("logistic_rnd: r must be a positive integer"); - endif - if (! (isscalar (c) && (c > 0) && (c == round (c)))) - error ("logistic_rnd: c must be a positive integer"); - endif - - rnd = - log (1 ./ rand (r, c) - 1); + rnd = - log (1 ./ rand (sz) - 1); endfunction
--- a/scripts/statistics/distributions/lognormal_cdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/lognormal_cdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -47,33 +47,32 @@ ## cdf = normal_cdf (log (x), log (a), v); ## Hence ... - [retval, x, a, v] = common_size (x, a, v); - if (retval > 0) - error ("lognormal_cdf: x, a and v must be of common size or scalars"); + if (!isscalar (a) || !isscalar (v)) + [retval, x, a, v] = common_size (x, a, v); + if (retval > 0) + error ("lognormal_cdf: x, a and v must be of common size or scalars"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - a = reshape (a, 1, s); - v = reshape (v, 1, s); - cdf = zeros (1, s); + cdf = zeros (size (x)); k = find (isnan (x) | !(a > 0) | !(a < Inf) | !(v > 0) | !(v < Inf)); if (any (k)) - cdf(k) = NaN * ones (1, length (k)); + cdf(k) = NaN; endif k = find ((x == Inf) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf)); if (any (k)) - cdf(k) = ones (1, length (k)); + cdf(k) = 1; endif k = find ((x > 0) & (x < Inf) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf)); if (any (k)) - cdf(k) = stdnormal_cdf ((log (x(k)) - log (a(k))) ./ sqrt (v(k))); + if (isscalar (a) && isscalar (v)) + cdf(k) = stdnormal_cdf ((log (x(k)) - log (a)) / sqrt (v)); + else + cdf(k) = stdnormal_cdf ((log (x(k)) - log (a(k))) ./ sqrt (v(k))); + endif endif - cdf = reshape (cdf, r, c); - endfunction
--- a/scripts/statistics/distributions/lognormal_inv.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/lognormal_inv.m Thu Apr 08 23:52:45 2004 +0000 @@ -47,34 +47,33 @@ ## inv = exp (normal_inv (x, log (a), v)); ## Hence ... - [retval, x, a, v] = common_size (x, a, v); - if (retval > 0) - error ("lognormal_inv: x, a and v must be of common size or scalars"); + if (!isscalar (a) || !isscalar (v)) + [retval, x, a, v] = common_size (x, a, v); + if (retval > 0) + error ("lognormal_inv: x, a and v must be of common size or scalars"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - a = reshape (a, 1, s); - v = reshape (v, 1, s); - inv = zeros (1, s); + inv = zeros (size (x)); k = find (!(x >= 0) | !(x <= 1) | !(a > 0) | !(a < Inf) | !(v > 0) | !(v < Inf)); if (any (k)) - inv(k) = NaN * ones (1, length (k)); + inv(k) = NaN; endif k = find ((x == 1) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf)); if (any (k)) - inv(k) = Inf * ones (1, length (k)); + inv(k) = Inf; endif k = find ((x > 0) & (x < 1) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf)); if (any (k)) - inv(k) = a(k) .* exp (sqrt (v(k)) .* stdnormal_inv (x(k))); + if (isscalar (a) && isscalar (v)) + inv(k) = a .* exp (sqrt (v) .* stdnormal_inv (x(k))); + else + inv(k) = a(k) .* exp (sqrt (v(k)) .* stdnormal_inv (x(k))); + endif endif - inv = reshape (inv, r, c); - endfunction
--- a/scripts/statistics/distributions/lognormal_pdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/lognormal_pdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -47,28 +47,27 @@ ## pdf = (x > 0) ./ x .* normal_pdf (log (x), log (a), v); ## Hence ... - [retval, x, a, v] = common_size (x, a, v); - if (retval > 0) - error ("lognormal_pdf: x, a and v must be of common size or scalars"); + if (!isscalar (a) || !isscalar (v)) + [retval, x, a, v] = common_size (x, a, v); + if (retval > 0) + error ("lognormal_pdf: x, a and v must be of common size or scalars"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - a = reshape (a, 1, s); - v = reshape (v, 1, s); - pdf = zeros (1, s); + pdf = zeros (size (x)); k = find (isnan (x) | !(a > 0) | !(a < Inf) | !(v > 0) | !(v < Inf)); if (any (k)) - pdf(k) = NaN * ones (1, length (k)); + pdf(k) = NaN; endif k = find ((x > 0) & (x < Inf) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf)); if (any (k)) - pdf(k) = normal_pdf (log (x(k)), log (a(k)), v(k)) ./ x(k); + if (isscalar (a) && isscalar (v)) + pdf(k) = normal_pdf (log (x(k)), log (a), v) ./ x(k); + else + pdf(k) = normal_pdf (log (x(k)), log (a(k)), v(k)) ./ x(k); + endif endif - pdf = reshape (pdf, r, c); - endfunction
--- a/scripts/statistics/distributions/lognormal_rnd.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/lognormal_rnd.m Thu Apr 08 23:52:45 2004 +0000 @@ -19,9 +19,11 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} lognormal_rnd (@var{a}, @var{v}, @var{r}, @var{c}) +## @deftypefnx {Function File} {} lognormal_rnd (@var{a}, @var{v}, @var{sz}) ## Return an @var{r} by @var{c} matrix of random samples from the ## lognormal distribution with parameters @var{a} and @var{v}. Both ## @var{a} and @var{v} must be scalar or of size @var{r} by @var{c}. +## Or if @var{sz} is a vector, create a matrix of size @var{sz}. ## ## If @var{r} and @var{c} are omitted, the size of the result matrix is ## the common size of @var{a} and @var{v}. @@ -32,6 +34,15 @@ function rnd = lognormal_rnd (a, v, r, c) + if (nargin > 1) + if (!isscalar(a) || !isscalar(v)) + [retval, a, v] = common_size (a, v); + if (retval > 0) + error ("lognormal_rnd: a and v must be of common size or scalar"); + endif + endif + endif + if (nargin == 4) if (! (isscalar (r) && (r > 0) && (r == round (r)))) error ("lognormal_rnd: r must be a positive integer"); @@ -39,35 +50,51 @@ if (! (isscalar (c) && (c > 0) && (c == round (c)))) error ("lognormal_rnd: c must be a positive integer"); endif - [retval, a, v] = common_size (a, v, zeros (r, c)); - if (retval > 0) - error ("lognormal_rnd: a and v must be scalar or of size %d by %d", r, c); + sz = [r, c]; + + if (any (size (a) != 1) && + ((length (size (a)) != length (sz)) || any (size (a) != sz))) + error ("lognormal_rnd: a and b must be scalar or of size [r, c]"); + endif + + elseif (nargin == 3) + if (isscalar (r) && (r > 0)) + sz = [r, r]; + elseif (isvector(r) && all (r > 0)) + sz = r(:)'; + else + error ("lognormal_rnd: r must be a postive integer or vector"); + endif + + if (any (size (a) != 1) && + ((length (size (a)) != length (sz)) || any (size (a) != sz))) + error ("lognormal_rnd: a and b must be scalar or of size sz"); endif elseif (nargin == 2) - [retval, a, v] = common_size (a, v); - if (retval > 0) - error ("lognormal_rnd: a and v must be of common size or scalar"); - endif + sz = size(a); else usage ("lognormal_rnd (a, v, r, c)"); endif - [r, c] = size (a); - s = r * c; - a = reshape (a, 1, s); - v = reshape (v, 1, s); - rnd = zeros (1, s); + if (isscalar (a) && isscalar (v)) + if (!(a > 0) | !(a < Inf) | !(v > 0) | !(v < Inf)) + rnd = NaN * ones (sz); + elseif find ((a > 0) & (a < Inf) & (v > 0) & (v < Inf)); + rnd = a * exp (sqrt (v) .* randn (sz)); + else + rnd = zeros (sz); + endif + else + rnd = zeros (sz); + k = find (!(a > 0) | !(a < Inf) | !(v > 0) | !(v < Inf)); + if (any (k)) + rnd(k) = NaN * ones (1, length (k)); + endif - k = find (!(a > 0) | !(a < Inf) | !(v > 0) | !(v < Inf)); - if (any (k)) - rnd(k) = NaN * ones (1, length (k)); + k = find ((a > 0) & (a < Inf) & (v > 0) & (v < Inf)); + if (any (k)) + rnd(k) = a(k) .* exp (sqrt (v(k)) .* randn (1, length (k))); + endif endif - k = find ((a > 0) & (a < Inf) & (v > 0) & (v < Inf)); - if (any (k)) - rnd(k) = a(k) .* exp (sqrt (v(k)) .* randn (1, length (k))); - endif - - rnd = reshape (rnd, r, c); - -endfunction \ No newline at end of file +endfunction
--- a/scripts/statistics/distributions/pascal_cdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/pascal_cdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -36,51 +36,58 @@ usage ("pascal_cdf (x, n, p)"); endif - [retval, x, n, p] = common_size (x, n, p); - if (retval > 0) - error ("pascal_cdf: x, n and p must be of common size or scalar"); + if (!isscalar(n) || !isscalar(p)) + [retval, x, n, p] = common_size (x, n, p); + if (retval > 0) + error ("pascal_cdf: x, n and p must be of common size or scalar"); + endif endif - - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - n = reshape (n, 1, s); - p = reshape (p, 1, s); - cdf = zeros (1, s); + + cdf = zeros (size (x)); k = find (isnan (x) | (n < 1) | (n == Inf) | (n != round (n)) | (p < 0) | (p > 1)); if (any (k)) - cdf(k) = NaN * ones (1, length (k)); + cdf(k) = NaN; endif k = find ((x == Inf) & (n > 0) & (n < Inf) & (n == round (n)) & (p >= 0) & (p <= 1)); if (any (k)) - cdf(k) = ones (1, length (k)); + cdf(k) = 1; endif k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (n > 0) & (n < Inf) & (n == round (n)) & (p > 0) & (p <= 1)); if (any (k)) ## Does anyone know a better way to do the summation? - m = zeros (1, length (k)); + m = zeros (size (k)); x = floor (x(k)); - n = n(k); - p = p(k); y = cdf(k); - while (1) - l = find (m <= x); - if (any (l)) - y(l) = y(l) + pascal_pdf (m(l), n(l), p(l)); - m(l) = m(l) + 1; - else - break; - endif - endwhile + if (isscalar (n) && isscalar (p)) + while (1) + l = find (m <= x); + if (any (l)) + y(l) = y(l) + pascal_pdf (m(l), n, p); + m(l) = m(l) + 1; + else + break; + endif + endwhile + else + n = n(k); + p = p(k); + while (1) + l = find (m <= x); + if (any (l)) + y(l) = y(l) + pascal_pdf (m(l), n(l), p(l)); + m(l) = m(l) + 1; + else + break; + endif + endwhile + endif cdf(k) = y; endif - cdf = reshape (cdf, r, c); - endfunction
--- a/scripts/statistics/distributions/pascal_inv.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/pascal_inv.m Thu Apr 08 23:52:45 2004 +0000 @@ -37,50 +37,58 @@ usage ("pascal_inv (x, n, p)"); endif - [retval, x, n, p] = common_size (x, n, p); - if (retval > 0) - error ("pascal_inv: x, n and p must be of common size or scalar"); + if (!isscalar(n) || !isscalar(p)) + [retval, x, n, p] = common_size (x, n, p); + if (retval > 0) + error ("pascal_inv: x, n and p must be of common size or scalar"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - n = reshape (n, 1, s); - p = reshape (p, 1, s); - inv = zeros (1, s); + inv = zeros (size (x)); k = find (isnan (x) | (x < 0) | (x > 1) | (n < 1) | (n == Inf) | (n != round (n)) | (p < 0) | (p > 1)); if (any (k)) - inv(k) = NaN * ones (1, length (k)); + inv(k) = NaN; endif k = find ((x == 1) & (n > 0) & (n < Inf) & (n == round (n)) & (p >= 0) & (p <= 1)); if (any (k)) - inv(k) = Inf * ones (1, length (k)); + inv(k) = Inf; endif k = find ((x >= 0) & (x < 1) & (n > 0) & (n < Inf) & (n == round (n)) & (p > 0) & (p <= 1)); if (any (k)) + m = zeros (size (k)); x = x(k); - n = n(k); - p = p(k); - m = zeros (1, length (k)); - s = p .^ n; - while (1) - l = find (s < x); - if (any (l)) - m(l) = m(l) + 1; - s(l) = s(l) + pascal_pdf (m(l), n(l), p(l)); - else - break; - endif - endwhile + if (isscalar (n) && isscalar (p)) + s = p ^ n * ones (size(k)); + while (1) + l = find (s < x); + if (any (l)) + m(l) = m(l) + 1; + s(l) = s(l) + pascal_pdf (m(l), n, p); + else + break; + endif + endwhile + else + n = n(k); + p = p(k); + s = p .^ n; + while (1) + l = find (s < x); + if (any (l)) + m(l) = m(l) + 1; + s(l) = s(l) + pascal_pdf (m(l), n(l), p(l)); + else + break; + endif + endwhile + endif inv(k) = m; endif - inv = reshape (inv, r, c); - endfunction
--- a/scripts/statistics/distributions/pascal_pdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/pascal_pdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -37,37 +37,36 @@ usage ("pascal_pdf (x, n, p)"); endif - [retval, x, n, p] = common_size (x, n, p); - if (retval > 0) - error ("pascal_pdf: x, n and p must be of common size or scalar"); + if (!isscalar(n) || !isscalar(p)) + [retval, x, n, p] = common_size (x, n, p); + if (retval > 0) + error ("pascal_pdf: x, n and p must be of common size or scalar"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - n = reshape (n, 1, s); - p = reshape (p, 1, s); - cdf = zeros (1, s); + pdf = zeros (size (x)); k = find (isnan (x) | (n < 1) | (n == Inf) | (n != round (n)) | (p < 0) | (p > 1)); if (any (k)) - pdf(k) = NaN * ones (1, length (k)); + pdf(k) = NaN; endif ## Just for the fun of it ... k = find ((x == Inf) & (n > 0) & (n < Inf) & (n == round (n)) & (p == 0)); if (any (k)) - pdf(k) = ones (1, length (k)); + pdf(k) = 1; endif k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (n > 0) & (n < Inf) & (n == round (n)) & (p > 0) & (p <= 1)); if (any (k)) - pdf(k) = bincoeff (-n(k), x(k)) .* (p(k) .^ n(k)) .* ((p(k) - 1) .^ x(k)); + if (isscalar (n) && isscalar (p)) + pdf(k) = bincoeff (-n, x(k)) .* (p ^ n) .* ((p - 1) .^ x(k)); + else + pdf(k) = bincoeff (-n(k), x(k)) .* (p(k) .^ n(k)) .* ((p(k) - 1) .^ x(k)); + endif endif - pdf = reshape (pdf, r, c); - endfunction
--- a/scripts/statistics/distributions/pascal_rnd.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/pascal_rnd.m Thu Apr 08 23:52:45 2004 +0000 @@ -19,12 +19,14 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} pascal_rnd (@var{n}, @var{p}, @var{r}, @var{c}) +## @deftypefnx {Function File} {} pascal_rnd (@var{n}, @var{p}, @var{sz}) ## Return an @var{r} by @var{c} matrix of random samples from the Pascal ## (negative binomial) distribution with parameters @var{n} and @var{p}. -## Both @var{n} and @var{p} must be scalar or of size @var{r} by @var{c}. +## Both @var{n} and @var{p} must be scalar or of size @var{r} by @var{c}. ## ## If @var{r} and @var{c} are omitted, the size of the result matrix is -## the common size of @var{n} and @var{p}. +## the common size of @var{n} and @var{p}. Or if @var{sz} is a vector, +## create a matrix of size @var{sz}. ## @end deftypefn ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> @@ -32,6 +34,15 @@ function rnd = pascal_rnd (n, p, r, c) + if (nargin > 1) + if (!isscalar(n) || !isscalar(p)) + [retval, n, p] = common_size (n, p); + if (retval > 0) + error ("pascal_rnd: n and p must be of common size or scalar"); + endif + endif + endif + if (nargin == 4) if (! (isscalar (r) && (r > 0) && (r == round (r)))) error ("pascal_rnd: r must be a positive integer"); @@ -39,39 +50,71 @@ if (! (isscalar (c) && (c > 0) && (c == round (c)))) error ("pascal_rnd: c must be a positive integer"); endif - [retval, n, p] = common_size (n, p, zeros (r, c)); - if (retval > 0) - error ("pascal_rnd: n and p must be scalar or of size %d by %d", r, c); + sz = [r, c]; + + if (any (size (n) != 1) && + ((length (size (n)) != length (sz)) || any (size (n) != sz))) + error ("pascal_rnd: n and p must be scalar or of size [r, c]"); + endif + + elseif (nargin == 3) + if (isscalar (r) && (r > 0)) + sz = [r, r]; + elseif (isvector(r) && all (r > 0)) + sz = r(:)'; + else + error ("pascal_rnd: r must be a postive integer or vector"); + endif + + if (any (size (n) != 1) && + ((length (size (n)) != length (sz)) || any (size (n) != sz))) + error ("pascal_rnd: n and p must be scalar or of size sz"); endif elseif (nargin == 2) - [retval, n, p] = common_size (n, p); - if (retval > 0) - error ("pascal_rnd: n and p must be of common size or scalar"); - endif + sz = size(n); else usage ("pascal_rnd (n, p, r, c)"); endif - [r, c] = size (n); - s = r * c; - n = reshape (n, 1, s); - p = reshape (p, 1, s); - rnd = zeros (1, s); + if (isscalar (n) && isscalar (p)) + if ((n < 1) || (n == Inf) || (n != round (n)) || (p <= 0) || (p > 1)); + rnd = NaN * ones (sz) + elseif ((n > 0) && (n < Inf) && (n == round (n)) && + (p > 0) && (p <= 1)) + L = prod (sz); + tmp = floor (log (rand (n, L)) / log (1 - p)); + if (n == 1) + rnd = tmp; + else + ind = (1 : n)' * ones (1, L); + rnd = sum (tmp .* (ind <= n)); + endif + else + rnd = zeros (sz); + endif + else + rnd = zeros (sz); - k = find (!(n > 0) | !(n < Inf) | !(n == round (n)) | !(p <= 0) | !(p >= 1)); - if (any (k)) - rnd(k) = NaN * ones (1, length (k)); + k = find ((n < 1) || (n == Inf) || (n != round (n)) || + (p <= 0) || (p > 1)); + if (any (k)) + rnd(k) = NaN; + endif + + k = find ((n > 0) & (n < Inf) & (n == round (n)) & (p > 0) & (p <= 1)); + if (any (k)) + n = reshape (n, 1, prod (sz)); + p = reshape (p, 1, prod (sz)); + N = max (n(k)); + L = length (k); + tmp = floor (log (rand (N, L)) ./ (ones (N, 1) * log (1 - p(k)))); + if (N == 1) + rnd(k) = tmp; + else + ind = (1 : N)' * ones (1, L); + rnd(k) = sum (tmp .* (ind <= ones (N, 1) * n(k))); + endif + endif endif - k = find ((n > 0) & (n < Inf) & (n == round (n)) & (p >= 0) & (p <= 1)); - if (any (k)) - N = max (n(k)); - L = length (k); - tmp = floor (log (rand (N, L)) ./ (ones (N, 1) * log (1 - p(k)))); - ind = (1 : N)' * ones (1, L); - rnd(k) = sum (tmp .* (ind <= ones (N, 1) * n(k))); - endif - - rnd = reshape (rnd, r, c); - -endfunction \ No newline at end of file +endfunction
--- a/scripts/statistics/distributions/poisson_cdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/poisson_cdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -33,32 +33,32 @@ usage ("poisson_cdf (x, lambda)"); endif - [retval, x, l] = common_size (x, l); - if (retval > 0) - error ("poisson_cdf: x and lambda must be of common size or scalar"); + if (!isscalar (l)) + [retval, x, l] = common_size (x, l); + if (retval > 0) + error ("poisson_cdf: x and lambda must be of common size or scalar"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - l = reshape (l, 1, s); - cdf = zeros (1, s); + cdf = zeros (size (x)); k = find (isnan (x) | !(l > 0)); if (any (k)) - cdf(k) = NaN * ones (1, length (k)); + cdf(k) = NaN; endif k = find ((x == Inf) & (l > 0)); if (any (k)) - cdf(k) = ones (1, length (k)); + cdf(k) = 1; endif k = find ((x >= 0) & (x < Inf) & (l > 0)); if (any (k)) - cdf(k) = 1 - gammainc (l(k), floor (x(k)) + 1); + if (isscalar (l)) + cdf(k) = 1 - gammainc (l, floor (x(k)) + 1); + else + cdf(k) = 1 - gammainc (l(k), floor (x(k)) + 1); + endif endif - cdf = reshape (cdf, r, c); - endfunction
--- a/scripts/statistics/distributions/poisson_inv.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/poisson_inv.m Thu Apr 08 23:52:45 2004 +0000 @@ -33,41 +33,45 @@ usage ("poisson_inv (x, lambda)"); endif - [retval, x, l] = common_size (x, l); - if (retval > 0) - error ("poisson_inv: x and lambda must be of common size or scalar"); + if (!isscalar (l)) + [retval, x, l] = common_size (x, l); + if (retval > 0) + error ("poisson_inv: x and lambda must be of common size or scalar"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - l = reshape (l, 1, s); - inv = zeros (1, s); + inv = zeros (size (x)); k = find ((x < 0) | (x > 1) | isnan (x) | !(l > 0)); if (any (k)) - inv(k) = NaN * ones (1, length (k)); + inv(k) = NaN; endif k = find ((x == 1) & (l > 0)); if (any (k)) - inv(k) = Inf * ones (1, length (k)); + inv(k) = Inf; endif k = find ((x > 0) & (x < 1) & (l > 0)); if (any (k)) - cdf = exp (-l(k)); + if (isscalar (l)) + cdf = exp (-l) * ones (size (k)); + else + cdf = exp (-l(k)); + endif while (1) m = find (cdf < x(k)); if (any (m)) inv(k(m)) = inv(k(m)) + 1; - cdf(m) = cdf(m) + poisson_pdf (inv(k(m)), l(k(m))); + if (isscalar (l)) + cdf(m) = cdf(m) + poisson_pdf (inv(k(m)), l); + else + cdf(m) = cdf(m) + poisson_pdf (inv(k(m)), l(k(m))); + endif else break; endif endwhile endif - inv = reshape (inv, r, c); - -endfunction \ No newline at end of file +endfunction
--- a/scripts/statistics/distributions/poisson_pdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/poisson_pdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -32,27 +32,27 @@ usage ("poisson_pdf (x, lambda)"); endif - [retval, x, l] = common_size (x, l); - if (retval > 0) - error ("poisson_pdf: x and lambda must be of common size or scalar"); + if (!isscalar (l)) + [retval, x, l] = common_size (x, l); + if (retval > 0) + error ("poisson_pdf: x and lambda must be of common size or scalar"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - l = reshape (l, 1, s); - pdf = zeros (1, s); + pdf = zeros (size (x)); k = find (!(l > 0) | isnan (x)); if (any (k)) - pdf(k) = NaN * ones (1, length (k)); + pdf(k) = NaN; endif k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (l > 0)); if (any (k)) - pdf(k) = exp (x(k) .* log (l(k)) - l(k) - gammaln (x(k) + 1)); + if (isscalar (l)) + pdf(k) = exp (x(k) .* log (l) - l - gammaln (x(k) + 1)); + else + pdf(k) = exp (x(k) .* log (l(k)) - l(k) - gammaln (x(k) + 1)); + endif endif - pdf = reshape (pdf, r, c); - endfunction
--- a/scripts/statistics/distributions/poisson_rnd.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/poisson_rnd.m Thu Apr 08 23:52:45 2004 +0000 @@ -39,43 +39,76 @@ if (! (isscalar (c) && (c > 0) && (c == round (c)))) error ("poisson_rnd: c must be a positive integer"); endif - [retval, l] = common_size (l, zeros (r, c)); - if (retval > 0) - error ("poisson_rnd: lambda must be scalar or of size %d by %d", r, c); + sz = [r, c]; + + if (any (size (l) != 1) && + ((length (size (l)) != length (sz)) || any (size (l) != sz))) + error ("poisson_rnd: lambda must be scalar or of size [r, c]"); endif - elseif (nargin != 1) + elseif (nargin == 2) + if (isscalar (r) && (r > 0)) + sz = [r, r]; + elseif (isvector(r) && all (r > 0)) + sz = r(:)'; + else + error ("poisson_rnd: r must be a postive integer or vector"); + endif + + if (any (size (l) != 1) && + ((length (size (l)) != length (sz)) || any (size (l) != sz))) + error ("poisson_rnd: lambda must be scalar or of size sz"); + endif + elseif (nargin == 1) + sz = size (l); + else usage ("poisson_rnd (lambda, r, c)"); endif - [r, c] = size (l); - s = r * c; - l = reshape (l, 1, s); - rnd = zeros (1, s); + if (isscalar (l)) + + if (!(l > 0) | !(l < Inf)) + rnd = NaN * ones (sz); + elseif ((l > 0) & (l < Inf)) + num = zeros (sz); + sum = - log (1 - rand (sz)) ./ l; + while (1) + ind = find (sum < 1); + if (any (ind)) + sum(ind) = (sum(ind) - log (1 - rand (size (ind))) / l); + num(ind) = num(ind) + 1; + else + break; + endif + endwhile + rnd = num; + else + rnd = zeros (sz); + endif + else + rnd = zeros (sz); - k = find (!(l > 0) | !(l < Inf)); - if (any (k)) - rnd(k) = NaN * ones (1, length (k)); + k = find (!(l > 0) | !(l < Inf)); + if (any (k)) + rnd(k) = NaN; + endif + + k = find ((l > 0) & (l < Inf)); + if (any (k)) + l = l(k); + num = zeros (size (k)); + sum = - log (1 - rand (size (k))) ./ l; + while (1) + ind = find (sum < 1); + if (any (ind)) + sum(ind) = (sum(ind) + - log (1 - rand (1, length (ind))) ./ l(ind)); + num(ind) = num(ind) + 1; + else + break; + endif + endwhile + rnd(k) = num; + endif endif - k = find ((l > 0) & (l < Inf)); - if (any (k)) - l = l(k); - len = length (k); - num = zeros (1, len); - sum = - log (1 - rand (1, len)) ./ l; - while (1) - ind = find (sum < 1); - if (any (ind)) - sum(ind) = (sum(ind) - - log (1 - rand (1, length (ind))) ./ l(ind)); - num(ind) = num(ind) + 1; - else - break; - endif - endwhile - rnd(k) = num; - endif - - rnd = reshape (rnd, r, c); - endfunction
--- a/scripts/statistics/distributions/t_cdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/t_cdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -33,36 +33,36 @@ usage ("t_cdf (x, n)"); endif - [retval, x, n] = common_size (x, n); - if (retval > 0) - error ("t_cdf: x and n must be of common size or scalar"); + if (!isscalar (n)) + [retval, x, n] = common_size (x, n); + if (retval > 0) + error ("t_cdf: x and n must be of common size or scalar"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - n = reshape (n, 1, s); - cdf = zeros (1, s); + cdf = zeros (size (x)); k = find (isnan (x) | !(n > 0)); if (any (k)) - cdf(k) = NaN * ones (1, length (k)); + cdf(k) = NaN; endif k = find ((x == Inf) & (n > 0)); if (any (k)) - cdf(k) = ones (1, length (k)); + cdf(k) = 1; endif k = find ((x > -Inf) & (x < Inf) & (n > 0)); if (any (k)) - cdf(k) = betainc (1 ./ (1 + x(k) .^ 2 ./ n(k)), n(k) / 2, 1 / 2) / 2; + if (isscalar (n)) + cdf(k) = betainc (1 ./ (1 + x(k) .^ 2 ./ n), n / 2, 1 / 2) / 2; + else + cdf(k) = betainc (1 ./ (1 + x(k) .^ 2 ./ n(k)), n(k) / 2, 1 / 2) / 2; + endif ind = find (x(k) > 0); if (any (ind)) cdf(k(ind)) = 1 - cdf(k(ind)); endif endif - cdf = reshape (cdf, r, c); - endfunction
--- a/scripts/statistics/distributions/t_inv.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/t_inv.m Thu Apr 08 23:52:45 2004 +0000 @@ -37,37 +37,41 @@ usage ("t_inv (x, n)"); endif - [retval, x, n] = common_size (x, n); - if (retval > 0) - error ("t_inv: x and n must be of common size or scalar"); + if (!isscalar (n)) + [retval, x, n] = common_size (x, n); + if (retval > 0) + error ("t_inv: x and n must be of common size or scalar"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - n = reshape (n, 1, s); - inv = zeros (1, s); + inv = zeros (size (x)); k = find ((x < 0) | (x > 1) | isnan (x) | !(n > 0)); if (any (k)) - inv(k) = NaN * ones (1, length (k)); + inv(k) = NaN; endif k = find ((x == 0) & (n > 0)); if (any (k)) - inv(k) = (-Inf) * ones (1, length (k)); + inv(k) = -Inf; endif k = find ((x == 1) & (n > 0)); if (any (k)) - inv(k) = Inf * ones (1, length (k)); + inv(k) = Inf; endif k = find ((x > 0) & (x < 1) & (n > 0) & (n < 10000)); if (any (k)) - inv(k) = (sign (x(k) - 1/2) - .* sqrt (n(k) .* (1 ./ beta_inv (2*min (x(k), 1 - x(k)), - n(k)/2, 1/2) - 1))); + if (isscalar (n)) + inv(k) = (sign (x(k) - 1/2) + .* sqrt (n .* (1 ./ beta_inv (2*min (x(k), 1 - x(k)), + n/2, 1/2) - 1))); + else + inv(k) = (sign (x(k) - 1/2) + .* sqrt (n(k) .* (1 ./ beta_inv (2*min (x(k), 1 - x(k)), + n(k)/2, 1/2) - 1))); + endif endif ## For large n, use the quantiles of the standard normal @@ -76,6 +80,4 @@ inv(k) = stdnormal_inv (x(k)); endif - inv = reshape (inv, r, c); - -endfunction \ No newline at end of file +endfunction
--- a/scripts/statistics/distributions/t_pdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/t_pdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -33,28 +33,29 @@ usage ("t_pdf (x, n)"); endif - [retval, x, n] = common_size (x, n); - if (retval > 0) - error ("t_pdf: x and n must be of common size or scalar"); + if (!isscalar (n)) + [retval, x, n] = common_size (x, n); + if (retval > 0) + error ("t_pdf: x and n must be of common size or scalar"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - n = reshape (n, 1, s); - pdf = zeros (1, s); + pdf = zeros (size (x)); k = find (isnan (x) | !(n > 0) | !(n < Inf)); if (any (k)) - pdf(k) = NaN * ones (1, length (k)); + pdf(k) = NaN; endif k = find (!isinf (x) & !isnan (x) & (n > 0) & (n < Inf)); if (any (k)) - pdf(k) = (exp (- (n(k) + 1) .* log (1 + x(k) .^ 2 ./ n(k))/2) - ./ (sqrt (n(k)) .* beta (n(k)/2, 1/2))); + if (isscalar (n)) + pdf(k) = (exp (- (n + 1) .* log (1 + x(k) .^ 2 ./ n)/2) + / (sqrt (n) * beta (n/2, 1/2))); + else + pdf(k) = (exp (- (n(k) + 1) .* log (1 + x(k) .^ 2 ./ n(k))/2) + ./ (sqrt (n(k)) .* beta (n(k)/2, 1/2))); + endif endif - pdf = reshape (pdf, r, c); - endfunction
--- a/scripts/statistics/distributions/t_rnd.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/t_rnd.m Thu Apr 08 23:52:45 2004 +0000 @@ -19,9 +19,11 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} t_rnd (@var{n}, @var{r}, @var{c}) +## @deftypefnx {Function File} {} t_rnd (@var{n}, @var{sz}) ## Return an @var{r} by @var{c} matrix of random samples from the t ## (Student) distribution with @var{n} degrees of freedom. @var{n} must -## be a scalar or of size @var{r} by @var{c}. +## be a scalar or of size @var{r} by @var{c}. Or if @va{sz} is a +## vector create a matrix of size @var{sz}. ## ## If @var{r} and @var{c} are omitted, the size of the result matrix is ## the size of @var{n}. @@ -39,29 +41,51 @@ if (! (isscalar (c) && (c > 0) && (c == round (c)))) error ("t_rnd: c must be a positive integer"); endif - [retval, n] = common_size (n, zeros (r, c)); - if (retval > 0) - error ("t_rnd: n must be scalar or of size %d by %d", r, c); + sz = [r, c]; + + if (any (size (n) != 1) && + ((length (size (n)) != length (sz)) || any (size (n) != sz))) + error ("t_rnd: n must be scalar or of size sz"); endif - elseif (nargin != 1) + elseif (nargin == 2) + if (isscalar (r) && (r > 0)) + sz = [r, r]; + elseif (isvector(r) && all (r > 0)) + sz = r(:)'; + else + error ("t_rnd: r must be a postive integer or vector"); + endif + + if (any (size (n) != 1) && + ((length (size (n)) != length (sz)) || any (size (n) != sz))) + error ("t_rnd: n must be scalar or of size sz"); + endif + elseif (nargin == 1) + sz = size (n); + else usage ("t_rnd (n, r, c)"); endif - [r, c] = size (n); - s = r * c; - n = reshape (n, 1, s); - rnd = zeros (1, s); + if (isscalar (n)) + if (!(n > 0) || !(n < Inf)) + rnd = NaN * ones (sz); + elseif ((n > 0) && (n < Inf)) + rnd = t_inv (rand (sz), n); + else + rnd = zeros (size (n)); + endif + else + rnd = zeros (size (n)); - k = find (!(n > 0) | !(n < Inf)); - if (any (k)) - rnd(k) = NaN * ones (1, length (k)); + k = find (!(n > 0) | !(n < Inf)); + if (any (k)) + rnd(k) = NaN; + endif + + k = find ((n > 0) & (n < Inf)); + if (any (k)) + rnd(k) = t_inv (rand (size (k)), n(k)); + endif endif - k = find ((n > 0) & (n < Inf)); - if (any (k)) - rnd(k) = t_inv (rand (1, length (k)), n(k)); - endif - - rnd = reshape (rnd, r, c); - endfunction
--- a/scripts/statistics/distributions/weibull_cdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/weibull_cdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -40,36 +40,34 @@ usage ("weibull_cdf (x, alpha, sigma)"); endif - [retval, x, shape, scale] = common_size (x, shape, scale); - if (retval > 0) - error ("weibull_cdf: x, alpha and sigma must be of common size or scalar"); + if (!isscalar (shape) || !isscalar (scale)) + [retval, x, shape, scale] = common_size (x, shape, scale); + if (retval > 0) + error ("weibull_cdf: x, alpha and sigma must be of common size or scalar"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - shape = reshape (shape, 1, s); - scale = reshape (scale, 1, s); - - cdf = NaN * ones (1, s); + cdf = NaN * ones (size (x)); ok = ((shape > 0) & (shape < Inf) & (scale > 0) & (scale < Inf)); k = find ((x <= 0) & ok); if (any (k)) - cdf(k) = zeros (1, length (k)); + cdf(k) = 0; endif k = find ((x > 0) & (x < Inf) & ok); if (any (k)) - cdf(k) = 1 - exp (- (x(k) ./ scale(k)) .^ shape(k)); + if (isscalar (shape) && isscalar (scale)) + cdf(k) = 1 - exp (- (x(k) / scale) .^ shape); + else + cdf(k) = 1 - exp (- (x(k) ./ scale(k)) .^ shape(k)); + endif endif k = find ((x == Inf) & ok); if (any (k)) - cdf(k) = ones (1, length (k)); + cdf(k) = 1; endif - cdf = reshape (cdf, r, c); - endfunction
--- a/scripts/statistics/distributions/weibull_inv.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/weibull_inv.m Thu Apr 08 23:52:45 2004 +0000 @@ -33,35 +33,34 @@ usage ("weibull_inv (x, alpha, sigma)"); endif - [retval, x, shape, scale] = common_size (x, shape, scale); - if (retval > 0) - error ("weibull_inv: x, alpha and sigma must be of common size or scalar"); + if (!isscalar (shape) || !isscalar (scale)) + [retval, x, shape, scale] = common_size (x, shape, scale); + if (retval > 0) + error ("weibull_inv: x, alpha and sigma must be of common size or scalar"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - shape = reshape (shape, 1, s); - scale = reshape (scale, 1, s); + inv = NaN * ones (size (x)); - inv = NaN * ones (1, s); ok = ((shape > 0) & (shape < Inf) & (scale > 0) & (scale < Inf)); k = find ((x == 0) & ok); if (any (k)) - inv(k) = -Inf * ones (1, length (k)); + inv(k) = -Inf; endif k = find ((x > 0) & (x < 1) & ok); if (any (k)) - inv(k) = scale(k) .* (- log (1 - x(k))) .^ (1 ./ shape(k)); + if (isscalar (shape) && isscalar (scale)) + inv(k) = scale * (- log (1 - x(k))) .^ (1 / shape); + else + inv(k) = scale(k) .* (- log (1 - x(k))) .^ (1 ./ shape(k)); + endif endif k = find ((x == 1) & ok); if (any (k)) - inv(k) = Inf * ones (1, length (k)); + inv(k) = Inf; endif - inv = reshape (inv, r, c); - endfunction
--- a/scripts/statistics/distributions/weibull_pdf.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/weibull_pdf.m Thu Apr 08 23:52:45 2004 +0000 @@ -40,32 +40,32 @@ usage ("weibull_pdf (x, alpha, sigma)"); endif - [retval, x, shape, scale] = common_size (x, shape, scale); - if (retval > 0) - error ("weibull_pdf: x, alpha and sigma must be of common size or scalar"); + if (!isscalar (shape) || !isscalar (scale)) + [retval, x, shape, scale] = common_size (x, shape, scale); + if (retval > 0) + error ("weibull_pdf: x, alpha and sigma must be of common size or scalar"); + endif endif - [r, c] = size (x); - s = r * c; - x = reshape (x, 1, s); - shape = reshape (shape, 1, s); - scale = reshape (scale, 1, s); - - pdf = NaN * ones (1, s); + pdf = NaN * ones (size (x)); ok = ((shape > 0) & (shape < Inf) & (scale > 0) & (scale < Inf)); k = find ((x > -Inf) & (x <= 0) & ok); if (any (k)) - pdf(k) = zeros (1, length (k)); + pdf(k) = 0; endif k = find ((x > 0) & (x < Inf) & ok); if (any (k)) - pdf(k) = (shape(k) .* (scale(k) .^ -shape(k)) - .* (x(k) .^ (shape(k) - 1)) - .* exp(- (x(k) ./ scale(k)) .^ shape(k))); + if (isscalar (shape) && isscalar (scale)) + pdf(k) = (shape .* (scale .^ -shape) + .* (x(k) .^ (shape - 1)) + .* exp(- (x(k) / scale) .^ shape)); + else + pdf(k) = (shape(k) .* (scale(k) .^ -shape(k)) + .* (x(k) .^ (shape(k) - 1)) + .* exp(- (x(k) ./ scale(k)) .^ shape(k))); + endif endif - pdf = reshape (pdf, r, c); - endfunction
--- a/scripts/statistics/distributions/weibull_rnd.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/weibull_rnd.m Thu Apr 08 23:52:45 2004 +0000 @@ -19,9 +19,11 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} weibull_rnd (@var{alpha}, @var{sigma}, @var{r}, @var{c}) +## @deftypefnx {Function File} {} weibull_rnd (@var{alpha}, @var{sigma}, @var{sz}) ## Return an @var{r} by @var{c} matrix of random samples from the ## Weibull distribution with parameters @var{alpha} and @var{sigma} -## which must be scalar or of size @var{r} by @var{c}. +## which must be scalar or of size @var{r} by @var{c}. Or if @var{sz} +## is a vector return a matrix of size @var{sz}. ## ## If @var{r} and @var{c} are omitted, the size of the result matrix is ## the common size of @var{alpha} and @var{sigma}. @@ -32,6 +34,15 @@ function rnd = weibull_rnd (shape, scale, r, c) + if (nargin > 1) + if (!isscalar(shape) || !isscalar(scale)) + [retval, shape, scale] = common_size (shape, scale); + if (retval > 0) + error ("weibull_rnd: shape and scale must be of common size or scalar"); + endif + endif + endif + if (nargin == 4) if (! (isscalar (r) && (r > 0) && (r == round (r)))) error ("weibull_rnd: r must be a positive integer"); @@ -39,32 +50,46 @@ if (! (isscalar (c) && (c > 0) && (c == round (c)))) error ("weibull_rnd: c must be a positive integer"); endif - [retval, shape, scale] = common_size (shape, scale, zeros (r, c)); - if (retval > 0) - error ("weibull_rnd: alpha and sigma must be scalar or of size %d by %d", - r, c); + sz = [r, c]; + + if (any (size (scale) != 1) && + ((length (size (scale)) != length (sz)) || + any (size (scale) != sz))) + error ("weilbull_rnd: shape and scale must be scalar or of size [r, c]"); + endif + elseif (nargin == 3) + if (isscalar (r) && (r > 0)) + sz = [r, r]; + elseif (isvector(r) && all (r > 0)) + sz = r(:)'; + else + error ("weibull_rnd: r must be a postive integer or vector"); + endif + + if (any (size (scale) != 1) && + ((length (size (scale)) != length (sz)) || + any (size (scale) != sz))) + error ("weibull_rnd: shape and scale must be scalar or of size sz"); endif elseif (nargin == 2) - [retval, shape, scale] = common_size (shape, scale); - if (retval > 0) - error ("weibull_rnd: alpha and sigma must be of common size or scalar"); - endif + sz = size(shape); else usage ("weibull_rnd (alpha, sigma, r, c)"); endif - [r, c] = size (shape); - s = r * c; - shape = reshape (shape, 1, s); - scale = reshape (scale, 1, s); - - rnd = NaN * ones (1, s); - k = find ((shape > 0) & (shape < Inf) & (scale > 0) & (scale < Inf)); - if (any (k)) - rnd(k) = (scale(k) - .* (- log (1 - rand (1, length (k)))) .^ (1 ./ shape(k))); + if (isscalar (shape) && isscalar (scale)) + if ((shape > 0) & (shape < Inf) & (scale > 0) & (scale < Inf)) + rnd = (scale * (- log (1 - rand (sz))) .^ (1 / shape)); + else + rnd = NaN * ones (sz); + endif + else + rnd = NaN * ones (sz); + k = find ((shape > 0) & (shape < Inf) & (scale > 0) & (scale < Inf)); + if (any (k)) + rnd(k) = (scale(k) + .* (- log (1 - rand (size (k)))) .^ (1 ./ shape(k))); + endif endif - rnd = reshape (rnd, r, c); - -endfunction \ No newline at end of file +endfunction
--- a/scripts/statistics/distributions/wiener_rnd.m Thu Apr 08 18:00:43 2004 +0000 +++ b/scripts/statistics/distributions/wiener_rnd.m Thu Apr 08 23:52:45 2004 +0000 @@ -40,7 +40,11 @@ elseif (nargin == 2) n = 1000; elseif (nargin > 3) - usage ("wiener_rnd (t, d,n)"); + usage ("wiener_rnd (t, d, n)"); + endif + + if (!isscalar (t) || !isscalar (d) || !isscalar (n)) + error ("wiener_rnd: t, d and n must all be positive integers"); endif retval = randn (n * t, d);