1 ## Copyright (C) 2012 Rik Wehbring
2 ## Copyright (C) 1997-2012 Kurt Hornik
4 ## This file is part of Octave.
6 ## Octave is free software; you can redistribute it and/or modify it
7 ## under the terms of the GNU General Public License as published by
8 ## the Free Software Foundation; either version 3 of the License, or (at
9 ## your option) any later version.
11 ## Octave is distributed in the hope that it will be useful, but
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 ## General Public License for more details.
16 ## You should have received a copy of the GNU General Public License
17 ## along with Octave; see the file COPYING. If not, see
18 ## <http://www.gnu.org/licenses/>.
21 ## @deftypefn {Function File} {} hygernd (@var{t}, @var{m}, @var{n})
22 ## @deftypefnx {Function File} {} hygernd (@var{t}, @var{m}, @var{n}, @var{r})
23 ## @deftypefnx {Function File} {} hygernd (@var{t}, @var{m}, @var{n}, @var{r}, @var{c}, @dots{})
24 ## @deftypefnx {Function File} {} hygernd (@var{t}, @var{m}, @var{n}, [@var{sz}])
25 ## Return a matrix of random samples from the hypergeometric distribution
26 ## with parameters @var{t}, @var{m}, and @var{n}.
28 ## The parameters @var{t}, @var{m}, and @var{n} must be positive integers
29 ## with @var{m} and @var{n} not greater than @var{t}.
31 ## When called with a single size argument, return a square matrix with
32 ## the dimension specified. When called with more than one scalar argument the
33 ## first two arguments are taken as the number of rows and columns and any
34 ## further arguments specify additional matrix dimensions. The size may also
35 ## be specified with a vector of dimensions @var{sz}.
37 ## If no size arguments are given then the result matrix is the common size of
38 ## @var{t}, @var{m}, and @var{n}.
41 function rnd = hygernd (t, m, n, varargin)
47 if (! isscalar (t) || ! isscalar (m) || ! isscalar (n))
48 [retval, t, m, n] = common_size (t, m, n);
50 error ("hygernd: T, M, and N must be of common size or scalars");
54 if (iscomplex (t) || iscomplex (m) || iscomplex (n))
55 error ("hygernd: T, M, and N must not be complex");
61 if (isscalar (varargin{1}) && varargin{1} >= 0)
62 sz = [varargin{1}, varargin{1}];
63 elseif (isrow (varargin{1}) && all (varargin{1} >= 0))
66 error ("hygernd: dimension vector must be row vector of non-negative integers");
69 if (any (cellfun (@(x) (!isscalar (x) || x < 0), varargin)))
70 error ("hygernd: dimensions must be non-negative integers");
75 if (!isscalar (t) && !isequal (size (t), sz))
76 error ("hygernd: T, M, and N must be scalar or of size SZ");
79 if (isa (t, "single") || isa (m, "single") || isa (n, "single"))
85 ok = ((t >= 0) & (m >= 0) & (n > 0) & (m <= t) & (n <= t) &
86 (t == fix (t)) & (m == fix (m)) & (n == fix (n)));
91 p = hygepdf (v, t, m, n);
92 rnd = v(lookup (cumsum (p(1:end-1)) / sum (p), rand (sz)) + 1);
93 rnd = reshape (rnd, sz);
94 if (strcmp (cls, "single"))
103 for i = find (ok(:)') # Must be row vector arg to for loop
105 p = hygepdf (v, t(i), m(i), n(i));
106 rnd(i) = v(lookup (cumsum (p(1 : end-1)) / sum (p), rn(i)) + 1);
113 %!assert(size (hygernd (4,2,2)), [1, 1]);
114 %!assert(size (hygernd (4*ones(2,1), 2,2)), [2, 1]);
115 %!assert(size (hygernd (4*ones(2,2), 2,2)), [2, 2]);
116 %!assert(size (hygernd (4, 2*ones(2,1), 2)), [2, 1]);
117 %!assert(size (hygernd (4, 2*ones(2,2), 2)), [2, 2]);
118 %!assert(size (hygernd (4, 2, 2*ones(2,1))), [2, 1]);
119 %!assert(size (hygernd (4, 2, 2*ones(2,2))), [2, 2]);
120 %!assert(size (hygernd (4, 2, 2, 3)), [3, 3]);
121 %!assert(size (hygernd (4, 2, 2, [4 1])), [4, 1]);
122 %!assert(size (hygernd (4, 2, 2, 4, 1)), [4, 1]);
124 %!assert(class (hygernd (4,2,2)), "double");
125 %!assert(class (hygernd (single(4),2,2)), "single");
126 %!assert(class (hygernd (single([4 4]),2,2)), "single");
127 %!assert(class (hygernd (4,single(2),2)), "single");
128 %!assert(class (hygernd (4,single([2 2]),2)), "single");
129 %!assert(class (hygernd (4,2,single(2))), "single");
130 %!assert(class (hygernd (4,2,single([2 2]))), "single");
132 %% Test input validation
135 %!error hygernd (1,2)
136 %!error hygernd (ones(3),ones(2),ones(2), 2)
137 %!error hygernd (ones(2),ones(3),ones(2), 2)
138 %!error hygernd (ones(2),ones(2),ones(3), 2)
139 %!error hygernd (i, 2, 2)
140 %!error hygernd (2, i, 2)
141 %!error hygernd (2, 2, i)
142 %!error hygernd (4,2,2, -1)
143 %!error hygernd (4,2,2, ones(2))
144 %!error hygernd (4,2,2, [2 -1 2])
145 %!error hygernd (4*ones(2),2,2, 3)
146 %!error hygernd (4*ones(2),2,2, [3, 2])
147 %!error hygernd (4*ones(2),2,2, 3, 2)