Code covered by the BSD License

# Generation of Random Variates

### James Huntley (view profile)

generates random variates from over 870 univariate distributions

gamma_pdf(x,a,m,b)
```function y = gamma_pdf(x,a,m,b)
%GAMMA_PDF Gamma probability density function.
%   Y = GAMMA_PDF(X,A,B) returns the gamma probability density function
%   with parameters A and B, at the values in X.
%
%   The size of Y is the common size of the input arguments. A scalar input
%   functions as a constant matrix of the same size as the other inputs.
%
%   Some references refer to the gamma distribution with a single
%   parameter. This corresponds to the default of B = 1.

%   References:
%      [1]  L. Devroye, "Non-Uniform Random Variate Generation",
%      Springer-Verlag, 1986, pages 401-402.

%   Copyright 1993-2000 The MathWorks, Inc.
%   \$Revision: 2.9 \$  \$Date: 2000/05/26 18:52:53 \$

if nargin < 4,
b = 1;
end

if nargin < 3,
error('Requires at least two input arguments');
end

x = x - m;
[errorcode x a b] = distchck(3,x,a,b);

if errorcode > 0
error('Requires non-scalar arguments to match in size.');
end

% Initialize Y to zero.
y = zeros(size(x));

%   Return NaN if the arguments are outside their respective limits.
y(a <= 0 | b <= 0) = NaN;

k=find(x > 0 & ~(a <= 0 | b <= 0));
if any(k)
y(k) = (a(k) - 1) .* log(x(k)) - (x(k) ./ b(k)) - gammaln(a(k)) - a(k) .* log(b(k));
y(k) = exp(y(k));
end
y(x == 0 & a < 1) = Inf;
k2 = find(x == 0 & a == 1);
if any(k2)
y(k2) = (1./b(k2));
end

return
```