[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Low hanging fruit - Accelerated random distribution functions
From: |
David Bateman |
Subject: |
Re: Low hanging fruit - Accelerated random distribution functions |
Date: |
Fri, 23 Feb 2007 11:21:30 +0100 |
User-agent: |
Thunderbird 1.5.0.7 (X11/20060921) |
Paul Kienzle wrote:
> On Feb 22, 2007, at 7:29 PM, Paul Kienzle wrote:
>
>> On Feb 22, 2007, at 11:46 AM, David Bateman wrote:
>>
>>> pascal_rnd,
>>
>> Needs work.
>>
>
> Looking at the heart of pascal, it is a poisson
> variate which a lambda chosen at random from a
> beta distribution:
>
> function r = nbinrnd(n,p,varargin)
> r = randp( (1-p)./p .* randg(n,varargin{:}) );
> endfunction
>
> I renamed it nbinrnd for compatibility.
>
> - Paul
>
>
Ok, then consider the attached patch to have replace pascal_* with nbin*
and adapt nbinrnd accordingly. Note that nbinrnd returns NaN for invalid
n and p so it needs to be slightly more complex that the code you proposed..
D.
--
David Bateman address@hidden
Motorola Labs - Paris +33 1 69 35 48 04 (Ph)
Parc Les Algorithmes, Commune de St Aubin +33 6 72 01 06 33 (Mob)
91193 Gif-Sur-Yvette FRANCE +33 1 69 35 77 01 (Fax)
The information contained in this communication has been classified as:
[x] General Business Information
[ ] Motorola Internal Use Only
[ ] Motorola Confidential Proprietary
*** ./scripts/deprecated/pascal_inv.m.orig44 2007-02-23 11:00:04.841353812
+0100
--- ./scripts/deprecated/pascal_inv.m 2007-02-23 10:59:21.634102189 +0100
***************
*** 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 <address@hidden>
+ ## Description: Quantile function of the Pascal distribution
+
+ function inv = pascal_inv (varargin)
+
+ inv = nbininv(varargin{:});
+
+ endfunction
*** ./scripts/deprecated/pascal_pdf.m.orig44 2007-02-23 11:00:04.841353812
+0100
--- ./scripts/deprecated/pascal_pdf.m 2007-02-23 10:58:41.260683789 +0100
***************
*** 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 <address@hidden>
+ ## Description: PDF of the Pascal (negative binomial) distribution
+
+ function pdf = pascal_pdf (varargin)
+
+ pdf = nbinpdf (varargin{:});
+
+ endfunction
*** ./scripts/deprecated/pascal_rnd.m.orig44 2007-02-23 11:00:04.841353812
+0100
--- ./scripts/deprecated/pascal_rnd.m 2007-02-23 10:56:43.502298583 +0100
***************
*** 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 <address@hidden>
+ ## Description: Random deviates from the Pascal distribution
+
+ function rnd = pascal_rnd (varargin)
+
+ rnd = nbinrnd (varargin{:});
+
+ endfunction
*** ./scripts/deprecated/pascal_cdf.m.orig44 2007-02-23 10:59:46.055547005
+0100
--- ./scripts/deprecated/pascal_cdf.m 2007-02-23 10:56:51.743761115 +0100
***************
*** 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 <address@hidden>
+ ## Description: CDF of the Pascal (negative binomial) distribution
+
+ function cdf = pascal_cdf (varargin)
+
+ cdf = nbincdf(varargin{:});
+
+ endfunction
*** ./scripts/statistics/distributions/pascal_inv.m.orig44 2006-10-10
18:10:30.000000000 +0200
--- ./scripts/statistics/distributions/pascal_inv.m 2007-02-23
10:50:53.366762623 +0100
***************
*** 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} {} 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 <address@hidden>
- ## 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
--- 0 ----
*** ./scripts/statistics/distributions/nbinpdf.m.orig44 2007-02-23
10:51:30.314378940 +0100
--- ./scripts/statistics/distributions/nbinpdf.m 2007-02-23
10:54:03.019927867 +0100
***************
*** 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 <address@hidden>
+ ## 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
*** ./scripts/statistics/distributions/nbininv.m.orig44 2007-02-23
10:50:53.366762623 +0100
--- ./scripts/statistics/distributions/nbininv.m 2007-02-23
10:54:15.869064159 +0100
***************
*** 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 <address@hidden>
+ ## 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
*** ./scripts/statistics/distributions/nbincdf.m.orig44 2007-02-23
10:49:21.146782780 +0100
--- ./scripts/statistics/distributions/nbincdf.m 2007-02-23
10:53:46.154065326 +0100
***************
*** 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 <address@hidden>
+ ## 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
*** ./scripts/statistics/distributions/pascal_rnd.m.orig44 2006-10-10
18:10:30.000000000 +0200
--- ./scripts/statistics/distributions/pascal_rnd.m 2007-02-23
10:49:43.806293556 +0100
***************
*** 1,120 ****
- ## 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 <address@hidden>
- ## 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
--- 0 ----
*** ./scripts/statistics/distributions/nbinrnd.m.orig44 2007-02-23
10:50:04.832917661 +0100
--- ./scripts/statistics/distributions/nbinrnd.m 2007-02-23
11:18:18.841197853 +0100
***************
*** 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 <address@hidden>
+ ## 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
*** ./scripts/statistics/distributions/pascal_cdf.m.orig44 2006-10-10
18:10:30.000000000 +0200
--- ./scripts/statistics/distributions/pascal_cdf.m 2007-02-23
10:50:29.588304867 +0100
***************
*** 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} {} 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 <address@hidden>
- ## 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
--- 0 ----
*** ./scripts/statistics/distributions/pascal_pdf.m.orig44 2006-10-10
18:10:30.000000000 +0200
--- ./scripts/statistics/distributions/pascal_pdf.m 2007-02-23
10:51:30.314378940 +0100
***************
*** 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} {} 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 <address@hidden>
- ## 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
--- 0 ----
- Re: Low hanging fruit - Accelerated random distribution functions, (continued)
- Re: Low hanging fruit - Accelerated random distribution functions, Paul Kienzle, 2007/02/22
- Re: Low hanging fruit - Accelerated random distribution functions, Paul Kienzle, 2007/02/22
- Re: Low hanging fruit - Accelerated random distribution functions, David Bateman, 2007/02/23
- Re: Low hanging fruit - Accelerated random distribution functions, Paul Kienzle, 2007/02/23
- Re: Low hanging fruit - Accelerated random distribution functions, David Bateman, 2007/02/23
- Re: Low hanging fruit - Accelerated random distribution functions, Paul Kienzle, 2007/02/23
- Re: Low hanging fruit - Accelerated random distribution functions, David Bateman, 2007/02/23
- Re: Low hanging fruit - Accelerated random distribution functions, David Bateman, 2007/02/23
Re: Low hanging fruit - Accelerated random distribution functions, Paul Kienzle, 2007/02/22
Re: Low hanging fruit - Accelerated random distribution functions, David Bateman, 2007/02/23
- Re: Low hanging fruit - Accelerated random distribution functions, Daniel J Sebald, 2007/02/23
- Re: Low hanging fruit - Accelerated random distribution functions, David Bateman, 2007/02/23
- Re: Low hanging fruit - Accelerated random distribution functions, Daniel J Sebald, 2007/02/23
- Re: Low hanging fruit - Accelerated random distribution functions, David Bateman, 2007/02/24
- Re: Low hanging fruit - Accelerated random distribution functions, Daniel J Sebald, 2007/02/24
- Re: Low hanging fruit - Accelerated random distribution functions, David Bateman, 2007/02/24
Re: Low hanging fruit - Accelerated random distribution functions, Daniel J Sebald, 2007/02/24