sgvm_gamma_variate

SGVM_GAMMA_VARIATE generate gamma distributed samples
   RVS = SGVM_GAMMA_VARIATE(N, A) (default scale is 1)
   RVS = SGVM_GAMMA_VARIATE(N, A, B)

   Generate N samples from the gamma distribution with shape A and scale
   B. The gamma distribution is:
           f(x| a, b) = 1/(gamma(k)*b^k) * x^(k-1) * exp(-x/b)
   This sampling follows the algorithm described here:
   https://en.wikipedia.org/wiki/Gamma_distribution#Generating_gamma-distributed_random_variables

0001 function rvs = sgvm_gamma_variate(n, a, b)
0002 %SGVM_GAMMA_VARIATE generate gamma distributed samples
0003 %   RVS = SGVM_GAMMA_VARIATE(N, A) (default scale is 1)
0004 %   RVS = SGVM_GAMMA_VARIATE(N, A, B)
0005 %
0006 %   Generate N samples from the gamma distribution with shape A and scale
0007 %   B. The gamma distribution is:
0008 %           f(x| a, b) = 1/(gamma(k)*b^k) * x^(k-1) * exp(-x/b)
0009 %   This sampling follows the algorithm described here:
0010 %   https://en.wikipedia.org/wiki/Gamma_distribution#Generating_gamma-distributed_random_variables
0011 
0012 %   SynGrid
0013 %   Copyright (c) 2018, Power Systems Engineering Research Center (PSERC)
0014 %   by Eran Schweitzer, Arizona State University
0015 %
0016 %   This file is part of SynGrid.
0017 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0018 
0019 if nargin < 3
0020     b = 1;
0021 end
0022 
0023 if ~(a > 0)
0024     error('sgvm_gamma_variate: shape parameter ''a'' must be greater than 0.')
0025 elseif~(b > 0)
0026     error('sgvm_gamma_variate: scale parameter ''b'' must be greater than 0.')
0027 end
0028 
0029 x = floor(a);
0030 delta = a - x;
0031 
0032 gamma_n1 = -sum(log(rand(n, x)), 2);
0033 
0034 if delta ~= 0
0035     flag = false(n,1);
0036     
0037     zeta = zeros(n,1);
0038     eta  = zeros(n,1);
0039     
0040     while ~all(flag)
0041         idx = find(~flag);
0042         n = length(idx);
0043         U = rand(n,1);
0044         V = rand(n,1);
0045         W = rand(n,1);
0046         
0047         test = U <= exp(1)/(exp(1) + delta);
0048         
0049         zeta(idx(test)) = V(test).^(1/delta);
0050         eta(idx(test))  = W(test).*zeta(idx(test)).^(delta - 1);
0051         
0052         zeta(idx(~test)) = 1 - log(V(~test));
0053         eta(idx(~test))  = W(~test).*exp(-zeta(idx(~test)));
0054         
0055         flag = eta <= zeta.^(delta - 1) .* exp(-zeta);
0056     end
0057 else
0058     zeta = 0;
0059 end
0060 
0061 rvs = b*(zeta + gamma_n1);