# HG changeset patch # User dbateman # Date 1172240442 0 # Node ID 859f7aaea2540532fd4827ac2fb5ed6fcf67d171 # Parent 1f5de98984c341130be3a050469d0f58eb42c31f [project @ 2007-02-23 14:20:42 by dbateman] diff -r 1f5de98984c3 -r 859f7aaea254 scripts/ChangeLog --- a/scripts/ChangeLog Fri Feb 23 05:24:56 2007 +0000 +++ b/scripts/ChangeLog Fri Feb 23 14:20:42 2007 +0000 @@ -1,3 +1,30 @@ +2007-02-23 David Bateman + + * statistics/distributions/discrete_rnd.m, + statistics/distributions/geornd.m, + statistics/distributions/lognnd.m, + statistics/distributions/nbinrnd.m, + statistics/distributions/wblrnd.m: Accelerate distributions. + + * statistics/distributions/unidcdf.m, + statistics/distributions/unidinv.m, + statistics/distributions/unidpdf.m, + statistics/distributions/unidrnd.m: New functions based on + discrete_cdf, etc. + + * statistics/distributions/pascal_cdf.m, + statistics/distributions/pascal_inv.m, + statistics/distributions/pascal_pdf.m, + statistics/distributions/pascal_rnd.m: Remove. + * statistics/distributions/nbincdf.m, + statistics/distributions/nbininv.m, + statistics/distributions/nbinpdf.m, + statistics/distributions/nbinrnd.m: Replace with matlab + compatible functions. + * deprecated/pascal_cdf.m, deprecated/pascal_inv.m, + deprecated/pascal_pdf.m, deprecated/pascal_rnd.m: Use the new + nbincdf, etc functions to implement these. + 2007-02-22 Daniel J Sebald * plot/__uiobject_draw_axes__.m: Insert newline between plot diff -r 1f5de98984c3 -r 859f7aaea254 scripts/deprecated/pascal_cdf.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/deprecated/pascal_cdf.m Fri Feb 23 14:20:42 2007 +0000 @@ -0,0 +1,37 @@ +## 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} {} pascal_cdf (@var{x}, @var{n}, @var{p}) +## For each element of @var{x}, compute the CDF at x of the Pascal +## (negative binomial) distribution with parameters @var{n} and @var{p}. +## +## The number of failures in a Bernoulli experiment with success +## probability @var{p} before the @var{n}-th success follows this +## distribution. +## @end deftypefn + +## Author: KH +## Description: CDF of the Pascal (negative binomial) distribution + +function cdf = pascal_cdf (varargin) + + cdf = nbincdf(varargin{:}); + +endfunction diff -r 1f5de98984c3 -r 859f7aaea254 scripts/deprecated/pascal_inv.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/deprecated/pascal_inv.m Fri Feb 23 14:20:42 2007 +0000 @@ -0,0 +1,38 @@ +## 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} {} pascal_inv (@var{x}, @var{n}, @var{p}) +## For each element of @var{x}, compute the quantile at @var{x} of the +## Pascal (negative binomial) distribution with parameters @var{n} and +## @var{p}. +## +## The number of failures in a Bernoulli experiment with success +## probability @var{p} before the @var{n}-th success follows this +## distribution. +## @end deftypefn + +## Author: KH +## Description: Quantile function of the Pascal distribution + +function inv = pascal_inv (varargin) + + inv = nbininv(varargin{:}); + +endfunction diff -r 1f5de98984c3 -r 859f7aaea254 scripts/deprecated/pascal_pdf.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/deprecated/pascal_pdf.m Fri Feb 23 14:20:42 2007 +0000 @@ -0,0 +1,38 @@ +## 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} {} pascal_pdf (@var{x}, @var{n}, @var{p}) +## For each element of @var{x}, compute the probability density function +## (PDF) at @var{x} of the Pascal (negative binomial) distribution with +## parameters @var{n} and @var{p}. +## +## The number of failures in a Bernoulli experiment with success +## probability @var{p} before the @var{n}-th success follows this +## distribution. +## @end deftypefn + +## Author: KH +## Description: PDF of the Pascal (negative binomial) distribution + +function pdf = pascal_pdf (varargin) + + pdf = nbinpdf (varargin{:}); + +endfunction diff -r 1f5de98984c3 -r 859f7aaea254 scripts/deprecated/pascal_rnd.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/deprecated/pascal_rnd.m Fri Feb 23 14:20:42 2007 +0000 @@ -0,0 +1,39 @@ +## 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} {} pascal_rnd (@var{n}, @var{p}, @var{r}, @var{c}) +## @deftypefnx {Function File} {} pascal_rnd (@var{n}, @var{p}, @var{sz}) +## Return an @var{r} by @var{c} matrix of random samples from the Pascal +## (negative binomial) distribution with parameters @var{n} and @var{p}. +## Both @var{n} and @var{p} must be scalar or of size @var{r} by @var{c}. +## +## If @var{r} and @var{c} are omitted, the size of the result matrix is +## the common size of @var{n} and @var{p}. Or if @var{sz} is a vector, +## create a matrix of size @var{sz}. +## @end deftypefn + +## Author: KH +## Description: Random deviates from the Pascal distribution + +function rnd = pascal_rnd (varargin) + + rnd = nbinrnd (varargin{:}); + +endfunction diff -r 1f5de98984c3 -r 859f7aaea254 scripts/statistics/distributions/nbincdf.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/nbincdf.m Fri Feb 23 14:20:42 2007 +0000 @@ -0,0 +1,93 @@ +## 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} {} nbincdf (@var{x}, @var{n}, @var{p}) +## For each element of @var{x}, compute the CDF at x of the Pascal +## (negative binomial) distribution with parameters @var{n} and @var{p}. +## +## The number of failures in a Bernoulli experiment with success +## probability @var{p} before the @var{n}-th success follows this +## distribution. +## @end deftypefn + +## Author: KH +## Description: CDF of the Pascal (negative binomial) distribution + +function cdf = nbincdf (x, n, p) + + if (nargin != 3) + print_usage (); + endif + + if (!isscalar(n) || !isscalar(p)) + [retval, x, n, p] = common_size (x, n, p); + if (retval > 0) + error ("nbincdf: x, n and p must be of common size or scalar"); + endif + endif + + cdf = zeros (size (x)); + + k = find (isnan (x) | (n < 1) | (n == Inf) | (n != round (n)) + | (p < 0) | (p > 1)); + if (any (k)) + cdf(k) = NaN; + endif + + k = find ((x == Inf) & (n > 0) & (n < Inf) & (n == round (n)) + & (p >= 0) & (p <= 1)); + if (any (k)) + cdf(k) = 1; + endif + + k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (n > 0) + & (n < Inf) & (n == round (n)) & (p > 0) & (p <= 1)); + if (any (k)) + ## Does anyone know a better way to do the summation? + m = zeros (size (k)); + x = floor (x(k)); + y = cdf(k); + if (isscalar (n) && isscalar (p)) + while (1) + l = find (m <= x); + if (any (l)) + y(l) = y(l) + nbinpdf (m(l), n, p); + m(l) = m(l) + 1; + else + break; + endif + endwhile + else + n = n(k); + p = p(k); + while (1) + l = find (m <= x); + if (any (l)) + y(l) = y(l) + nbinpdf (m(l), n(l), p(l)); + m(l) = m(l) + 1; + else + break; + endif + endwhile + endif + cdf(k) = y; + endif + +endfunction diff -r 1f5de98984c3 -r 859f7aaea254 scripts/statistics/distributions/nbininv.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/nbininv.m Fri Feb 23 14:20:42 2007 +0000 @@ -0,0 +1,94 @@ +## 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} {} nbininv (@var{x}, @var{n}, @var{p}) +## For each element of @var{x}, compute the quantile at @var{x} of the +## Pascal (negative binomial) distribution with parameters @var{n} and +## @var{p}. +## +## The number of failures in a Bernoulli experiment with success +## probability @var{p} before the @var{n}-th success follows this +## distribution. +## @end deftypefn + +## Author: KH +## Description: Quantile function of the Pascal distribution + +function inv = nbininv (x, n, p) + + if (nargin != 3) + print_usage (); + endif + + if (!isscalar(n) || !isscalar(p)) + [retval, x, n, p] = common_size (x, n, p); + if (retval > 0) + error ("nbininv: x, n and p must be of common size or scalar"); + endif + endif + + inv = zeros (size (x)); + + k = find (isnan (x) | (x < 0) | (x > 1) | (n < 1) | (n == Inf) + | (n != round (n)) | (p < 0) | (p > 1)); + if (any (k)) + inv(k) = NaN; + endif + + k = find ((x == 1) & (n > 0) & (n < Inf) & (n == round (n)) + & (p >= 0) & (p <= 1)); + if (any (k)) + inv(k) = Inf; + endif + + k = find ((x >= 0) & (x < 1) & (n > 0) & (n < Inf) + & (n == round (n)) & (p > 0) & (p <= 1)); + if (any (k)) + m = zeros (size (k)); + x = x(k); + if (isscalar (n) && isscalar (p)) + s = p ^ n * ones (size(k)); + while (1) + l = find (s < x); + if (any (l)) + m(l) = m(l) + 1; + s(l) = s(l) + nbinpdf (m(l), n, p); + else + break; + endif + endwhile + else + n = n(k); + p = p(k); + s = p .^ n; + while (1) + l = find (s < x); + if (any (l)) + m(l) = m(l) + 1; + s(l) = s(l) + nbinpdf (m(l), n(l), p(l)); + else + break; + endif + endwhile + endif + inv(k) = m; + endif + +endfunction diff -r 1f5de98984c3 -r 859f7aaea254 scripts/statistics/distributions/nbinpdf.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/nbinpdf.m Fri Feb 23 14:20:42 2007 +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} {} nbinpdf (@var{x}, @var{n}, @var{p}) +## For each element of @var{x}, compute the probability density function +## (PDF) at @var{x} of the Pascal (negative binomial) distribution with +## parameters @var{n} and @var{p}. +## +## The number of failures in a Bernoulli experiment with success +## probability @var{p} before the @var{n}-th success follows this +## distribution. +## @end deftypefn + +## Author: KH +## Description: PDF of the Pascal (negative binomial) distribution + +function pdf = nbinpdf (x, n, p) + + if (nargin != 3) + print_usage (); + endif + + if (!isscalar(n) || !isscalar(p)) + [retval, x, n, p] = common_size (x, n, p); + if (retval > 0) + error ("nbinpdf: x, n and p must be of common size or scalar"); + endif + endif + + pdf = zeros (size (x)); + + k = find (isnan (x) | (n < 1) | (n == Inf) | (n != round (n)) + | (p < 0) | (p > 1)); + if (any (k)) + pdf(k) = NaN; + endif + + ## Just for the fun of it ... + k = find ((x == Inf) & (n > 0) & (n < Inf) & (n == round (n)) + & (p == 0)); + if (any (k)) + pdf(k) = 1; + endif + + k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (n > 0) + & (n < Inf) & (n == round (n)) & (p > 0) & (p <= 1)); + if (any (k)) + if (isscalar (n) && isscalar (p)) + pdf(k) = bincoeff (-n, x(k)) .* (p ^ n) .* ((p - 1) .^ x(k)); + else + pdf(k) = bincoeff (-n(k), x(k)) .* (p(k) .^ n(k)) .* ((p(k) - 1) .^ x(k)); + endif + endif + +endfunction diff -r 1f5de98984c3 -r 859f7aaea254 scripts/statistics/distributions/nbinrnd.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/statistics/distributions/nbinrnd.m Fri Feb 23 14:20:42 2007 +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} {} nbinrnd (@var{n}, @var{p}, @var{r}, @var{c}) +## @deftypefnx {Function File} {} nbinrnd (@var{n}, @var{p}, @var{sz}) +## Return an @var{r} by @var{c} matrix of random samples from the Pascal +## (negative binomial) distribution with parameters @var{n} and @var{p}. +## Both @var{n} and @var{p} must be scalar or of size @var{r} by @var{c}. +## +## If @var{r} and @var{c} are omitted, the size of the result matrix is +## the common size of @var{n} and @var{p}. Or if @var{sz} is a vector, +## create a matrix of size @var{sz}. +## @end deftypefn + +## Author: KH +## Description: Random deviates from the Pascal distribution + +function rnd = nbinrnd (n, p, r, c) + + if (nargin > 1) + if (!isscalar(n) || !isscalar(p)) + [retval, n, p] = common_size (n, p); + if (retval > 0) + error ("nbinrnd: 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 ("nbinrnd: r must be a positive integer"); + endif + if (! (isscalar (c) && (c > 0) && (c == round (c)))) + error ("nbinrnd: 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 ("nbinrnd: n and p must be scalar or of size [r, c]"); + endif + + elseif (nargin == 3) + if (isscalar (r) && (r > 0)) + sz = [r, r]; + elseif (isvector(r) && all (r > 0)) + sz = r(:)'; + else + error ("nbinrnd: r must be a postive integer or vector"); + endif + + if (any (size (n) != 1) && + ((length (size (n)) != length (sz)) || any (size (n) != sz))) + error ("nbinrnd: n and p must be scalar or of size sz"); + endif + elseif (nargin == 2) + sz = size(n); + else + print_usage (); + endif + + if (isscalar (n) && isscalar (p)) + if ((n < 1) || (n == Inf) || (n != round (n)) || (p <= 0) || (p > 1)); + rnd = NaN * ones (sz); + elseif ((n > 0) && (n < Inf) && (n == round (n)) && + (p > 0) && (p <= 1)) + rnd = randp ((1 - p) ./ p .* randg (n, sz)); + else + rnd = zeros (sz); + endif + else + rnd = zeros (sz); + + k = find ((n < 1) || (n == Inf) || (n != round (n)) || + (p <= 0) || (p > 1)); + if (any (k)) + rnd(k) = NaN; + endif + + k = find ((n > 0) & (n < Inf) & (n == round (n)) & (p > 0) & (p <= 1)); + if (any (k)) + rnd(k) = randp ((1 - p(k)) ./ p(k) .* randg (n(k), size(k))); + endif + endif + +endfunction diff -r 1f5de98984c3 -r 859f7aaea254 scripts/statistics/distributions/pascal_cdf.m --- a/scripts/statistics/distributions/pascal_cdf.m Fri Feb 23 05:24:56 2007 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,93 +0,0 @@ -## 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} {} pascal_cdf (@var{x}, @var{n}, @var{p}) -## For each element of @var{x}, compute the CDF at x of the Pascal -## (negative binomial) distribution with parameters @var{n} and @var{p}. -## -## The number of failures in a Bernoulli experiment with success -## probability @var{p} before the @var{n}-th success follows this -## distribution. -## @end deftypefn - -## Author: KH -## Description: CDF of the Pascal (negative binomial) distribution - -function cdf = pascal_cdf (x, n, p) - - if (nargin != 3) - print_usage (); - endif - - if (!isscalar(n) || !isscalar(p)) - [retval, x, n, p] = common_size (x, n, p); - if (retval > 0) - error ("pascal_cdf: x, n and p must be of common size or scalar"); - endif - endif - - cdf = zeros (size (x)); - - k = find (isnan (x) | (n < 1) | (n == Inf) | (n != round (n)) - | (p < 0) | (p > 1)); - if (any (k)) - cdf(k) = NaN; - endif - - k = find ((x == Inf) & (n > 0) & (n < Inf) & (n == round (n)) - & (p >= 0) & (p <= 1)); - if (any (k)) - cdf(k) = 1; - endif - - k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (n > 0) - & (n < Inf) & (n == round (n)) & (p > 0) & (p <= 1)); - if (any (k)) - ## Does anyone know a better way to do the summation? - m = zeros (size (k)); - x = floor (x(k)); - y = cdf(k); - if (isscalar (n) && isscalar (p)) - while (1) - l = find (m <= x); - if (any (l)) - y(l) = y(l) + pascal_pdf (m(l), n, p); - m(l) = m(l) + 1; - else - break; - endif - endwhile - else - n = n(k); - p = p(k); - while (1) - l = find (m <= x); - if (any (l)) - y(l) = y(l) + pascal_pdf (m(l), n(l), p(l)); - m(l) = m(l) + 1; - else - break; - endif - endwhile - endif - cdf(k) = y; - endif - -endfunction diff -r 1f5de98984c3 -r 859f7aaea254 scripts/statistics/distributions/pascal_inv.m --- a/scripts/statistics/distributions/pascal_inv.m Fri Feb 23 05:24:56 2007 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,94 +0,0 @@ -## 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} {} pascal_inv (@var{x}, @var{n}, @var{p}) -## For each element of @var{x}, compute the quantile at @var{x} of the -## Pascal (negative binomial) distribution with parameters @var{n} and -## @var{p}. -## -## The number of failures in a Bernoulli experiment with success -## probability @var{p} before the @var{n}-th success follows this -## distribution. -## @end deftypefn - -## Author: KH -## Description: Quantile function of the Pascal distribution - -function inv = pascal_inv (x, n, p) - - if (nargin != 3) - print_usage (); - endif - - if (!isscalar(n) || !isscalar(p)) - [retval, x, n, p] = common_size (x, n, p); - if (retval > 0) - error ("pascal_inv: x, n and p must be of common size or scalar"); - endif - endif - - inv = zeros (size (x)); - - k = find (isnan (x) | (x < 0) | (x > 1) | (n < 1) | (n == Inf) - | (n != round (n)) | (p < 0) | (p > 1)); - if (any (k)) - inv(k) = NaN; - endif - - k = find ((x == 1) & (n > 0) & (n < Inf) & (n == round (n)) - & (p >= 0) & (p <= 1)); - if (any (k)) - inv(k) = Inf; - endif - - k = find ((x >= 0) & (x < 1) & (n > 0) & (n < Inf) - & (n == round (n)) & (p > 0) & (p <= 1)); - if (any (k)) - m = zeros (size (k)); - x = x(k); - if (isscalar (n) && isscalar (p)) - s = p ^ n * ones (size(k)); - while (1) - l = find (s < x); - if (any (l)) - m(l) = m(l) + 1; - s(l) = s(l) + pascal_pdf (m(l), n, p); - else - break; - endif - endwhile - else - n = n(k); - p = p(k); - s = p .^ n; - while (1) - l = find (s < x); - if (any (l)) - m(l) = m(l) + 1; - s(l) = s(l) + pascal_pdf (m(l), n(l), p(l)); - else - break; - endif - endwhile - endif - inv(k) = m; - endif - -endfunction diff -r 1f5de98984c3 -r 859f7aaea254 scripts/statistics/distributions/pascal_pdf.m --- a/scripts/statistics/distributions/pascal_pdf.m Fri Feb 23 05:24:56 2007 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,72 +0,0 @@ -## 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} {} pascal_pdf (@var{x}, @var{n}, @var{p}) -## For each element of @var{x}, compute the probability density function -## (PDF) at @var{x} of the Pascal (negative binomial) distribution with -## parameters @var{n} and @var{p}. -## -## The number of failures in a Bernoulli experiment with success -## probability @var{p} before the @var{n}-th success follows this -## distribution. -## @end deftypefn - -## Author: KH -## Description: PDF of the Pascal (negative binomial) distribution - -function pdf = pascal_pdf (x, n, p) - - if (nargin != 3) - print_usage (); - endif - - if (!isscalar(n) || !isscalar(p)) - [retval, x, n, p] = common_size (x, n, p); - if (retval > 0) - error ("pascal_pdf: x, n and p must be of common size or scalar"); - endif - endif - - pdf = zeros (size (x)); - - k = find (isnan (x) | (n < 1) | (n == Inf) | (n != round (n)) - | (p < 0) | (p > 1)); - if (any (k)) - pdf(k) = NaN; - endif - - ## Just for the fun of it ... - k = find ((x == Inf) & (n > 0) & (n < Inf) & (n == round (n)) - & (p == 0)); - if (any (k)) - pdf(k) = 1; - endif - - k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (n > 0) - & (n < Inf) & (n == round (n)) & (p > 0) & (p <= 1)); - if (any (k)) - if (isscalar (n) && isscalar (p)) - pdf(k) = bincoeff (-n, x(k)) .* (p ^ n) .* ((p - 1) .^ x(k)); - else - pdf(k) = bincoeff (-n(k), x(k)) .* (p(k) .^ n(k)) .* ((p(k) - 1) .^ x(k)); - endif - endif - -endfunction diff -r 1f5de98984c3 -r 859f7aaea254 scripts/statistics/distributions/pascal_rnd.m --- a/scripts/statistics/distributions/pascal_rnd.m Fri Feb 23 05:24:56 2007 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,120 +0,0 @@ -## 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} {} pascal_rnd (@var{n}, @var{p}, @var{r}, @var{c}) -## @deftypefnx {Function File} {} pascal_rnd (@var{n}, @var{p}, @var{sz}) -## Return an @var{r} by @var{c} matrix of random samples from the Pascal -## (negative binomial) distribution with parameters @var{n} and @var{p}. -## Both @var{n} and @var{p} must be scalar or of size @var{r} by @var{c}. -## -## If @var{r} and @var{c} are omitted, the size of the result matrix is -## the common size of @var{n} and @var{p}. Or if @var{sz} is a vector, -## create a matrix of size @var{sz}. -## @end deftypefn - -## Author: KH -## Description: Random deviates from the Pascal distribution - -function rnd = pascal_rnd (n, p, r, c) - - if (nargin > 1) - if (!isscalar(n) || !isscalar(p)) - [retval, n, p] = common_size (n, p); - if (retval > 0) - error ("pascal_rnd: n and p must be of common size or scalar"); - endif - endif - endif - - if (nargin == 4) - if (! (isscalar (r) && (r > 0) && (r == round (r)))) - error ("pascal_rnd: r must be a positive integer"); - endif - if (! (isscalar (c) && (c > 0) && (c == round (c)))) - error ("pascal_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 ("pascal_rnd: n and p must be scalar or of size [r, c]"); - endif - - elseif (nargin == 3) - if (isscalar (r) && (r > 0)) - sz = [r, r]; - elseif (isvector(r) && all (r > 0)) - sz = r(:)'; - else - error ("pascal_rnd: r must be a postive integer or vector"); - endif - - if (any (size (n) != 1) && - ((length (size (n)) != length (sz)) || any (size (n) != sz))) - error ("pascal_rnd: n and p must be scalar or of size sz"); - endif - elseif (nargin == 2) - sz = size(n); - else - print_usage (); - endif - - if (isscalar (n) && isscalar (p)) - if ((n < 1) || (n == Inf) || (n != round (n)) || (p <= 0) || (p > 1)); - rnd = NaN * ones (sz) - elseif ((n > 0) && (n < Inf) && (n == round (n)) && - (p > 0) && (p <= 1)) - L = prod (sz); - tmp = floor (log (rand (n, L)) / log (1 - p)); - if (n == 1) - rnd = tmp; - else - ind = (1 : n)' * ones (1, L); - rnd = sum (tmp .* (ind <= n)); - endif - else - rnd = zeros (sz); - endif - else - rnd = zeros (sz); - - k = find ((n < 1) || (n == Inf) || (n != round (n)) || - (p <= 0) || (p > 1)); - if (any (k)) - rnd(k) = NaN; - endif - - k = find ((n > 0) & (n < Inf) & (n == round (n)) & (p > 0) & (p <= 1)); - if (any (k)) - n = reshape (n, 1, prod (sz)); - p = reshape (p, 1, prod (sz)); - N = max (n(k)); - L = length (k); - tmp = floor (log (rand (N, L)) ./ (ones (N, 1) * log (1 - p(k)))); - if (N == 1) - rnd(k) = tmp; - else - ind = (1 : N)' * ones (1, L); - rnd(k) = sum (tmp .* (ind <= ones (N, 1) * n(k))); - endif - endif - endif - -endfunction