Mercurial > octave

comparison scripts/statistics/distributions/binomial_rnd.m @ 3426:f8dde1807dee

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2000-01-13 08:40:00 by jwe]
author jwe
date Thu, 13 Jan 2000 08:40:53 +0000
parents e4f4b2d26ee9
children 434790acb067
comparison
equal deleted inserted replaced
3425:8625164a0a39 3426:f8dde1807dee
1 ## Copyright (C) 1995, 1996, 1997 Kurt Hornik 1 ## Copyright (C) 1995, 1996, 1997 Kurt Hornik
2 ## 2 ##
3 ## This program is free software; you can redistribute it and/or modify 3 ## This program is free software; you can redistribute it and/or modify
4 ## it under the terms of the GNU General Public License as published by 4 ## it under the terms of the GNU General Public License as published by
5 ## the Free Software Foundation; either version 2, or (at your option) 5 ## the Free Software Foundation; either version 2, or (at your option)
6 ## any later version. 6 ## any later version.
7 ## 7 ##
8 ## This program is distributed in the hope that it will be useful, but 8 ## This program is distributed in the hope that it will be useful, but
9 ## WITHOUT ANY WARRANTY; without even the implied warranty of 9 ## WITHOUT ANY WARRANTY; without even the implied warranty of
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 ## General Public License for more details. 11 ## General Public License for more details.
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: binomial_rnd (n, p [, r, c]) 17 ## usage: binomial_rnd (n, p [, r, c])
21 ## matrix is the common size of n and p. 21 ## matrix is the common size of n and p.
22 ## 22 ##
23 ## binomial_rnd (n, p, r, c) returns an r by c matrix of random samples 23 ## binomial_rnd (n, p, r, c) returns an r by c matrix of random samples
24 ## from the binomial distribution with parameters n and p. Both n and p 24 ## from the binomial distribution with parameters n and p. Both n and p
25 ## must be scalar or of size r by c. 25 ## must be scalar or of size r by c.
26 26
27 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at> 27 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
28 ## Description: Random deviates from the binomial distribution 28 ## Description: Random deviates from the binomial distribution
29 29
30 function rnd = binomial_rnd (n, p, r, c) 30 function rnd = binomial_rnd (n, p, r, c)
31 31
37 error ("binomial_rnd: c must be a positive integer"); 37 error ("binomial_rnd: c must be a positive integer");
38 endif 38 endif
39 [retval, n, p] = common_size (n, p, zeros (r, c)); 39 [retval, n, p] = common_size (n, p, zeros (r, c));
40 if (retval > 0) 40 if (retval > 0)
41 error (strcat("binomial_rnd: ", 41 error (strcat("binomial_rnd: ",
42 "n and p must be scalar or of size ", 42 "n and p must be scalar or of size ",
43 sprintf ("%d by %d", r, c))); 43 sprintf ("%d by %d", r, c)));
44 endif 44 endif
45 elseif (nargin == 2) 45 elseif (nargin == 2)
46 [retval, n, p] = common_size (n, p); 46 [retval, n, p] = common_size (n, p);
47 if (retval > 0) 47 if (retval > 0)
48 error ("binomial_rnd: n and p must be of common size or scalar"); 48 error ("binomial_rnd: n and p must be of common size or scalar");
54 [r, c] = size (n); 54 [r, c] = size (n);
55 s = r * c; 55 s = r * c;
56 n = reshape (n, 1, s); 56 n = reshape (n, 1, s);
57 p = reshape (p, 1, s); 57 p = reshape (p, 1, s);
58 rnd = zeros (1, s); 58 rnd = zeros (1, s);
59 59
60 k = find (!(n > 0) | !(n < Inf) | !(n == round (n)) | 60 k = find (!(n > 0) | !(n < Inf) | !(n == round (n)) |
61 !(p <= 0) | !(p >= 1)); 61 !(p <= 0) | !(p >= 1));
62 if any (k) 62 if any (k)
63 rnd(k) = NaN * ones (1, length (k)); 63 rnd(k) = NaN * ones (1, length (k));
64 endif 64 endif
65 65
66 k = find ((n > 0) & (n < Inf) & (n == round (n)) & (p >= 0) & (p <= 1)); 66 k = find ((n > 0) & (n < Inf) & (n == round (n)) & (p >= 0) & (p <= 1));
67 if any (k) 67 if any (k)
68 N = max (n(k)); 68 N = max (n(k));
69 L = length (k); 69 L = length (k);
70 tmp = rand (N, L); 70 tmp = rand (N, L);
71 ind = (1 : N)' * ones (1, L); 71 ind = (1 : N)' * ones (1, L);
72 rnd(k) = sum ((tmp < ones (N, 1) * p(k)) & 72 rnd(k) = sum ((tmp < ones (N, 1) * p(k)) &
73 (ind <= ones (N, 1) * n(k))); 73 (ind <= ones (N, 1) * n(k)));
74 endif 74 endif
75 75
76 rnd = reshape (rnd, r, c); 76 rnd = reshape (rnd, r, c);
77 77
78 endfunction 78 endfunction