wolffd@0: function r = gamrnd(a,b,m,n);
wolffd@0: %GAMRND Random matrices from gamma distribution.
wolffd@0: %   R = GAMRND(A,B) returns a matrix of random numbers chosen   
wolffd@0: %   from the gamma 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 = GAMRND(A,B,M,N) returns an M by N matrix. 
wolffd@0: %
wolffd@0: %   Some references refer to the gamma distribution
wolffd@0: %   with a single parameter. This corresponds to GAMRND
wolffd@0: %   with B = 1. (See Devroye, pages 401-402.)
wolffd@0: 
wolffd@0: %   GAMRND uses a rejection or an inversion method depending on the
wolffd@0: %   value of A. 
wolffd@0: 
wolffd@0: %   References:
wolffd@0: %      [1]  L. Devroye, "Non-Uniform Random Variate Generation", 
wolffd@0: %      Springer-Verlag, 1986
wolffd@0: 
wolffd@0: %   B.A. Jones 2-1-93
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: 
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: % Initialize r to zero.
wolffd@0: lth = rows*columns;
wolffd@0: r = zeros(lth,1);
wolffd@0: a = a(:); b = b(:);
wolffd@0: 
wolffd@0: scalara = (length(a) == 1);
wolffd@0: if scalara 
wolffd@0:    a = a*ones(lth,1);
wolffd@0: end
wolffd@0: 
wolffd@0: scalarb = (length(b) == 1);
wolffd@0: if scalarb 
wolffd@0:    b = b*ones(lth,1);
wolffd@0: end
wolffd@0: 
wolffd@0: % If a == 1, then gamma is exponential. (Devroye, page 405).
wolffd@0: k = find(a == 1);
wolffd@0: if any(k)
wolffd@0:    r(k) = -b(k) .* log(rand(size(k)));
wolffd@0: end 
wolffd@0: 
wolffd@0: 
wolffd@0: k = find(a < 1 & a > 0);
wolffd@0: % (Devroye, page 418 Johnk's generator)
wolffd@0: if any(k)
wolffd@0:   c = zeros(lth,1);
wolffd@0:   d = zeros(lth,1);
wolffd@0:   c(k) = 1 ./ a(k);
wolffd@0:   d(k) = 1 ./ (1 - a(k));
wolffd@0:   accept = k;
wolffd@0:   while ~isempty(accept)
wolffd@0:     u = rand(size(accept));
wolffd@0:     v = rand(size(accept));
wolffd@0:     x = u .^ c(accept);
wolffd@0:     y = v .^ d(accept);
wolffd@0:     k1 = find((x + y) <= 1); 
wolffd@0:     if ~isempty(k1)
wolffd@0:       e = -log(rand(size(k1))); 
wolffd@0:       r(accept(k1)) = e .* x(k1) ./ (x(k1) + y(k1));
wolffd@0:       accept(k1) = [];
wolffd@0:     end
wolffd@0:   end
wolffd@0:   r(k) = r(k) .* b(k);
wolffd@0: end
wolffd@0: 
wolffd@0: % Use a rejection method for a > 1.
wolffd@0: k = find(a > 1);
wolffd@0: % (Devroye, page 410 Best's algorithm)
wolffd@0: bb = zeros(size(a));
wolffd@0: c  = bb;
wolffd@0: if any(k)
wolffd@0:   bb(k) = a(k) - 1;
wolffd@0:   c(k) = 3 * a(k) - 3/4;
wolffd@0:   accept = k; 
wolffd@0:   count = 1;
wolffd@0:   while ~isempty(accept)
wolffd@0:     m = length(accept);
wolffd@0:     u = rand(m,1);
wolffd@0:     v = rand(m,1);
wolffd@0:     w = u .* (1 - u);
wolffd@0:     y = sqrt(c(accept) ./ w) .* (u - 0.5);
wolffd@0:     x = bb(accept) + y;
wolffd@0:     k1 = find(x >= 0);
wolffd@0:     if ~isempty(k1)
wolffd@0:       z = 64 * (w .^ 3) .* (v .^ 2);
wolffd@0:       k2 = (z(k1) <= (1 - 2 * (y(k1) .^2) ./ x(k1)));
wolffd@0:       k3 = k1(find(k2));
wolffd@0:       r(accept(k3)) = x(k3); 
wolffd@0:       k4 = k1(find(~k2));
wolffd@0:       k5 = k4(find(log(z(k4)) <= (2*(bb(accept(k4)).*log(x(k4)./bb(accept(k4)))-y(k4)))));
wolffd@0:       r(accept(k5)) = x(k5);
wolffd@0:       omit = [k3; k5];
wolffd@0:       accept(omit) = [];
wolffd@0:     end
wolffd@0:   end
wolffd@0:   r(k) = r(k) .* b(k);
wolffd@0: end
wolffd@0: 
wolffd@0: % Return NaN if a or b is not positive.
wolffd@0: r(b <= 0 | a <= 0) = NaN;
wolffd@0: 
wolffd@0: r = reshape(r,rows,columns);