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