octave-maintainers
[Top][All Lists]
Advanced

[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 ----

reply via email to

[Prev in Thread] Current Thread [Next in Thread]