Code covered by the BSD License

# Generation of Random Variates

### James Huntley (view profile)

generates random variates from over 870 univariate distributions

gamma(x,a)
```function res = gamma(x,a)
%GAMMA  Gamma function.
%   Y = GAMMA(X) evaluates the gamma function for each element of X.
%   X must be real.  The gamma function is defined as:
%
%      gamma(x) = integral from 0 to inf of t^(x-1) exp(-t) dt.
%
%   The gamma function interpolates the factorial function.  For
%   integer n, gamma(n+1) = n! (n factorial) = prod(1:n).
%

%   C. B. Moler, 5-7-91, 11-4-92.
%   Ref: Abramowitz & Stegun, Handbook of Mathematical Functions, sec. 6.1.
%   Copyright 1984-2000 The MathWorks, Inc.
%   \$Revision: 5.14 \$  \$Date: 2000/07/28 18:59:48 \$

%   This is based on a FORTRAN program by W. J. Cody,
%   Argonne National Laboratory, NETLIB/SPECFUN, October 12, 1989.
%
% References: "An Overview of Software Development for Special
%              Functions", W. J. Cody, Lecture Notes in Mathematics,
%              506, Numerical Analysis Dundee, 1975, G. A. Watson
%              (ed.), Springer Verlag, Berlin, 1976.
%
%              Computer Approximations, Hart, Et. Al., Wiley and
%              sons, New York, 1968.
%
%#mex

if nargin == 2
% For backward compatibility with MATLAB 3.5
warning(sprintf([ ...
'This usage of gamma(a,x) is obsolete and will be eliminated\n' ...
res = gammainc(a,x);
return
end

p = [-1.71618513886549492533811e+0; 2.47656508055759199108314e+1;
-3.79804256470945635097577e+2; 6.29331155312818442661052e+2;
8.66966202790413211295064e+2; -3.14512729688483675254357e+4;
-3.61444134186911729807069e+4; 6.64561438202405440627855e+4];
q = [-3.08402300119738975254353e+1; 3.15350626979604161529144e+2;
-1.01515636749021914166146e+3; -3.10777167157231109440444e+3;
2.25381184209801510330112e+4; 4.75584627752788110767815e+3;
-1.34659959864969306392456e+5; -1.15132259675553483497211e+5];
c = [-1.910444077728e-03; 8.4171387781295e-04;
-5.952379913043012e-04; 7.93650793500350248e-04;
-2.777777777777681622553e-03; 8.333333333333333331554247e-02;
5.7083835261e-03];

if ~isreal(x)
error('Input argument must be real.')
end
res = zeros(size(x));
xn = zeros(size(x));
%
%  Catch negative x.
%
kneg = find(x <= 0);
if ~isempty(kneg)
y = -x(kneg);
y1 = fix(y);
res(kneg) = y - y1;
fact = -pi ./ sin(pi*res(kneg)) .* (1 - 2*rem(y1,2));
x(kneg) = y + 1;
end
%
%  x is now positive.
%  Map x in interval [0,1] to [1,2]
%
k1 = find(x < 1);
x1 = x(k1);
x(k1) = x1 + 1;
%
%  Map x in interval [1,12] to [1,2]
%
k = find(x < 12);
xn(k) = fix(x(k)) - 1;
x(k) = x(k) - xn(k);
%
%  Evaluate approximation for 1 < x < 2
%
if ~isempty(k)
z = x(k) - 1;
xnum = 0*z;
xden = xnum + 1;
for i = 1:8
xnum = (xnum + p(i)) .* z;
xden = xden .* z + q(i);
end
res(k) = xnum ./ xden + 1;
end
%
%  Adjust result for case  0.0 < x < 1.0
%
res(k1) = res(k1) ./ x1;
%
%  Adjust result for case  2.0 < x < 12.0
%
for j = 1:max(xn(:))
k = find(xn);
res(k) = res(k) .* x(k);
x(k) = x(k) + 1;
xn(k) = xn(k) - 1;
end
%
%  Evaluate approximation for x >= 12
%
k = find(x >= 12);
if ~isempty(k)
y = x(k);
ysq = y .* y;
sum = c(7);
for i = 1:6
sum = sum ./ ysq + c(i);
end
spi = 0.9189385332046727417803297;
sum = sum ./ y - y + spi;
sum = sum + (y-0.5).*log(y);
res(k) = exp(sum);
end
%