Mercurial > octave-nkf
changeset 5410:56e066f5efc1
[project @ 2005-07-13 17:43:35 by jwe]
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/betacdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,65 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} beta_cdf (@var{x}, @var{a}, @var{b}) +## For each element of @var{x}, returns the CDF at @var{x} of the beta +## distribution with parameters @var{a} and @var{b}, i.e., +## PROB (beta (@var{a}, @var{b}) <= @var{x}). +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the Beta distribution + +function cdf = beta_cdf (x, a, b) + + if (nargin != 3) + usage ("beta_cdf (a, b, x)"); + endif + + if (!isscalar (a) || !isscalar(b)) + [retval, x, a, b] = common_size (x, a, b); + if (retval > 0) + error ("beta_cdf: x, a and b must be of common size or scalar"); + endif + endif + + sz = size(x); + cdf = zeros (sz); + + k = find (!(a > 0) | !(b > 0) | isnan (x)); + if (any (k)) + cdf (k) = NaN; + endif + + k = find ((x >= 1) & (a > 0) & (b > 0)); + if (any (k)) + cdf (k) = 1; + endif + + k = find ((x > 0) & (x < 1) & (a > 0) & (b > 0)); + if (any (k)) + if (isscalar (a) && isscalar(b)) + cdf (k) = betainc (x(k), a, b); + else + cdf (k) = betainc (x(k), a(k), b(k)); + endif + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/betainv.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,97 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} beta_inv (@var{x}, @var{a}, @var{b}) +## For each component of @var{x}, compute the quantile (the inverse of +## the CDF) at @var{x} of the Beta distribution with parameters @var{a} +## and @var{b}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the Beta distribution + +function inv = beta_inv (x, a, b) + + if (nargin != 3) + usage ("beta_inv (x, a, b)"); + endif + + if (!isscalar (a) || !isscalar(b)) + [retval, x, a, b] = common_size (x, a, b); + if (retval > 0) + error ("beta_inv: x, a and b must be of common size or scalars"); + endif + endif + + sz = size (x); + inv = zeros (sz); + + k = find ((x < 0) | (x > 1) | !(a > 0) | !(b > 0) | isnan (x)); + if (any (k)) + inv (k) = NaN; + endif + + k = find ((x == 1) & (a > 0) & (b > 0)); + if (any (k)) + inv (k) = 1; + endif + + k = find ((x > 0) & (x < 1) & (a > 0) & (b > 0)); + if (any (k)) + if (!isscalar(a) || !isscalar(b)) + a = a (k); + b = b (k); + y = a ./ (a + b); + else + y = a / (a + b) * ones (size (k)); + endif + x = x (k); + l = find (y < eps); + if (any (l)) + y(l) = sqrt (eps) * ones (length (l), 1); + endif + l = find (y > 1 - eps); + if (any (l)) + y(l) = 1 - sqrt (eps) * ones (length (l), 1); + endif + + y_old = y; + for i = 1 : 10000 + h = (beta_cdf (y_old, a, b) - x) ./ beta_pdf (y_old, a, b); + y_new = y_old - h; + ind = find (y_new <= eps); + if (any (ind)) + y_new (ind) = y_old (ind) / 10; + endif + ind = find (y_new >= 1 - eps); + if (any (ind)) + y_new (ind) = 1 - (1 - y_old (ind)) / 10; + endif + h = y_old - y_new; + if (max (abs (h)) < sqrt (eps)) + break; + endif + y_old = y_new; + endfor + + inv (k) = y_new; + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/betapdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,61 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} beta_pdf (@var{x}, @var{a}, @var{b}) +## For each element of @var{x}, returns the PDF at @var{x} of the beta +## distribution with parameters @var{a} and @var{b}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the Beta distribution + +function pdf = beta_pdf (x, a, b) + + if (nargin != 3) + usage ("beta_pdf (a, b, x)"); + endif + + if (!isscalar (a) || !isscalar(b)) + [retval, x, a, b] = common_size (x, a, b); + if (retval > 0) + error ("beta_pdf: x, a and b must be of common size or scalar"); + endif + endif + + sz = size (x); + pdf = zeros (sz); + + k = find (!(a > 0) | !(b > 0) | isnan (x)); + if (any (k)) + pdf (k) = NaN; + endif + + k = find ((x > 0) & (x < 1) & (a > 0) & (b > 0)); + if (any (k)) + if (isscalar(a) && isscalar(b)) + pdf(k) = exp ((a - 1) .* log (x(k)) + + (b - 1) .* log (1 - x(k))) ./ beta (a, b); + else + pdf(k) = exp ((a(k) - 1) .* log (x(k)) + + (b(k) - 1) .* log (1 - x(k))) ./ beta (a(k), b(k)); + endif + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/betarnd.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,98 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} beta_rnd (@var{a}, @var{b}, @var{r}, @var{c}) +## @deftypefnx {Function File} {} beta_rnd (@var{a}, @var{b}, @var{sz}) +## Return an @var{r} by @var{c} or @code{size (@var{sz})} matrix of +## random samples from the Beta distribution with parameters @var{a} and +## @var{b}. Both @var{a} and @var{b} 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{a} and @var{b}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the Beta distribution + +function rnd = beta_rnd (a, b, r, c) + + if (nargin > 1) + if (!isscalar(a) || !isscalar(b)) + [retval, a, b] = common_size (a, b); + if (retval > 0) + error ("beta_rnd: a and b must be of common size or scalar"); + endif + endif + endif + + if (nargin == 4) + if (! (isscalar (r) && (r > 0) && (r == round (r)))) + error ("beta_rnd: r must be a positive integer"); + endif + if (! (isscalar (c) && (c > 0) && (c == round (c)))) + error ("beta_rnd: c must be a positive integer"); + endif + sz = [r, c]; + + if (any (size (a) != 1) + && (length (size (a)) != length (sz) || any (size (a) != sz))) + error ("beta_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 ("beta_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 ("beta_rnd: a and b must be scalar or of size sz"); + endif + elseif (nargin == 2) + sz = size(a); + else + usage ("beta_rnd (a, b, r, c)"); + endif + + if (isscalar(a) && isscalar(b)) + if (find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf))) + rnd = NaN * ones (sz); + else + rnd = beta_inv (rand(sz), a, b); + endif + else + rnd = zeros (sz); + + k = find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf)); + if (any (k)) + rnd(k) = NaN * ones (size (k)); + endif + + k = find ((a > 0) & (a < Inf) & (b > 0) & (b < Inf)); + if (any (k)) + rnd(k) = beta_inv (rand (size (k)), a(k), b(k)); + endif + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/binocdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,68 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} binomial_cdf (@var{x}, @var{n}, @var{p}) +## For each element of @var{x}, compute the CDF at @var{x} of the +## binomial distribution with parameters @var{n} and @var{p}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the binomial distribution + +function cdf = binomial_cdf (x, n, p) + + if (nargin != 3) + usage ("binomial_cdf (x, n, p)"); + endif + + if (!isscalar (n) || !isscalar (p)) + [retval, x, n, p] = common_size (x, n, p); + if (retval > 0) + error ("binomial_cdf: x, n and p must be of common size or scalar"); + endif + endif + + sz = size (x); + cdf = zeros (sz); + + k = find (isnan (x) | !(n >= 0) | (n != round (n)) + | !(p >= 0) | !(p <= 1)); + if (any (k)) + cdf(k) = NaN; + endif + + k = find ((x >= n) & (n >= 0) & (n == round (n)) + & (p >= 0) & (p <= 1)); + if (any (k)) + cdf(k) = 1; + endif + + k = find ((x >= 0) & (x < n) & (n == round (n)) + & (p >= 0) & (p <= 1)); + if (any (k)) + tmp = floor (x(k)); + if (isscalar (n) && isscalar (p)) + cdf(k) = 1 - betainc (p, tmp + 1, n - tmp); + else + cdf(k) = 1 - betainc (p(k), tmp + 1, n(k) - tmp); + endif + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/binoinv.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,79 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} binomial_inv (@var{x}, @var{n}, @var{p}) +## For each element of @var{x}, compute the quantile at @var{x} of the +## binomial distribution with parameters @var{n} and @var{p}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the binomial distribution + +function inv = binomial_inv (x, n, p) + + if (nargin != 3) + usage ("binomial_inv (x, n, p)"); + endif + + if (!isscalar (n) || !isscalar (p)) + [retval, x, n, p] = common_size (x, n, p); + if (retval > 0) + error ("binomial_inv: x, n and p must be of common size or scalars"); + endif + endif + + sz = size (x); + inv = zeros (sz); + + k = find (!(x >= 0) | !(x <= 1) | !(n >= 0) | (n != round (n)) + | !(p >= 0) | !(p <= 1)); + if (any (k)) + inv(k) = NaN; + endif + + k = find ((x >= 0) & (x <= 1) & (n >= 0) & (n == round (n)) + & (p >= 0) & (p <= 1)); + if (any (k)) + if (isscalar (n) && isscalar (p)) + cdf = binomial_pdf (0, n, p) * ones (size(k)); + while (any (inv(k) < n)) + m = find (cdf < x(k)); + if (any (m)) + inv(k(m)) = inv(k(m)) + 1; + cdf(m) = cdf(m) + binomial_pdf (inv(k(m)), n, p); + else + break; + endif + endwhile + else + cdf = binomial_pdf (0, n(k), p(k)); + while (any (inv(k) < n(k))) + m = find (cdf < x(k)); + if (any (m)) + inv(k(m)) = inv(k(m)) + 1; + cdf(m) = cdf(m) + binomial_pdf (inv(k(m)), n(k(m)), p(k(m))); + else + break; + endif + endwhile + endif + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/binopdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,61 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} binomial_pdf (@var{x}, @var{n}, @var{p}) +## For each element of @var{x}, compute the probability density function +## (PDF) at @var{x} of the binomial distribution with parameters @var{n} +## and @var{p}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the binomial distribution + +function pdf = binomial_pdf (x, n, p) + + if (nargin != 3) + usage ("binomial_pdf (x, n, p)"); + endif + + if (! isscalar (n) || ! isscalar (p)) + [retval, x, n, p] = common_size (x, n, p); + if (retval > 0) + error ("binomial_pdf: x, n and p must be of common size or scalar"); + endif + endif + + k = ((x >= 0) & (x <= n) + & (x == round (x)) & (n == round (n)) + & (p >= 0) & (p <= 1)); + + pdf = zeros (size (x)); + pdf(! k) = NaN; + if (any (k(:))) + x = x(k); + if (! isscalar (n)) + n = n(k); + endif + if (! isscalar (p)) + p = p(k); + endif + z = gammaln(n+1) - gammaln(x+1) - gammaln(n-x+1) + x.*log(p) + (n-x).*log(1-p); + pdf(k) = exp (z); + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/binornd.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,108 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} binomial_rnd (@var{n}, @var{p}, @var{r}, @var{c}) +## @deftypefnx {Function File} {} binomial_rnd (@var{n}, @var{p}, @var{sz}) +## Return an @var{r} by @var{c} or a @code{size (@var{sz})} matrix of +## random samples from the 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}. +## +## If @var{r} and @var{c} are omitted, the size of the result matrix is +## the common size of @var{n} and @var{p}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the binomial distribution + +function rnd = binomial_rnd (n, p, r, c) + + if (nargin > 1) + if (!isscalar(n) || !isscalar(p)) + [retval, n, p] = common_size (n, p); + if (retval > 0) + error ("binomial_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 ("binomial_rnd: r must be a positive integer"); + endif + if (! (isscalar (c) && (c > 0) && (c == round (c)))) + error ("binomial_rnd: c must be a positive integer"); + endif + sz = [r, c]; + + if (any (size (n) != 1) + && (length (size (n)) != length (sz) || any (size (n) != sz))) + error ("binomial_rnd: n and 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 ("binomial_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 ("binomial_rnd: n and must be scalar or of size sz"); + endif + elseif (nargin == 2) + sz = size(n); + else + usage ("binomial_rnd (n, p, r, c)"); + endif + + if (isscalar (n) && isscalar (p)) + if (find (!(n > 0) | !(n < Inf) | !(n == round (n)) | + !(p >= 0) | !(p <= 1))) + rnd = NaN * ones (sz); + else + nel = prod (sz); + tmp = rand (n, nel); + rnd = sum(tmp < ones (n, nel) * p, 1); + rnd = reshape(rnd, 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; + 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 = rand (N, L); + ind = (1 : N)' * ones (1, L); + rnd(k) = sum ((tmp < ones (N, 1) * p(k)(:)') & + (ind <= ones (N, 1) * n(k)(:)'),1); + endif + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/chi2cdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,45 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} chisquare_cdf (@var{x}, @var{n}) +## For each element of @var{x}, compute the cumulative distribution +## function (CDF) at @var{x} of the chisquare distribution with @var{n} +## degrees of freedom. +## @end deftypefn + +## Author: TT <Teresa.Twaroch@ci.tuwien.ac.at> +## Description: CDF of the chi-square distribution + +function cdf = chisquare_cdf (x, n) + + if (nargin != 2) + usage ("chisquare_cdf (x, n)"); + endif + + if (!isscalar (n)) + [retval, x, n] = common_size (x, n); + if (retval > 0) + error ("chisquare_cdf: x and n must be of common size or scalar"); + endif + endif + + cdf = gamma_cdf (x, n / 2, 1 / 2); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/chi2inv.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,45 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} chisquare_inv (@var{x}, @var{n}) +## For each element of @var{x}, compute the quantile (the inverse of the +## CDF) at @var{x} of the chisquare distribution with @var{n} degrees of +## freedom. +## @end deftypefn + +## Author: TT <Teresa.Twaroch@ci.tuwien.ac.at> +## Description: Quantile function of the chi-square distribution + +function inv = chisquare_inv (x, n) + + if (nargin != 2) + usage ("chisquare_inv (x, n)"); + endif + + if (!isscalar (n)) + [retval, x, n] = common_size (x, n); + if (retval > 0) + error ("chisquare_inv: x and n must be of common size or scalar"); + endif + endif + + inv = gamma_inv (x, n / 2, 1 / 2); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/chi2pdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,45 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} chisquare_pdf (@var{x}, @var{n}) +## For each element of @var{x}, compute the probability density function +## (PDF) at @var{x} of the chisquare distribution with @var{k} degrees +## of freedom. +## @end deftypefn + +## Author: TT <Teresa.Twaroch@ci.tuwien.ac.at> +## Description: PDF of the chi-sqaure distribution + +function pdf = chisquare_pdf (x, n) + + if (nargin != 2) + usage ("chisquare_pdf (x, n)"); + endif + + if (!isscalar (n)) + [retval, x, n] = common_size (x, n); + if (retval > 0) + error ("chisquare_pdf: x and n must be of common size or scalar"); + endif + endif + + pdf = gamma_pdf (x, n / 2, 1 / 2); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/chi2rnd.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,92 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} chisquare_rnd (@var{n}, @var{r}, @var{c}) +## @deftypefnx {Function File} {} chisquare_rnd (@var{n}, @var{sz}) +## Return an @var{r} by @var{c} or a @code{size (@var{sz})} matrix of +## random samples from the chisquare distribution with @var{n} degrees +## of freedom. @var{n} 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{n}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the chi-square distribution + +function rnd = chisquare_rnd (n, r, c) + + if (nargin == 3) + if (! (isscalar (r) && (r > 0) && (r == round (r)))) + error ("chisquare_rnd: r must be a positive integer"); + endif + if (! (isscalar (c) && (c > 0) && (c == round (c)))) + error ("chisquare_rnd: c must be a positive integer"); + endif + sz = [r, c]; + + if (any (size (n) != 1) + && (length (size (n)) != length (sz) || any (size (n) != sz))) + error ("chisquare_rnd: n 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 ("chisquare_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 ("chisquare_rnd: n must be scalar or of size sz"); + endif + elseif (nargin == 1) + sz = size(n); + else + usage ("chisquare_rnd (n, r, c)"); + endif + + if (isscalar (n)) + if (find (!(n > 0) | !(n < Inf))) + rnd = NaN * ones (sz); + else + rnd = chisquare_inv (rand (sz), n); + endif + else + [retval, n, dummy] = common_size (n, ones (sz)); + if (retval > 0) + error ("chisquare_rnd: a and b must be of common size or scalar"); + endif + + rnd = zeros (sz); + k = find (!(n > 0) | !(n < Inf)); + if (any (k)) + rnd(k) = NaN; + endif + + k = find ((n > 0) & (n < Inf)); + if (any (k)) + rnd(k) = chisquare_inv (rand (size (k)), n(k)); + endif + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/expcdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,74 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} exponential_cdf (@var{x}, @var{lambda}) +## For each element of @var{x}, compute the cumulative distribution +## function (CDF) at @var{x} of the exponential distribution with +## parameter @var{lambda}. +## +## The arguments can be of common size or scalar. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the exponential distribution + +function cdf = exponential_cdf (x, l) + + if (nargin != 2) + usage ("exponential_cdf (x, lambda)"); + endif + + 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 + + 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; + endif + + k = find ((x == Inf) & (l > 0)); + if (any (k)) + cdf(k) = 1; + endif + + k = find ((x > 0) & (x < Inf) & (l > 0)); + if (any (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 + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/expinv.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,72 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} exponential_inv (@var{x}, @var{lambda}) +## For each element of @var{x}, compute the quantile (the inverse of the +## CDF) at @var{x} of the exponential distribution with parameter +## @var{lambda}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the exponential distribution + +function inv = exponential_inv (x, l) + + if (nargin != 2) + usage ("exponential_inv (x, lambda)"); + endif + + 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 + + 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; + endif + + k = find ((x == 1) & (l > 0)); + if (any (k)) + inv(k) = Inf; + endif + + k = find ((x > 0) & (x < 1) & (l > 0)); + if (any (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 + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/exppdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,65 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} exponential_pdf (@var{x}, @var{lambda}) +## For each element of @var{x}, compute the probability density function +## (PDF) of the exponential distribution with parameter @var{lambda}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the exponential distribution + +function pdf = exponential_pdf (x, l) + + if (nargin != 2) + usage ("exponential_pdf (x, lambda)"); + endif + + 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 + + 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; + endif + + k = find ((x > 0) & (x < Inf) & (l > 0)); + if (any (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 + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/exprnd.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,88 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} exponential_rnd (@var{lambda}, @var{r}, @var{c}) +## @deftypefnx {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}. 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}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the exponential distribution + +function rnd = exponential_rnd (l, r, c) + + if (nargin == 3) + if (! (isscalar (r) && (r > 0) && (r == round (r)))) + error ("exponential_rnd: r must be a positive integer"); + endif + if (! (isscalar (c) && (c > 0) && (c == round (c)))) + error ("exponential_rnd: c must be a positive integer"); + endif + 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 == 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 + + + 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 + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/fcdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,66 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} f_cdf (@var{x}, @var{m}, @var{n}) +## For each element of @var{x}, compute the CDF at @var{x} of the F +## distribution with @var{m} and @var{n} degrees of freedom, i.e., +## PROB (F (@var{m}, @var{n}) <= @var{x}). +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the F distribution + +function cdf = f_cdf (x, m, n) + + if (nargin != 3) + usage ("f_cdf (x, m, n)"); + endif + + 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 + + sz = size (x); + cdf = zeros (sz); + + k = find (!(m > 0) | !(n > 0) | isnan (x)); + if (any (k)) + cdf(k) = NaN; + endif + + k = find ((x == Inf) & (m > 0) & (n > 0)); + if (any (k)) + cdf(k) = 1; + endif + + k = find ((x > 0) & (x < Inf) & (m > 0) & (n > 0)); + if (any (k)) + 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 + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/finv.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,66 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} f_inv (@var{x}, @var{m}, @var{n}) +## For each component of @var{x}, compute the quantile (the inverse of +## the CDF) at @var{x} of the F distribution with parameters @var{m} and +## @var{n}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the F distribution + +function inv = f_inv (x, m, n) + + if (nargin != 3) + usage ("f_inv (x, m, n)"); + endif + + 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 + + sz = size (x); + inv = zeros (sz); + + k = find ((x < 0) | (x > 1) | isnan (x) | !(m > 0) | !(n > 0)); + if (any (k)) + inv(k) = NaN; + endif + + k = find ((x == 1) & (m > 0) & (n > 0)); + if (any (k)) + inv(k) = Inf; + endif + + k = find ((x > 0) & (x < 1) & (m > 0) & (n > 0)); + if (any (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 + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/fpdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,66 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} f_pdf (@var{x}, @var{m}, @var{n}) +## For each element of @var{x}, compute the probability density function +## (PDF) at @var{x} of the F distribution with @var{m} and @var{n} +## degrees of freedom. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the F distribution + +function pdf = f_pdf (x, m, n) + + if (nargin != 3) + usage ("f_pdf (x, m, n)"); + endif + + 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 + + sz = size (x); + pdf = zeros (sz); + + k = find (isnan (x) | !(m > 0) | !(n > 0)); + if (any (k)) + pdf(k) = NaN; + endif + + k = find ((x > 0) & (x < Inf) & (m > 0) & (n > 0)); + if (any (k)) + 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 + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/frnd.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,103 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- 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}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the F distribution + +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"); + endif + if (! (isscalar (c) && (c > 0) && (c == round (c)))) + error ("f_rnd: c must be a positive integer"); + endif + 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) + sz = size(a); + else + usage ("f_rnd (m, n, r, c)"); + endif + + + 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; + 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 + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/gamcdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,60 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} gamma_cdf (@var{x}, @var{a}, @var{b}) +## For each element of @var{x}, compute the cumulative distribution +## function (CDF) at @var{x} of the Gamma distribution with parameters +## @var{a} and @var{b}. +## @end deftypefn + +## Author: TT <Teresa.Twaroch@ci.tuwien.ac.at> +## Description: CDF of the Gamma distribution + +function cdf = gamma_cdf (x, a, b) + + if (nargin != 3) + usage ("gamma_cdf (x, a, b)"); + endif + + if (!isscalar (a) || !isscalar(b)) + [retval, x, a, b] = common_size (x, a, b); + if (retval > 0) + error ("gamma_cdf: x, a and b must be of common size or scalars"); + endif + endif + + sz = size (x); + cdf = zeros (sz); + + k = find (!(a > 0) | !(b > 0) | isnan (x)); + if (any (k)) + cdf (k) = NaN; + endif + + k = find ((x > 0) & (a > 0) & (b > 0)); + if (any (k)) + if (isscalar (a) && isscalar(b)) + cdf (k) = gammainc (b * x(k), a); + else + cdf (k) = gammainc (b(k) .* x(k), a(k)); + endif + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/gaminv.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,89 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} gamma_inv (@var{x}, @var{a}, @var{b}) +## For each component of @var{x}, compute the quantile (the inverse of +## the CDF) at @var{x} of the Gamma distribution with parameters @var{a} +## and @var{b}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the Gamma distribution + +function inv = gamma_inv (x, a, b) + + if (nargin != 3) + usage ("gamma_inv (x, a, b)"); + endif + + if (!isscalar (a) || !isscalar(b)) + [retval, x, a, b] = common_size (x, a, b); + if (retval > 0) + error ("gamma_inv: x, a and b must be of common size or scalars"); + endif + endif + + sz = size (x); + inv = zeros (sz); + + k = find ((x < 0) | (x > 1) | isnan (x) | !(a > 0) | !(b > 0)); + if (any (k)) + inv (k) = NaN; + endif + + k = find ((x == 1) & (a > 0) & (b > 0)); + if (any (k)) + inv (k) = Inf; + endif + + k = find ((x > 0) & (x < 1) & (a > 0) & (b > 0)); + if (any (k)) + if (!isscalar(a) || !isscalar(b)) + a = a (k); + b = b (k); + y = a ./ b; + else + y = a / b * ones (size (k)); + endif + x = x (k); + l = find (x < eps); + if (any (l)) + y(l) = sqrt (eps) * ones (length (l), 1); + endif + + y_old = y; + for i = 1 : 100 + h = (gamma_cdf (y_old, a, b) - x) ./ gamma_pdf (y_old, a, b); + y_new = y_old - h; + ind = find (y_new <= eps); + if (any (ind)) + y_new (ind) = y_old (ind) / 10; + h = y_old - y_new; + endif + if (max (abs (h)) < sqrt (eps)) + break; + endif + y_old = y_new; + endfor + + inv (k) = y_new; + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/gampdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,73 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} gamma_pdf (@var{x}, @var{a}, @var{b}) +## For each element of @var{x}, return the probability density function +## (PDF) at @var{x} of the Gamma distribution with parameters @var{a} +## and @var{b}. +## @end deftypefn + +## Author: TT <Teresa.Twaroch@ci.tuwien.ac.at> +## Description: PDF of the Gamma distribution + +function pdf = gamma_pdf (x, a, b) + + if (nargin != 3) + usage ("gamma_pdf (x, a, b)"); + endif + + if (!isscalar (a) || !isscalar(b)) + [retval, x, a, b] = common_size (x, a, b); + if (retval > 0) + error ("gamma_pdf: x, a and b must be of common size or scalars"); + endif + endif + + sz = size(x); + pdf = zeros (sz); + + k = find (!(a > 0) | !(b > 0) | isnan (x)); + if (any (k)) + pdf (k) = NaN; + endif + + k = find ((x > 0) & (a > 0) & (a <= 1) & (b > 0)); + if (any (k)) + if (isscalar(a) && isscalar(b)) + pdf(k) = ((b .^ a) .* (x(k) .^ (a - 1)) + .* exp(-b .* x(k)) ./ gamma (a)); + else + pdf(k) = ((b(k) .^ a(k)) .* (x(k) .^ (a(k) - 1)) + .* exp(-b(k) .* x(k)) ./ gamma (a(k))); + endif + endif + + k = find ((x > 0) & (a > 1) & (b > 0)); + if (any (k)) + if (isscalar(a) && isscalar(b)) + pdf(k) = exp (a .* log (b) + (a-1) .* log (x(k)) + - b .* x(k) - gammaln (a)); + else + pdf(k) = exp (a(k) .* log (b(k)) + (a(k)-1) .* log (x(k)) + - b(k) .* x(k) - gammaln (a(k))); + endif + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/gamrnd.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,97 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} gamma_rnd (@var{a}, @var{b}, @var{r}, @var{c}) +## @deftypefnx {Function File} {} gamma_rnd (@var{a}, @var{b}, @var{sz}) +## Return an @var{r} by @var{c} or a @code{size (@var{sz})} matrix of +## random samples from the Gamma distribution with parameters @var{a} +## and @var{b}. Both @var{a} and @var{b} 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{a} and @var{b}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the Gamma distribution + +function rnd = gamma_rnd (a, b, r, c) + + if (nargin > 1) + if (!isscalar(a) || !isscalar(b)) + [retval, a, b] = common_size (a, b); + if (retval > 0) + error ("gamma_rnd: a and b must be of common size or scalar"); + endif + endif + endif + + if (nargin == 4) + if (! (isscalar (r) && (r > 0) && (r == round (r)))) + error ("gamma_rnd: r must be a positive integer"); + endif + if (! (isscalar (c) && (c > 0) && (c == round (c)))) + error ("gamma_rnd: c must be a positive integer"); + endif + sz = [r, c]; + + if (any (size (a) != 1) + && (length (size (a)) != length (sz) || any (size (a) != sz))) + error ("gamma_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 ("uniform_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 ("gamma_rnd: a and b must be scalar or of size sz"); + endif + elseif (nargin == 2) + sz = size(a); + else + usage ("gamma_rnd (a, b, r, c)"); + endif + + rnd = zeros (sz); + + if (isscalar (a) && isscalar(b)) + if (find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf))) + rnd = NaN * ones (sz); + else + rnd = gamma_inv (rand (sz), a, b); + endif + else + k = find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf)); + if (any (k)) + rnd(k) = NaN; + endif + k = find ((a > 0) & (a < Inf) & (b > 0) & (b < Inf)); + if (any (k)) + rnd(k) = gamma_inv (rand (size (k)), a(k), b(k)); + endif + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/geocdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,65 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} geometric_cdf (@var{x}, @var{p}) +## For each element of @var{x}, compute the CDF at @var{x} of the +## geometric distribution with parameter @var{p}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the geometric distribution + +function cdf = geometric_cdf (x, p) + + if (nargin != 2) + usage ("geometric_cdf (x, p)"); + endif + + 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 + + cdf = zeros (size (x)); + + k = find (isnan (x) | !(p >= 0) | !(p <= 1)); + if (any (k)) + cdf(k) = NaN; + endif + + k = find ((x == Inf) & (p >= 0) & (p <= 1)); + if (any (k)) + cdf(k) = 1; + endif + + k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (p > 0) & (p <= 1)); + if (any (k)) + 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 + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/geoinv.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,65 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} geometric_inv (@var{x}, @var{p}) +## For each element of @var{x}, compute the quantile at @var{x} of the +## geometric distribution with parameter @var{p}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the geometric distribution + +function inv = geometric_inv (x, p) + + if (nargin != 2) + usage ("geometric_inv (x, p)"); + endif + + 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 + + inv = zeros (size (x)); + + k = find (!(x >= 0) | !(x <= 1) | !(p >= 0) | !(p <= 1)); + if (any (k)) + inv(k) = NaN; + endif + + k = find ((x == 1) & (p >= 0) & (p <= 1)); + if (any (k)) + inv(k) = Inf; + endif + + k = find ((x > 0) & (x < 1) & (p > 0) & (p <= 1)); + if (any (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 + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/geopdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,66 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} geometric_pdf (@var{x}, @var{p}) +## For each element of @var{x}, compute the probability density function +## (PDF) at @var{x} of the geometric distribution with parameter @var{p}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the geometric distribution + +function pdf = geometric_pdf (x, p) + + if (nargin != 2) + usage ("geometric_pdf (x, p)"); + endif + + 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 + + pdf = zeros (size (x)); + + k = find (isnan (x) | !(p >= 0) | !(p <= 1)); + if (any (k)) + pdf(k) = NaN; + endif + + ## Just for the fun of it ... + k = find ((x == Inf) & (p == 0)); + if (any (k)) + pdf(k) = 1; + endif + + k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (p > 0) & (p <= 1)); + if (any (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 + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/geornd.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,99 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- 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 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 the geometric distribution + +function rnd = geometric_rnd (p, r, c) + + if (nargin == 3) + if (! (isscalar (r) && (r > 0) && (r == round (r)))) + error ("geometric_rnd: r must be a positive integer"); + endif + if (! (isscalar (c) && (c > 0) && (c == round (c)))) + error ("geometric_rnd: c must be a positive integer"); + endif + 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 + + + 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)); + 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 + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/hygecdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,55 @@ +## Copyright (C) 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} hypergeometric_cdf (@var{x}, @var{m}, @var{t}, @var{n}) +## Compute the cumulative distribution function (CDF) at @var{x} of the +## hypergeometric distribution with parameters @var{m}, @var{t}, and +## @var{n}. This is the probability of obtaining not more than @var{x} +## marked items when randomly drawing a sample of size @var{n} without +## replacement from a population of total size @var{t} containing +## @var{m} marked items. +## +## 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 + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the hypergeometric distribution + +function cdf = hypergeometric_cdf (x, m, t, n) + + if (nargin != 4) + 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)) | + (t != round (t)) | (n != round (n)) | (m > t) | (n > t)) + cdf = NaN * ones (size (x)) + else + cdf = discrete_cdf (x, 0 : n, hypergeometric_pdf (0 : n, m, t, n)); + endif + +endfunction + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/hygeinv.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,50 @@ +## Copyright (C) 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} hypergeometric_inv (@var{x}, @var{m}, @var{t}, @var{n}) +## For each element of @var{x}, compute the quantile at @var{x} of the +## hypergeometric distribution with parameters @var{m}, @var{t}, and +## @var{n}. +## +## 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 + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the hypergeometric distribution + +function inv = hypergeometric_inv (x, m, t, n) + + if (nargin != 4) + 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)) + else + inv = discrete_inv (x, 0 : n, hypergeometric_pdf (0 : n, m, t, n)); + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/hygepdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,73 @@ +## Copyright (C) 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} hypergeometric_pdf (@var{x}, @var{m}, @var{t}, @var{n}) +## Compute the probability density function (PDF) at @var{x} of the +## hypergeometric distribution with parameters @var{m}, @var{t}, and +## @var{n}. This is the probability of obtaining @var{x} marked items +## when randomly drawing a sample of size @var{n} without replacement +## from a population of total size @var{t} containing @var{m} marked items. +## +## The arguments must be of common size or scalar. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the hypergeometric distribution + +function pdf = hypergeometric_pdf (x, m, t, n) + + if (nargin != 4) + usage ("hypergeometric_pdf (x, m, t, n)"); + endif + + 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 + + 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)); + ## everything in i2 gives 0 unless in i1 + i2 = ((x != round (x)) | (x < 0) | (x > m) | (n < x) | (n-x > t-m)); + k = find (i1); + if (any (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)) + 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/hygernd.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,79 @@ +## Copyright (C) 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} hypergeometric_rnd (@var{n_size}, @var{m}, @var{t}, @var{n}) +## @deftypefnx {Function File} {} hypergeometric_rnd (@var{m}, @var{t}, @var{n}, @var{r}, @var{c}) +## @deftypefnx {Function File} {} hypergeometric_rnd (@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 (m, t, n, r, c) + + 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 (sz) + else + rnd = discrete_rnd (0 : n, hypergeometric_pdf (0 : n, m, t, n), sz); + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/logncdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,78 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} lognormal_cdf (@var{x}, @var{a}, @var{v}) +## For each element of @var{x}, compute the cumulative distribution +## function (CDF) at @var{x} of the lognormal distribution with +## parameters @var{a} and @var{v}. If a random variable follows this +## distribution, its logarithm is normally distributed with mean +## @code{log (@var{a})} and variance @var{v}. +## +## Default values are @var{a} = 1, @var{v} = 1. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the log normal distribution + +function cdf = lognormal_cdf (x, a, v) + + if (! ((nargin == 1) || (nargin == 3))) + usage ("lognormal_cdf (x, a, v)"); + endif + + if (nargin == 1) + a = 1; + v = 1; + endif + + ## The following "straightforward" implementation unfortunately does + ## not work (because exp (Inf) -> NaN etc): + ## cdf = normal_cdf (log (x), log (a), v); + ## Hence ... + + 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 + + cdf = zeros (size (x)); + + k = find (isnan (x) | !(a > 0) | !(a < Inf) | !(v > 0) | !(v < Inf)); + if (any (k)) + cdf(k) = NaN; + endif + + k = find ((x == Inf) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf)); + if (any (k)) + cdf(k) = 1; + endif + + k = find ((x > 0) & (x < Inf) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf)); + if (any (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 + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/logninv.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,79 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} lognormal_inv (@var{x}, @var{a}, @var{v}) +## For each element of @var{x}, compute the quantile (the inverse of the +## CDF) at @var{x} of the lognormal distribution with parameters @var{a} +## and @var{v}. If a random variable follows this distribution, its +## logarithm is normally distributed with mean @code{log (@var{a})} and +## variance @var{v}. +## +## Default values are @var{a} = 1, @var{v} = 1. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the log normal distribution + +function inv = lognormal_inv (x, a, v) + + if (! ((nargin == 1) || (nargin == 3))) + usage ("lognormal_inv (x, a, v)"); + endif + + if (nargin == 1) + a = 1; + v = 1; + endif + + ## The following "straightforward" implementation unfortunately does + ## not work (because exp (Inf) -> NaN): + ## inv = exp (normal_inv (x, log (a), v)); + ## Hence ... + + 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 + + inv = zeros (size (x)); + + k = find (!(x >= 0) | !(x <= 1) | !(a > 0) | !(a < Inf) + | !(v > 0) | !(v < Inf)); + if (any (k)) + inv(k) = NaN; + endif + + k = find ((x == 1) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf)); + if (any (k)) + inv(k) = Inf; + endif + + k = find ((x > 0) & (x < 1) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf)); + if (any (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 + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/lognpdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,73 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} lognormal_pdf (@var{x}, @var{a}, @var{v}) +## For each element of @var{x}, compute the probability density function +## (PDF) at @var{x} of the lognormal distribution with parameters +## @var{a} and @var{v}. If a random variable follows this distribution, +## its logarithm is normally distributed with mean @code{log (@var{a})} +## and variance @var{v}. +## +## Default values are @var{a} = 1, @var{v} = 1. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the log normal distribution + +function pdf = lognormal_pdf (x, a, v) + + if (! ((nargin == 1) || (nargin == 3))) + usage ("lognormal_pdf (x, a, v)"); + endif + + if (nargin == 1) + a = 1; + v = 1; + endif + + ## The following "straightforward" implementation unfortunately does + ## not work for the special cases (Inf, ...) + ## pdf = (x > 0) ./ x .* normal_pdf (log (x), log (a), v); + ## Hence ... + + 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 + + pdf = zeros (size (x)); + + k = find (isnan (x) | !(a > 0) | !(a < Inf) | !(v > 0) | !(v < Inf)); + if (any (k)) + pdf(k) = NaN; + endif + + k = find ((x > 0) & (x < Inf) & (a > 0) & (a < Inf) & (v > 0) & (v < Inf)); + if (any (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 + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/lognrnd.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,100 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- 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}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the log normal distribution + +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"); + endif + if (! (isscalar (c) && (c > 0) && (c == round (c)))) + error ("lognormal_rnd: c must be a positive integer"); + endif + 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) + sz = size(a); + else + usage ("lognormal_rnd (a, v, r, c)"); + endif + + 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) = a(k) .* exp (sqrt (v(k)) .* randn (1, length (k))); + endif + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/normcdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,73 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} normal_cdf (@var{x}, @var{m}, @var{v}) +## For each element of @var{x}, compute the cumulative distribution +## function (CDF) at @var{x} of the normal distribution with mean +## @var{m} and variance @var{v}. +## +## Default values are @var{m} = 0, @var{v} = 1. +## @end deftypefn + +## Author: TT <Teresa.Twaroch@ci.tuwien.ac.at> +## Description: CDF of the normal distribution + +function cdf = normal_cdf (x, m, v) + + if (! ((nargin == 1) || (nargin == 3))) + usage ("normal_cdf (x, m, v)"); + endif + + if (nargin == 1) + m = 0; + v = 1; + endif + + if (!isscalar (m) || !isscalar(v)) + [retval, x, m, v] = common_size (x, m, v); + if (retval > 0) + error ("normal_cdf: x, m and v must be of common size or scalar"); + endif + endif + + sz = size (x); + cdf = zeros (sz); + + if (isscalar (m) && isscalar(v)) + if (find (isinf (m) | isnan (m) | !(v >= 0) | !(v < Inf))) + cdf = NaN * ones (sz); + else + cdf = stdnormal_cdf ((x - m) ./ sqrt (v)); + endif + else + k = find (isinf (m) | isnan (m) | !(v >= 0) | !(v < Inf)); + if (any (k)) + cdf(k) = NaN; + endif + + k = find (!isinf (m) & !isnan (m) & (v >= 0) & (v < Inf)); + if (any (k)) + cdf(k) = stdnormal_cdf ((x(k) - m(k)) ./ sqrt (v(k))); + endif + endif + + cdf((v == 0) & (x == m)) = 0.5; + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/norminv.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,79 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} normal_inv (@var{x}, @var{m}, @var{v}) +## For each element of @var{x}, compute the quantile (the inverse of the +## CDF) at @var{x} of the normal distribution with mean @var{m} and +## variance @var{v}. +## +## Default values are @var{m} = 0, @var{v} = 1. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the normal distribution + +function inv = normal_inv (x, m, v) + + if (nargin != 1 && nargin != 3) + usage ("normal_inv (x, m, v)"); + endif + + if (nargin == 1) + m = 0; + v = 1; + endif + + if (!isscalar (m) || !isscalar(v)) + [retval, x, m, v] = common_size (x, m, v); + if (retval > 0) + error ("normal_inv: x, m and v must be of common size or scalars"); + endif + endif + + sz = size (x); + inv = zeros (sz); + + if (isscalar (m) && isscalar(v)) + if (find (isinf (m) | isnan (m) | !(v > 0) | !(v < Inf))) + inv = NaN * ones (sz); + else + inv = m + sqrt (v) .* stdnormal_inv (x); + endif + else + k = find (isinf (m) | isnan (m) | !(v > 0) | !(v < Inf)); + if (any (k)) + inv(k) = NaN; + endif + + k = find (!isinf (m) & !isnan (m) & (v > 0) & (v < Inf)); + if (any (k)) + inv(k) = m(k) + sqrt (v(k)) .* stdnormal_inv (x(k)); + endif + endif + + k = find ((v == 0) & (x > 0) & (x < 1)); + if (any (k)) + inv(k) = m(k); + endif + + inv((v == 0) & (x == 0)) = -Inf; + inv((v == 0) & (x == 1)) = Inf; + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/normpdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,74 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} normal_pdf (@var{x}, @var{m}, @var{v}) +## For each element of @var{x}, compute the probability density function +## (PDF) at @var{x} of the normal distribution with mean @var{m} and +## variance @var{v}. +## +## Default values are @var{m} = 0, @var{v} = 1. +## @end deftypefn + +## Author: TT <Teresa.Twaroch@ci.tuwien.ac.at> +## Description: PDF of the normal distribution + +function pdf = normal_pdf (x, m, v) + + if (nargin != 1 && nargin != 3) + usage ("normal_pdf (x, m, v)"); + endif + + if (nargin == 1) + m = 0; + v = 1; + endif + + if (!isscalar (m) || !isscalar(v)) + [retval, x, m, v] = common_size (x, m, v); + if (retval > 0) + error ("normal_pdf: x, m and v must be of common size or scalars"); + endif + endif + + sz = size (x); + pdf = zeros (sz); + + if (isscalar (m) && isscalar(v)) + if (find (isinf (m) | isnan (m) | !(v >= 0) | !(v < Inf))) + pdf = NaN * ones (sz); + else + pdf = stdnormal_pdf ((x - m) ./ sqrt (v)) ./ sqrt (v); + endif + else + k = find (isinf (m) | isnan (m) | !(v >= 0) | !(v < Inf)); + if (any (k)) + pdf(k) = NaN; + endif + + k = find (!isinf (m) & !isnan (m) & (v >= 0) & (v < Inf)); + if (any (k)) + pdf(k) = stdnormal_pdf ((x(k) - m(k)) ./ sqrt (v(k))) ./ sqrt (v(k)); + endif + endif + + pdf((v == 0) & (x == m)) = Inf; + pdf((v == 0) & ((x < m) | (x > m))) = 0; + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/normrnd.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,92 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} normal_rnd (@var{m}, @var{v}, @var{r}, @var{c}) +## @deftypefnx {Function File} {} normal_rnd (@var{m}, @var{v}, @var{sz}) +## Return an @var{r} by @var{c} or @code{size (@var{sz})} matrix of +## random samples from the normal distribution with parameters @var{m} +## and @var{v}. Both @var{m} and @var{v} 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{m} and @var{v}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the normal distribution + +function rnd = normal_rnd (m, v, r, c) + + if (nargin > 1) + if (!isscalar(m) || !isscalar(v)) + [retval, m, v] = common_size (m, v); + if (retval > 0) + error ("normal_rnd: m and v must be of common size or scalar"); + endif + endif + endif + + if (nargin == 4) + if (! (isscalar (r) && (r > 0) && (r == round (r)))) + error ("normal_rnd: r must be a positive integer"); + endif + if (! (isscalar (c) && (c > 0) && (c == round (c)))) + error ("normal_rnd: c must be a positive integer"); + endif + sz = [r, c]; + + if (any (size (m) != 1) + && (length (size (m)) != length (sz) || any (size (m) != sz))) + error ("normal_rnd: m and v 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 ("normal_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 ("normal_rnd: m and v must be scalar or of size sz"); + endif + elseif (nargin == 2) + sz = size(m); + else + usage ("normal_rnd (m, v, r, c)"); + endif + + if (isscalar (m) && isscalar (v)) + if (find (isnan (m) | isinf (m) | !(v > 0) | !(v < Inf))) + rnd = NaN * ones (sz); + else + rnd = m + sqrt (v) .* randn (sz); + endif + else + rnd = m + sqrt (v) .* randn (sz); + k = find (isnan (m) | isinf (m) | !(v > 0) | !(v < Inf)); + if (any (k)) + rnd(k) = NaN; + endif + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/poisscdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,64 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} poisson_cdf (@var{x}, @var{lambda}) +## For each element of @var{x}, compute the cumulative distribution +## function (CDF) at @var{x} of the Poisson distribution with parameter +## lambda. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the Poisson distribution + +function cdf = poisson_cdf (x, l) + + if (nargin != 2) + usage ("poisson_cdf (x, lambda)"); + endif + + 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 + + cdf = zeros (size (x)); + + k = find (isnan (x) | !(l > 0)); + if (any (k)) + cdf(k) = NaN; + endif + + k = find ((x == Inf) & (l > 0)); + if (any (k)) + cdf(k) = 1; + endif + + k = find ((x >= 0) & (x < Inf) & (l > 0)); + if (any (k)) + 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 + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/poissinv.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,77 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} poisson_inv (@var{x}, @var{lambda}) +## For each component of @var{x}, compute the quantile (the inverse of +## the CDF) at @var{x} of the Poisson distribution with parameter +## @var{lambda}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the Poisson distribution + +function inv = poisson_inv (x, l) + + if (nargin != 2) + usage ("poisson_inv (x, lambda)"); + endif + + 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 + + inv = zeros (size (x)); + + k = find ((x < 0) | (x > 1) | isnan (x) | !(l > 0)); + if (any (k)) + inv(k) = NaN; + endif + + k = find ((x == 1) & (l > 0)); + if (any (k)) + inv(k) = Inf; + endif + + k = find ((x > 0) & (x < 1) & (l > 0)); + if (any (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; + 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 + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/poisspdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,58 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} poisson_pdf (@var{x}, @var{lambda}) +## For each element of @var{x}, compute the probability density function +## (PDF) at @var{x} of the poisson distribution with parameter @var{lambda}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the Poisson distribution + +function pdf = poisson_pdf (x, l) + + if (nargin != 2) + usage ("poisson_pdf (x, lambda)"); + endif + + 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 + + pdf = zeros (size (x)); + + k = find (!(l > 0) | isnan (x)); + if (any (k)) + pdf(k) = NaN; + endif + + k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (l > 0)); + if (any (k)) + 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 + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/poissrnd.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,114 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} poisson_rnd (@var{lambda}, @var{r}, @var{c}) +## Return an @var{r} by @var{c} matrix of random samples from the +## Poisson distribution with parameter @var{lambda}, 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{lambda}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the Poisson distribution + +function rnd = poisson_rnd (l, r, c) + + if (nargin == 3) + if (! (isscalar (r) && (r > 0) && (r == round (r)))) + error ("poisson_rnd: r must be a positive integer"); + endif + if (! (isscalar (c) && (c > 0) && (c == round (c)))) + error ("poisson_rnd: c must be a positive integer"); + endif + 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 == 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 + + 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; + 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 (size (ind))) ./ l(ind)); + num(ind) = num(ind) + 1; + else + break; + endif + endwhile + rnd(k) = num; + endif + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/tcdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,68 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} t_cdf (@var{x}, @var{n}) +## For each element of @var{x}, compute the CDF at @var{x} of the +## t (Student) distribution with @var{n} degrees of freedom, i.e., +## PROB (t(@var{n}) <= @var{x}). +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the t distribution + +function cdf = t_cdf (x, n) + + if (nargin != 2) + usage ("t_cdf (x, n)"); + endif + + 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 + + cdf = zeros (size (x)); + + k = find (isnan (x) | !(n > 0)); + if (any (k)) + cdf(k) = NaN; + endif + + k = find ((x == Inf) & (n > 0)); + if (any (k)) + cdf(k) = 1; + endif + + k = find ((x > -Inf) & (x < Inf) & (n > 0)); + if (any (k)) + 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 + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/tinv.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,83 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} t_inv (@var{x}, @var{n}) +## For each component of @var{x}, compute the quantile (the inverse of +## the CDF) at @var{x} of the t (Student) distribution with parameter +## @var{n}. +## @end deftypefn + +## For very large n, the "correct" formula does not really work well, +## and the quantiles of the standard normal distribution are used +## directly. + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the t distribution + +function inv = t_inv (x, n) + + if (nargin != 2) + usage ("t_inv (x, n)"); + endif + + 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 + + inv = zeros (size (x)); + + k = find ((x < 0) | (x > 1) | isnan (x) | !(n > 0)); + if (any (k)) + inv(k) = NaN; + endif + + k = find ((x == 0) & (n > 0)); + if (any (k)) + inv(k) = -Inf; + endif + + k = find ((x == 1) & (n > 0)); + if (any (k)) + inv(k) = Inf; + endif + + k = find ((x > 0) & (x < 1) & (n > 0) & (n < 10000)); + if (any (k)) + 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 + k = find ((x > 0) & (x < 1) & (n >= 10000)); + if (any (k)) + inv(k) = stdnormal_inv (x(k)); + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/tpdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,61 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} t_pdf (@var{x}, @var{n}) +## For each element of @var{x}, compute the probability density function +## (PDF) at @var{x} of the @var{t} (Student) distribution with @var{n} +## degrees of freedom. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the t distribution + +function pdf = t_pdf (x, n) + + if (nargin != 2) + usage ("t_pdf (x, n)"); + endif + + 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 + + pdf = zeros (size (x)); + + k = find (isnan (x) | !(n > 0) | !(n < Inf)); + if (any (k)) + pdf(k) = NaN; + endif + + k = find (!isinf (x) & !isnan (x) & (n > 0) & (n < Inf)); + if (any (k)) + 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 + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/trnd.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,91 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- 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}. 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{n}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the t distribution + +function rnd = t_rnd (n, r, c) + + if (nargin == 3) + if (! (isscalar (r) && (r > 0) && (r == round (r)))) + error ("t_rnd: r must be a positive integer"); + endif + if (! (isscalar (c) && (c > 0) && (c == round (c)))) + error ("t_rnd: c must be a positive integer"); + endif + 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 == 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 + + 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; + endif + + k = find ((n > 0) & (n < Inf)); + if (any (k)) + rnd(k) = t_inv (rand (size (k)), n(k)); + endif + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/unifcdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,71 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} uniform_cdf (@var{x}, @var{a}, @var{b}) +## Return the CDF at @var{x} of the uniform distribution on [@var{a}, +## @var{b}], i.e., PROB (uniform (@var{a}, @var{b}) <= x). +## +## Default values are @var{a} = 0, @var{b} = 1. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the uniform distribution + +function cdf = uniform_cdf (x, a, b) + + if (nargin != 1 && nargin != 3) + usage ("uniform_cdf (x, a, b)"); + endif + + if (nargin == 1) + a = 0; + b = 1; + endif + + if (!isscalar (a) || !isscalar(b)) + [retval, x, a, b] = common_size (x, a, b); + if (retval > 0) + error ("uniform_cdf: x, a and b must be of common size or scalar"); + endif + endif + + sz = size (x); + cdf = zeros (sz); + + k = find (isnan (x) | !(a < b)); + if (any (k)) + cdf(k) = NaN; + endif + + k = find ((x >= b) & (a < b)); + if (any (k)) + cdf(k) = 1; + endif + + k = find ((x > a) & (x < b)); + if (any (k)) + if (isscalar (a) && isscalar(b)) + cdf(k) = (x(k) < b) .* (x(k) - a) ./ (b - a); + else + cdf(k) = (x(k) < b(k)) .* (x(k) - a(k)) ./ (b(k) - a(k)); + endif + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/unifinv.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,66 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} uniform_inv (@var{x}, @var{a}, @var{b}) +## For each element of @var{x}, compute the quantile (the inverse of the +## CDF) at @var{x} of the uniform distribution on [@var{a}, @var{b}]. +## +## Default values are @var{a} = 0, @var{b} = 1. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the uniform distribution + +function inv = uniform_inv (x, a, b) + + if (nargin != 1 && nargin != 3) + usage ("uniform_inv (x, a, b)"); + endif + + if (nargin == 1) + a = 0; + b = 1; + endif + + if (!isscalar (a) || !isscalar(b)) + [retval, x, a, b] = common_size (x, a, b); + if (retval > 0) + error ("uniform_cdf: x, a and b must be of common size or scalar"); + endif + endif + + sz = size (x); + inv = zeros (sz); + + k = find ((x < 0) | (x > 1) | isnan (x) | !(a < b)); + if (any (k)) + inv(k) = NaN; + endif + + k = find ((x >= 0) & (x <= 1) & (a < b)); + if (any (k)) + if (isscalar (a) && isscalar(b)) + inv(k) = a + x(k) .* (b - a); + else + inv(k) = a(k) + x(k) .* (b(k) - a(k)); + endif + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/unifpdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,66 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} uniform_pdf (@var{x}, @var{a}, @var{b}) +## For each element of @var{x}, compute the PDF at @var{x} of the uniform +## distribution on [@var{a}, @var{b}]. +## +## Default values are @var{a} = 0, @var{b} = 1. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the uniform distribution + +function pdf = uniform_pdf (x, a, b) + + if (nargin != 1 && nargin != 3) + usage ("uniform_pdf (x, a, b)"); + endif + + if (nargin == 1) + a = 0; + b = 1; + endif + + if (!isscalar (a) || !isscalar(b)) + [retval, x, a, b] = common_size (x, a, b); + if (retval > 0) + error ("uniform_pdf: x, a and b must be of common size or scalars"); + endif + endif + + sz = size (x); + pdf = zeros (sz); + + k = find (isnan (x) | !(a < b)); + if (any (k)) + pdf(k) = NaN; + endif + + k = find ((x > a) & (x < b)); + if (any (k)) + if (isscalar (a) && isscalar(b)) + pdf(k) = 1 ./ (b - a); + else + pdf(k) = 1 ./ (b(k) - a(k)); + endif + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/unifrnd.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,92 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} uniform_rnd (@var{a}, @var{b}, @var{r}, @var{c}) +## @deftypefnx {Function File} {} uniform_rnd (@var{a}, @var{b}, @var{sz}) +## Return an @var{r} by @var{c} or a @code{size (@var{sz})} matrix of +## random samples from the uniform distribution on [@var{a}, @var{b}]. +## Both @var{a} and @var{b} 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{a} and @var{b}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Random deviates from the uniform distribution + +function rnd = uniform_rnd (a, b, r, c) + + if (nargin > 1) + if (!isscalar(a) || !isscalar(b)) + [retval, a, b] = common_size (a, b); + if (retval > 0) + error ("uniform_rnd: a and b must be of common size or scalar"); + endif + endif + endif + + if (nargin == 4) + if (! (isscalar (r) && (r > 0) && (r == round (r)))) + error ("uniform_rnd: r must be a positive integer"); + endif + if (! (isscalar (c) && (c > 0) && (c == round (c)))) + error ("uniform_rnd: c must be a positive integer"); + endif + sz = [r, c]; + + if (any (size (a) != 1) + && (length (size (a)) != length (sz) || any (size (a) != sz))) + error ("uniform_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 ("uniform_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 ("uniform_rnd: a and b must be scalar or of size sz"); + endif + elseif (nargin == 2) + sz = size(a); + else + usage ("uniform_rnd (a, b, r, c)"); + endif + + if (isscalar(a) && isscalar(b)) + if (find (!(-Inf < a) | !(a < b) | !(b < Inf))) + rnd = NaN * ones(sz); + else + rnd = a + (b - a) .* rand (sz); + endif + else + rnd = a + (b - a) .* rand (sz); + + k = find (!(-Inf < a) | !(a < b) | !(b < Inf)); + if (any (k)) + rnd(k) = NaN; + endif + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/weibcdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,73 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} weibull_cdf (@var{x}, @var{alpha}, @var{sigma}) +## Compute the cumulative distribution function (CDF) at @var{x} of the +## Weibull distribution with shape parameter @var{alpha} and scale +## parameter @var{sigma}, which is +## +## @example +## 1 - exp(-(x/sigma)^alpha) +## @end example +## +## @noindent +## for @var{x} >= 0. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: CDF of the Weibull distribution + +function cdf = weibull_cdf (x, shape, scale) + + if (nargin != 3) + usage ("weibull_cdf (x, alpha, sigma)"); + endif + + 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 + + cdf = NaN * ones (size (x)); + + ok = ((shape > 0) & (shape < Inf) & (scale > 0) & (scale < Inf)); + + k = find ((x <= 0) & ok); + if (any (k)) + cdf(k) = 0; + endif + + k = find ((x > 0) & (x < Inf) & ok); + if (any (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) = 1; + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/weibinv.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,66 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} weibull_inv (@var{x}, @var{lambda}, @var{alpha}) +## Compute the quantile (the inverse of the CDF) at @var{x} of the +## Weibull distribution with shape parameter @var{alpha} and scale +## parameter @var{sigma}. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: Quantile function of the Weibull distribution + +function inv = weibull_inv (x, shape, scale) + + if (nargin != 3) + usage ("weibull_inv (x, alpha, sigma)"); + endif + + 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 + + inv = NaN * ones (size (x)); + + ok = ((shape > 0) & (shape < Inf) & (scale > 0) & (scale < Inf)); + + k = find ((x == 0) & ok); + if (any (k)) + inv(k) = -Inf; + endif + + k = find ((x > 0) & (x < 1) & ok); + if (any (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; + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/weibpdf.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,71 @@ +## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} weibull_pdf (@var{x}, @var{alpha}, @var{sigma}) +## Compute the probability density function (PDF) at @var{x} of the +## Weibull distribution with shape parameter @var{alpha} and scale +## parameter @var{sigma} which is given by +## +## @example +## alpha * sigma^(-alpha) * x^(alpha-1) * exp(-(x/sigma)^alpha) +## @end example +## +## @noindent +## for @var{x} > 0. +## @end deftypefn + +## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> +## Description: PDF of the Weibull distribution + +function pdf = weibull_pdf (x, shape, scale) + + if (nargin != 3) + usage ("weibull_pdf (x, alpha, sigma)"); + endif + + 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 + + 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) = 0; + endif + + k = find ((x > 0) & (x < Inf) & ok); + if (any (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 + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/wienrnd.m Wed Jul 13 17:43:35 2005 +0000 @@ -0,0 +1,55 @@ +## Copyright (C) 1995, 1996, 1997 Friedrich Leisch +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} wiener_rnd (@var{t}, @var{d}, @var{n}) +## Return a simulated realization of the @var{d}-dimensional Wiener Process +## on the interval [0, @var{t}]. If @var{d} is omitted, @var{d} = 1 is +## used. The first column of the return matrix contains time, the +## remaining columns contain the Wiener process. +## +## The optional parameter @var{n} gives the number of summands used for +## simulating the process over an interval of length 1. If @var{n} is +## omitted, @var{n} = 1000 is used. +## @end deftypefn + +## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at> +## Description: Simulate a Wiener process + +function retval = wiener_rnd (t, d, n) + + if (nargin == 1) + d = 1; + n = 1000; + elseif (nargin == 2) + n = 1000; + elseif (nargin > 3) + 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); + retval = cumsum (retval) / sqrt (n); + + retval = [((1: n*t)' / n), retval]; + +endfunction