indxs2( log(u2(indxs2)) < (min_ab-1)*log(w(indxs2)/t) ) ] ) = 1;                              clear('u1');                              clear('u2');                              if a == min_ab                                   out( indxs(l) ) = w(l);                              else                                   out( indxs(l) ) = 1 - w(l);                              end                              indxs = indxs( ~l );                         end                    else                         % Atkinson's Algorithm                         if min_ab == 1                              t = 0.5;                              r = 0.5;                         else                              t = 1/(1+sqrt(max_ab*(1-max_ab)/(min_ab*(1-min_ab))));                              r = max_ab*t / (max_ab*t + min_ab*(1-t));                         end                         u1 = rand( sampleSize );                         out = zeros( sampleSize );                         w  = zeros( sampleSize );                         l1 = u1 < r;                         w(l1) = t*(u1(l1)/r).^(1/min_ab);                         l2 = u1 >= r;                         w(l2) = 1 - (1-t)*((1-u1(l2))/(1-r)).^(1/max_ab);                         if a == min_ab                              out = w;                         else                              out = 1 - w;                         end                         u2 = rand( sampleSize );                         indxs1 = find(l1);                         indxs2 = find(l2);                         indxs = [ indxs1( log(u2(l1)) >= (max_ab -1)*log((1-w(l1))/(1-t)) ), ...                              indxs2( log(u2(l2)) >= (min_ab -1) * log(w(l2)/t) ) ];                         clear('u2');                         while ~isempty( indxs )                              indxsSize = size( indxs );                              u1 = rand( indxsSize );                              w  = zeros( indxsSize );                              l1 = u1 < r;                              w(l1) = t*(u1(l1)/r).^(1/min_ab);                              l2 = u1 >= r;                              w(l2) = 1 - (1-t)*((1-u1(l2))/(1-r)).^(1/max_ab);                              u2 = rand( indxsSize );                              indxs1 = find(l1);                              indxs2 = find(l2);                              clear('u1');                              l = logical( zeros( indxsSize ) );                              l( [ indxs1(log(u2(l1)) < (max_ab -1)*log((1-w(l1))/(1-t))), ...                                   indxs2(log(u2(l2))< (min_ab -1) * log(w(l2)/t)) ] ) = 1;                              if a == min_ab                                   out(indxs(l)) = w(l);                              else                                   out(indxs(l)) = 1 - w(l);                              end                              indxs = indxs( ~l );                         end                    end                    out = m + (n-m) * out;                    reshape( out, sampleSizeIn );               case {'bessel'}                    % START bessel HELP                    %  THE BESSEL DISTRIBUTION                    %                    %  Bessel distribution arises in the theory of stochastic processes.                    %  Bessel(nu,a) is a discrete distribution on the non-negative integers with                    %  parameters nu > -1 and a > 0.                    %                    % pdf(y) = (a/2).^(2*y+nu) ./ (besseli(nu,a).*factorial(y).*gamma(y+nu+1));                    %                    % PARAMETERS:                    %   nu > -1, a > 0                    % SUPPORT:                    %   y = 0, 1, 2, 3, ...                    % CLASS:                    %   Discrete distributions                    %                    % USAGE:                    %   randraw('bessel', [nu, a], sampleSize) - generate sampleSize number                    %         of variates from the Bessel distribution with parameters                    %         'nu' and 'a'                    %   randraw('bessel') - help for the Bessel distribution;                    % EXAMPLES:                    %  1.   y = randraw('bessel', [2 0.9], [1 1e5]);                    %  2.   y = randraw('bessel', [0.6 3.2], 1, 1e5);                    %  3.   y = randraw('bessel', [-0.2 8.1], 1e5 );                    %  4.   y = randraw('bessel', [4 5.3], [1e5 1] );                    %  5.   randraw('bessel');                    % END bessel HELP                    % Method:                    %                    % We implemented Condensed Table-Lookup method suggested in                    %    George Marsaglia, "Fast Generation Of Discrete Random Variables,"                    %    Journal of Statistical Software, July 2004, Volume 11, Issue 4                    %                    % Reference:                    % L. Devroye, "Simulating Bessel random variables,"                    %  Statistics and Probability Letters, vol. 57, pp. 249-257, 2002.                    %                    checkParamsNum(funcName, 'Bessel', 'bessel', distribParams, [2]);                    nu = distribParams(1);                    a = distribParams(2);                    validateParam(funcName, 'Bessel', 'bessel', '[nu, a]', 'nu', nu, {'> -1'});                    validateParam(funcName, 'Bessel', 'bessel', '[nu, a]', 'a', a, {'> 0'});                    % mu = 0.5*a*besseli(nu+1,a)/besseli(nu,a);                    % chi2 = mu + 0.25*a^2*besseli(nu+1,a)/besseli(nu,a)*...                    %     (besseli(nu+2,a)/besseli(nu+1,a)-besseli(nu+1,a)/besseli(nu,a));                    besseliNuA =  besseli(nu, a);                    proceed = 1;                    if ~isfinite( besseliNuA )                         warnStr{1} = [upper(funcName), ' - Bessel Variates Generation: '];                         warnStr{2} = ['besseli(', num2str(nu), ', ' num2str(a), ') returns Inf.'];                         warnStr{3} = ['Unable to proceed, return zeros ...'];                         warning('%s\n  %s\n  %s',warnStr{1},warnStr{2},warnStr{3});                         %warning([upper(funcName), ' - Bessel Variates Generation: besseli(', num2str(nu), ', ' num2str(a), ') returns Inf. Unable to proceed, return zeros ...']);                         out = zeros( sampleSize );                         proceed = 0;                    end                    if besseliNuA == 0                         warnStr{1} = [upper(funcName), ' - Bessel Variates Generation: '];                         warnStr{2} = ['besseli(', num2str(nu), ', ' num2str(a), ') returns 0.'];                         warnStr{3} = ['Unable to proceed, return zeros ...'];                         warning('%s\n  %s\n  %s',warnStr{1},warnStr{2},warnStr{3});                         %warning([upper(funcName), '- Bessel Variates Generation: besseli(', num2str(nu), ', ' num2str(a), ') returns 0. Unable to proceed, return zeros ...']);                         out = zeros( sampleSize );                         proceed = 0;                    end                    if proceed                         p0 = exp( nu*log(a/2) - gammaln(nu+1) ) / besseliNuA;                         if p0 >= 5e-10                              t = p0;                              aa = (a/2)^2;                              nu1 = nu+1;                              i = 1;                              while t*2147483648 > 1                                   t = t * aa/((i)*(i+nu));                                   i = i + 1;                              end                              sizeP = i-1;                              offset = 0;                              P = round( 2^30*p0*cumprod([1, aa./((1:sizeP-1).*((1:sizeP-1)+nu))] ) );                         else % if p0 >= 5e-10                              m = floor(0.5*(sqrt(a^2+nu^2)-nu));                              pm = exp( (2*m+nu)*log(a/2) - log(besseliNuA) - ...                                   gammaln(m+1) - gammaln(m+nu+1) );                              aa = (a/2)^2;                              t = pm;                              i = m + 1;                              while t * 2147483648 > 1                                   t = t * aa/((i)*(i+nu));                                   i = i + 1;                              end                              last = i-2;                              t = pm;                              j = -1;                              for i = m-1:-1:0                                   t = t * (i+1)*(i+1+nu)/aa;                                   if t*2147483648 < 1                                        j=i;                                        break;                                   end                              end                              offset = j+1;                              sizeP = last-offset+1;                              P = zeros(1, sizeP);                              P(m-offset+1:last-offset+1) = ...                                   round( 2^30*pm*cumprod([1, aa./(((m+1):last).*(((m+1):last)+nu))] ) );                              P(m-offset:-1:1) = ...                                   round( 2^30*pm*cumprod((m:-1:(offset+1)).*((m:-1:(offset+1))+nu)/aa) );                         end % if p0 >= 5e-10, else ...                         out = randFrom5Tbls( P, offset, sampleSize);                    end % if proceed               case {'binom', 'binomial'}                    % START binom HELP START binomial HELP                    % THE BINOMIAL DISTRIBUTION                    %                    % pdf(y) = nchoosek(n,y)*p^y*(1-p)^(n-y) = ...                    %          exp( gammaln(n+1) - gammaln(n-y+1) - gammaln(y+1) + ...                    %               y*log(p) + (n-y)*log(1-p) );  0<p<1, n>1                    %                    %  Mean = n*p;                    %  Variance = n*p*(1-p);                    %  Mode = floor( (n+1)*p );                    %                    % PARAMETERS:                    %   p  - probability of success in a single trial; (0<p<1)                    %   n  - total number of trials; (n= 1, 2, 3, 4, ...)                    %                    % SUPPORT:                    %   y - number of success,  y = 0, 1, 2, 3 ...                    %                    % CLASS:                    %   Discrete distributions                    %                    % NOTES:                    %   Constructive definition:                    %    We consider a random experiment with n independent trials; in each trial                    %    a certain random event A can occur (the urn model with replacement is                    %    a special case of such an experiment). Let                    %      p  = probability of A in a single trial;                    %      n  = total number of trials;                    %      y  = number of successes (= number of trials where A occurs).                    %                    % USAGE:                    %   randraw('binom', [n, p], sampleSize) - generate sampleSize number                    %         of variates from the Binomial distribution with total number of trials                    %         'n' and probability of success in a single trial 'p'                    %   randraw('binom') - help for the Binomial distribution;                    %                    % EXAMPLES:                    %  1.   y = randraw('binom', [10 0.9], [1 1e5]);                    %  2.   y = randraw('binom', [100 0.15], 1, 1e5);                    %  3.   y = randraw('binom', [5 0.5], 1e5 );                    %  4.   y = randraw('binom', [1000 0.02], [1e5 1] );                    %  5.   randraw('binom');                    % END binom HELP END binomial HELP                    % Method:                    %                    % We implemented Condensed Table-Lookup method suggested in                    %    George Marsaglia, "Fast Generation Of Discrete Random Variables,"                    %    Journal of Statistical Software, July 2004, Volume 11, Issue 4                    checkParamsNum(funcName, 'Binomial', 'binomial', distribParams, [2]);                    n = distribParams(1);                    p = distribParams(2);                    validateParam(funcName, 'Binomial', 'binomial', '[n, p]', 'n', n, {'> 0','==integer'});                    validateParam(funcName, 'Binomial', 'binomial', '[n, p]', 'p', p, {'> 0','< 1'});                    if n*p > 1e3 & n > 1e3                         out = round( n*p + sqrt(n*p*(1-p))*randn( sampleSize ) );                    elseif p<1e-4 & n*p > 1 & n*p < 100                         out = feval(funcName,'poisson',n*p, sampleSize);                    else                         % if n large and p near 1, generate j=Binom(n,1-p), return n-j                         switchFlag = 0;                         if n>100 & p>0.99                              p = 1-p;                              switchFlag = 1;                         end                         mode = floor( (n+1)*p );                         q = 1 - p;                         h = p/q;                         pmode = exp( gammaln(n+1) - gammaln(n-mode+1) - gammaln(mode+1) + ...                              mode*log(p) + (n-mode)*log(1-p) );                         if( 1-pmode < 5e-10 )                              % "Fast Generation of Discrete Random Variables", p.3  citation:                              % "For ... discrete variates with an infinite number of probabilities, we select only                              % those for which, for a sample of size 2^31(10^9.33), the expected number of occurences exceeds                              % 0.5. The other probabilities are assumed zero. For those unusual situations where occurrences                              % with probability less than 5