Code covered by the BSD License

# Generation of Random Variates

### James Huntley (view profile)

generates random variates from over 870 univariate distributions

bin2ordk_pdf(n, k, p, nn)
```% bin2ordk_pdf.m - evaluates a Type II Binomial Order K Probability denisity.
%   See "Univariate Discrete Distributions", Johnson, Kemp, and Kotz,
%
%   Note:  Be sure that 'nn' > 2*'k' for normalization.
%
%   For nn = 40 (& p=0.5), 'distribs20.m' requires ~45 minutes to generate 1000 points!!!
%
%   Created by  J. Huntley,  10/17/07.
%
%

function [pdf] = bin2ordk_pdf(n, k, p, nn)

%persistent q qdp logqdp pn

%Initializations.
%if(isempty(q))
q = 1 - p;
qdp = q / p;
logqdp = log(qdp);
pn = p^nn;
%end

% Evaluate PDF.
if(n == 0)
pdf = pn;

elseif(n > 0)

% Main summantion loop over 'i'.
sumi = 0;
for ii = 1:nn
i = ii - 1;
n0 = nn - i;

% Fetch pre-stored partitions of 'n'. Frequencies returned in array, "d".
pname = ['partition' num2str(n0)];
spc = size(d,2);
d = double(d);

% Select rows of 'd' with only min(nn,spc) columns populated and store in array, 'xs'.
xlim = min(nn,spc);
xs = d(:,1:xlim);

% Find subsets of rows of 'xs' such that the sum of the last n-k elements equals
% 'n' (call this subset 'xs1') and 'n-1' (call this subset 'xs2').
sc = size(xs,2);
sr = size(xs,1);
ir1 = 0;
ir2 = 0;
xs1 = [];
xs2 = [];
for jr = 1:sr
if(sum(xs(jr,k+1:end)) == n-1)
ir1 = ir1 + 1;
xs1(ir1,:) = xs(jr,:);
end
if(sum(xs(jr,k+1:end)) == n-2)
ir2 = ir2 + 1;
xs2(ir2,:) = xs(jr,:);
end
end

% Calculate PDF.
% Sum over solutions up to order 'nn' for Diophantine Equation.
sum1 = 0;
if(ii <= k)
sr = size(xs1,1);
sc = size(xs1,2);
elseif(ii > k)
sr = size(xs2,1);
sc = size(xs2,2);
end
% Summation loop over common (intersecting) solutions.
for jr = 1:sr
if(ii <= k)
x(1:sc) = xs1(jr,1:sc);
elseif(ii > k)
x(1:sc) = xs2(jr,1:sc);
end
%fx(1:sc) = factorial(x(1:sc));
fx(1:sc) = gammaln(x(1:sc)+1);
sumx = sum(x);
%pfx = prod(fx);
pfx = sum(fx);
%sum1 = sum1 + factorial(sumx)*(qdp)^sumx / pfx;
sum1 = sum1 + exp(gammaln(sumx+1)+ sumx*logqdp - pfx);
end % summation loop over common solutions.
sumi = sumi + sum1;
end % main summation loop over i.
pdf = pn * sumi;

end % conditional for n > 0.

return

```