Mercurial > hg > octave-lyh

comparison scripts/statistics/distributions/pascal_pdf.m @ 3456:434790acb067

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2000-01-19 06:58:51 by jwe]
author jwe
date Wed, 19 Jan 2000 06:59:23 +0000
parents f8dde1807dee
children 38c61cbf086c
comparison
equal deleted inserted replaced
3455:f758be6e1730 3456:434790acb067
12 ## 12 ##
13 ## You should have received a copy of the GNU General Public License 13 ## You should have received a copy of the GNU General Public License
14 ## along with this file. If not, write to the Free Software Foundation, 14 ## along with this file. If not, write to the Free Software Foundation,
15 ## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. 15 ## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
16 16
17 ## usage: pascal_pdf (x, n, p) 17 ## -*- texinfo -*-
18 ## 18 ## @deftypefn {Function File} {} pascal_pdf (@var{x}, @var{n}, @var{p})
19 ## For each element of x, compute the probability density function (PDF) 19 ## For each element of @var{x}, compute the probability density function
20 ## at x of the Pascal (negative binomial) distribution with parameters n 20 ## (PDF) at @var{x} of the Pascal (negative binomial) distribution with
21 ## and p. 21 ## parameters @var{n} and @var{p}.
22 ## 22 ##
23 ## The number of failures in a Bernoulli experiment with success 23 ## The number of failures in a Bernoulli experiment with success
24 ## probability p before the n-th success follows this distribution. 24 ## probability @var{p} before the @var{n}-th success follows this
25 ## distribution.
26 ## @end deftypefn
25 27
26 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> 28 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
27 ## Description: PDF of the Pascal (negative binomial) distribution 29 ## Description: PDF of the Pascal (negative binomial) distribution
28 30
29 function pdf = pascal_pdf (x, n, p) 31 function pdf = pascal_pdf (x, n, p)
30 32
31 if (nargin != 3) 33 if (nargin != 3)
32 usage ("pascal_pdf (x, n, p)"); 34 usage ("pascal_pdf (x, n, p)");
33 endif 35 endif
34 36
35 [retval, x, n, p] = common_size (x, n, p); 37 [retval, x, n, p] = common_size (x, n, p);
36 if (retval > 0) 38 if (retval > 0)
37 error (["pascal_pdf: ", ... 39 error ("pascal_pdf: x, n and p must be of common size or scalar");
38 "x, n and p must be of common size or scalar"]);
39 endif 40 endif
40 41
41 [r, c] = size (x); 42 [r, c] = size (x);
42 s = r * c; 43 s = r * c;
43 x = reshape (x, 1, s); 44 x = reshape (x, 1, s);
44 n = reshape (n, 1, s); 45 n = reshape (n, 1, s);
45 p = reshape (p, 1, s); 46 p = reshape (p, 1, s);
46 cdf = zeros (1, s); 47 cdf = zeros (1, s);
47 48
48 k = find (isnan (x) | (n < 1) | (n == Inf) | (n != round (n)) ... 49 k = find (isnan (x) | (n < 1) | (n == Inf) | (n != round (n))
49 | (p < 0) | (p > 1)); 50 | (p < 0) | (p > 1));
50 if any (k) 51 if (any (k))
51 pdf(k) = NaN * ones (1, length (k)); 52 pdf(k) = NaN * ones (1, length (k));
52 endif 53 endif
53 54
54 ## Just for the fun of it ... 55 ## Just for the fun of it ...
55 k = find ((x == Inf) & (n > 0) & (n < Inf) & (n == round (n)) ... 56 k = find ((x == Inf) & (n > 0) & (n < Inf) & (n == round (n))
56 & (p == 0)); 57 & (p == 0));
57 if any (k) 58 if (any (k))
58 pdf(k) = ones (1, length (k)); 59 pdf(k) = ones (1, length (k));
59 endif 60 endif
60 61
61 k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (n > 0) ... 62 k = find ((x >= 0) & (x < Inf) & (x == round (x)) & (n > 0)
62 & (n < Inf) & (n == round (n)) & (p > 0) & (p <= 1)); 63 & (n < Inf) & (n == round (n)) & (p > 0) & (p <= 1));
63 if any (k) 64 if (any (k))
64 pdf(k) = bincoeff (-n(k), x(k)) .* (p(k) .^ n(k)) ... 65 pdf(k) = bincoeff (-n(k), x(k)) .* (p(k) .^ n(k)) .* ((p(k) - 1) .^ x(k));
65 .* ((p(k) - 1) .^ x(k));
66 endif 66 endif
67 67
68 pdf = reshape (pdf, r, c); 68 pdf = reshape (pdf, r, c);
69 69
70 endfunction 70 endfunction