1## Copyright (C) 2012 Rik Wehbring
2## Copyright (C) 1995-2016 Kurt Hornik
3##
4## This program is free software: you can redistribute it and/or
5## modify it under the terms of the GNU General Public License as
6## published by the Free Software Foundation, either version 3 of the
7## License, or (at your option) any later version.
8##
9## This program is distributed in the hope that it will be useful, but
10## WITHOUT ANY WARRANTY; without even the implied warranty of
11## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
12## General Public License for more details.
13##
14## You should have received a copy of the GNU General Public License
15## along with this program; see the file COPYING.  If not, see
16## <http://www.gnu.org/licenses/>.
17
18## -*- texinfo -*-
19## @deftypefn  {} {} betarnd (@var{a}, @var{b})
20## @deftypefnx {} {} betarnd (@var{a}, @var{b}, @var{r})
21## @deftypefnx {} {} betarnd (@var{a}, @var{b}, @var{r}, @var{c}, @dots{})
22## @deftypefnx {} {} betarnd (@var{a}, @var{b}, [@var{sz}])
23## Return a matrix of random samples from the Beta distribution with parameters
24## @var{a} and @var{b}.
25##
26## When called with a single size argument, return a square matrix with
27## the dimension specified.  When called with more than one scalar argument the
28## first two arguments are taken as the number of rows and columns and any
29## further arguments specify additional matrix dimensions.  The size may also
30## be specified with a vector of dimensions @var{sz}.
31##
32## If no size arguments are given then the result matrix is the common size of
33## @var{a} and @var{b}.
34## @end deftypefn
35
36## Author: KH <Kurt.Hornik@wu-wien.ac.at>
37## Description: Random deviates from the Beta distribution
38
39function rnd = betarnd (a, b, varargin)
40
41  if (nargin < 2)
42    print_usage ();
43  endif
44
45  if (! isscalar (a) || ! isscalar (b))
46    [retval, a, b] = common_size (a, b);
47    if (retval > 0)
48      error ("betarnd: A and B must be of common size or scalars");
49    endif
50  endif
51
52  if (iscomplex (a) || iscomplex (b))
53    error ("betarnd: A and B must not be complex");
54  endif
55
56  if (nargin == 2)
57    sz = size (a);
58  elseif (nargin == 3)
59    if (isscalar (varargin{1}) && varargin{1} >= 0)
60      sz = [varargin{1}, varargin{1}];
61    elseif (isrow (varargin{1}) && all (varargin{1} >= 0))
62      sz = varargin{1};
63    else
64      error ("betarnd: dimension vector must be row vector of non-negative integers");
65    endif
66  elseif (nargin > 3)
67    if (any (cellfun (@(x) (! isscalar (x) || x < 0), varargin)))
68      error ("betarnd: dimensions must be non-negative integers");
69    endif
70    sz = [varargin{:}];
71  endif
72
73  if (! isscalar (a) && ! isequal (size (a), sz))
74    error ("betarnd: A and B must be scalar or of size SZ");
75  endif
76
77  if (isa (a, "single") || isa (b, "single"))
78    cls = "single";
79  else
80    cls = "double";
81  endif
82
83  if (isscalar (a) && isscalar (b))
84    if ((a > 0) && (a < Inf) && (b > 0) && (b < Inf))
85      r = randg (a, sz, cls);
86      rnd = r ./ (r + randg (b, sz, cls));
87    else
88      rnd = NaN (sz, cls);
89    endif
90  else
91    rnd = NaN (sz, cls);
92
93    k = (a > 0) & (a < Inf) & (b > 0) & (b < Inf);
94    r = randg (a(k), cls);
95    rnd(k) = r ./ (r + randg (b(k), cls));
96  endif
97
98endfunction
99
100
101%!assert (size (betarnd (1,2)), [1, 1])
102%!assert (size (betarnd (ones (2,1), 2)), [2, 1])
103%!assert (size (betarnd (ones (2,2), 2)), [2, 2])
104%!assert (size (betarnd (1, 2*ones (2,1))), [2, 1])
105%!assert (size (betarnd (1, 2*ones (2,2))), [2, 2])
106%!assert (size (betarnd (1, 2, 3)), [3, 3])
107%!assert (size (betarnd (1, 2, [4 1])), [4, 1])
108%!assert (size (betarnd (1, 2, 4, 1)), [4, 1])
109
110## Test class of input preserved
111%!assert (class (betarnd (1, 2)), "double")
112%!assert (class (betarnd (single (1), 2)), "single")
113%!assert (class (betarnd (single ([1 1]), 2)), "single")
114%!assert (class (betarnd (1, single (2))), "single")
115%!assert (class (betarnd (1, single ([2 2]))), "single")
116
117## Test input validation
118%!error betarnd ()
119%!error betarnd (1)
120%!error betarnd (ones (3), ones (2))
121%!error betarnd (ones (2), ones (3))
122%!error betarnd (i, 2)
123%!error betarnd (2, i)
124%!error betarnd (1,2, -1)
125%!error betarnd (1,2, ones (2))
126%!error binornd (1,2, [2 -1 2])
127%!error betarnd (1,2, 1, ones (2))
128%!error betarnd (1,2, 1, -1)
129%!error betarnd (ones (2,2), 2, 3)
130%!error betarnd (ones (2,2), 2, [3, 2])
131%!error betarnd (ones (2,2), 2, 2, 3)
132