Code covered by the BSD License  

Highlights from
Generation of Random Variates

image thumbnail

Generation of Random Variates

by

James Huntley (view profile)

 

generates random variates from over 870 univariate distributions

gammaln(x)
function res = gammaln(x)
%GAMMALN Logarithm of gamma function.
%   Y = GAMMALN(X) computes the natural logarithm of the gamma
%   function for each element of X.  GAMMALN is defined as
%
%       LOG(GAMMA(X))
%
%   and is obtained without computing GAMMA(X).  Since the gamma
%   function can range over very large or very small values, its
%   logarithm is sometimes more useful.
%
%   See also GAMMA, GAMMAINC.

%   C. Moler, 2-1-91.
%   Copyright 1984-2001 The MathWorks, Inc. 
%   $Revision: 5.11 $  $Date: 2001/04/15 12:01:43 $

%   This is based on a FORTRAN program by W. J. Cody,
%   Argonne National Laboratory, NETLIB/SPECFUN, June 16, 1988.
%
% References:
%
%  1) W. J. Cody and K. E. Hillstrom, 'Chebyshev Approximations for
%     the Natural Logarithm of the Gamma Function,' Math. Comp. 21,
%     1967, pp. 198-203.
%
%  2) K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May,
%     1969.
% 
%  3) Hart, Et. Al., Computer Approximations, Wiley and sons, New
%     York, 1968.
%
%#mex

     d1 = -5.772156649015328605195174e-1;
     p1 = [4.945235359296727046734888e0; 2.018112620856775083915565e2; 
           2.290838373831346393026739e3; 1.131967205903380828685045e4; 
           2.855724635671635335736389e4; 3.848496228443793359990269e4; 
           2.637748787624195437963534e4; 7.225813979700288197698961e3];
     q1 = [6.748212550303777196073036e1; 1.113332393857199323513008e3; 
           7.738757056935398733233834e3; 2.763987074403340708898585e4; 
           5.499310206226157329794414e4; 6.161122180066002127833352e4; 
           3.635127591501940507276287e4; 8.785536302431013170870835e3];
     d2 = 4.227843350984671393993777e-1;
     p2 = [4.974607845568932035012064e0; 5.424138599891070494101986e2; 
           1.550693864978364947665077e4; 1.847932904445632425417223e5; 
           1.088204769468828767498470e6; 3.338152967987029735917223e6; 
           5.106661678927352456275255e6; 3.074109054850539556250927e6];
     q2 = [1.830328399370592604055942e2; 7.765049321445005871323047e3; 
           1.331903827966074194402448e5; 1.136705821321969608938755e6; 
           5.267964117437946917577538e6; 1.346701454311101692290052e7; 
           1.782736530353274213975932e7; 9.533095591844353613395747e6];
     d4 = 1.791759469228055000094023e0;
     p4 = [1.474502166059939948905062e4; 2.426813369486704502836312e6; 
           1.214755574045093227939592e8; 2.663432449630976949898078e9; 
           2.940378956634553899906876e10; 1.702665737765398868392998e11; 
           4.926125793377430887588120e11; 5.606251856223951465078242e11];
     q4 = [2.690530175870899333379843e3; 6.393885654300092398984238e5; 
           4.135599930241388052042842e7; 1.120872109616147941376570e9; 
           1.488613728678813811542398e10; 1.016803586272438228077304e11; 
           3.417476345507377132798597e11; 4.463158187419713286462081e11];
     c = [-1.910444077728e-03; 8.4171387781295e-04; 
          -5.952379913043012e-04; 7.93650793500350248e-04; 
          -2.777777777777681622553e-03; 8.333333333333333331554247e-02; 
           5.7083835261e-03];

   if ~isempty(x) & (any((imag(x) ~= 0) | (x < 0)))
      error('Input arguments must be real and nonnegative.')
   end
   res = x;
%
%  0 <= x <= eps
%
   k = find(x <= eps);
   if ~isempty(k)
      res(k) = -log(x(k));
   end
%
%  eps < x <= 0.5
%
   k = find((x > eps) & (x <= 0.5));
   if ~isempty(k)
      y = x(k);
      xden = ones(size(k));
      xnum = 0;
      for i = 1:8
         xnum = xnum .* y + p1(i);
         xden = xden .* y + q1(i);
      end
      res(k) = -log(y) + (y .* (d1 + y .* (xnum ./ xden)));
   end
%
%  0.5 < x <= 0.6796875
%
   k = find((x > 0.5) & (x <= 0.6796875));
   if ~isempty(k)
      xm1 = (x(k) - 0.5) - 0.5;
      xden = ones(size(k));
      xnum = 0;
      for i = 1:8
         xnum = xnum .* xm1 + p2(i);
         xden = xden .* xm1 + q2(i);
      end
      res(k) = -log(x(k)) + xm1 .* (d2 + xm1 .* (xnum ./ xden));
   end
%
%  0.6796875 < x <= 1.5
%
   k = find((x > 0.6796875) & (x <= 1.5));
   if ~isempty(k)
      xm1 = (x(k) - 0.5) - 0.5;
      xden = ones(size(k));
      xnum = 0;
      for i = 1:8
         xnum = xnum .* xm1 + p1(i);
         xden = xden .* xm1 + q1(i);
      end
      res(k) = xm1 .* (d1 + xm1 .* (xnum ./ xden));
   end
%
%  1.5 < x <= 4
%
   k = find((x > 1.5) & (x <= 4));
   if ~isempty(k)
      xm2 = x(k) - 2;
      xden = ones(size(k));
      xnum = 0;
      for i = 1:8
         xnum = xnum .* xm2 + p2(i);
         xden = xden .* xm2 + q2(i);
      end
      res(k) = xm2 .* (d2 + xm2 .* (xnum ./ xden));
   end
%
%  4 < x <= 12
%
   k = find((x > 4) & (x <= 12));
   if ~isempty(k)
      xm4 = x(k) - 4;
      xden = -ones(size(k));
      xnum = 0;
      for i = 1:8
         xnum = xnum .* xm4 + p4(i);
         xden = xden .* xm4 + q4(i);
      end
      res(k) = d4 + xm4 .* (xnum ./ xden);
   end
%
%  x > 12
%
   k = find(x > 12);
   if ~isempty(k)
      y = x(k);
      r = c(7)*ones(size(k));
      ysq = y .* y;
      for i = 1:6
         r = r ./ ysq + c(i);
      end
      r = r ./ y;
      corr = log(y);
      spi = 0.9189385332046727417803297;
      res(k) = r + spi - 0.5*corr + y .* (corr-1);
   end

Contact us