Daniel@0: function r = betarnd(a,b,m,n);
Daniel@0: %BETARND Random matrices from beta distribution.
Daniel@0: %   R = BETARND(A,B) returns a matrix of random numbers chosen   
Daniel@0: %   from the beta distribution with parameters A and B.
Daniel@0: %   The size of R is the common size of A and B if both are matrices.
Daniel@0: %   If either parameter is a scalar, the size of R is the size of the other
Daniel@0: %   parameter. Alternatively, R = BETARND(A,B,M,N) returns an M by N matrix. 
Daniel@0: 
Daniel@0: %   Reference:
Daniel@0: %      [1]  L. Devroye, "Non-Uniform Random Variate Generation", 
Daniel@0: %      Springer-Verlag, 1986
Daniel@0: 
Daniel@0: %   Copyright (c) 1993-98 by The MathWorks, Inc.
Daniel@0: %   $Revision: 1.1.1.1 $  $Date: 2005/04/26 02:29:18 $
Daniel@0: 
Daniel@0: if nargin < 2, 
Daniel@0:     error('Requires at least two input arguments'); 
Daniel@0: end 
Daniel@0: 
Daniel@0: if nargin == 2
Daniel@0:     [errorcode rows columns] = rndcheck(2,2,a,b);
Daniel@0: end
Daniel@0: 
Daniel@0: if nargin == 3
Daniel@0:     [errorcode rows columns] = rndcheck(3,2,a,b,m);
Daniel@0: end
Daniel@0: 
Daniel@0: if nargin == 4
Daniel@0:     [errorcode rows columns] = rndcheck(4,2,a,b,m,n);
Daniel@0: end
Daniel@0: 
Daniel@0: if errorcode > 0
Daniel@0:     error('Size information is inconsistent.');
Daniel@0: end
Daniel@0: 
Daniel@0: r = zeros(rows,columns);
Daniel@0: 
Daniel@0: % Use Theorem 4.1, case A (Devroye, page 430) to derive beta
Daniel@0: %   random numbers as a ratio of gamma random numbers.
Daniel@0: if prod(size(a)) == 1
Daniel@0:     a1 = a(ones(rows,1),ones(columns,1));
Daniel@0:     g1 = gamrnd(a1,1);
Daniel@0: else
Daniel@0:     g1 = gamrnd(a,1);
Daniel@0: end
Daniel@0: if prod(size(b)) == 1
Daniel@0:     b1 = b(ones(rows,1),ones(columns,1));
Daniel@0:     g2 = gamrnd(b1,1);
Daniel@0: else
Daniel@0:     g2 = gamrnd(b,1);
Daniel@0: end
Daniel@0: r = g1 ./ (g1 + g2);
Daniel@0: 
Daniel@0: % Return NaN if b is not positive.
Daniel@0: if any(any(b <= 0));
Daniel@0:     if prod(size(b) == 1)
Daniel@0:         tmp = NaN;
Daniel@0:         r = tmp(ones(rows,columns));
Daniel@0:     else
Daniel@0:         k = find(b <= 0);
Daniel@0:         tmp = NaN;
Daniel@0:         r(k) = tmp(ones(size(k)));
Daniel@0:     end
Daniel@0: end
Daniel@0: 
Daniel@0: % Return NaN if a is not positive.
Daniel@0: if any(any(a <= 0));
Daniel@0:     if prod(size(a) == 1)
Daniel@0:         tmp = NaN;
Daniel@0:         r = tmp(ones(rows,columns));
Daniel@0:     else
Daniel@0:         k = find(a <= 0);
Daniel@0:         tmp = NaN;
Daniel@0:         r(k) = tmp(ones(size(k)));
Daniel@0:     end
Daniel@0: end