Code covered by the BSD License

# Generation of Random Variates

### James Huntley (view profile)

generates random variates from over 870 univariate distributions

bin3ordk_pdf(n, k, p, nn)
```% bin3ordk_pdf.m - evaluates a Type III Binary order K Probability denisity.
%   See "Univariate Discrete Distributions", Johnson, Kemp, and Kotz,
%   J. Wiley, p.526, 2005.
%
%   Note sign error in text.
%
%   Created by  J. Huntley,  10/18/07.
%
%

function [pdf] = bin3ordk_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;
xs = [];

% 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;
xs1 = [];
for jr = 1:sr
sumj = 0;
for jc = k+1:sc
sumj = sumj + (jc-k)*xs(jr,jc);
end
%if(sumj == n-max(0,ii-k)-1)
if(sumj == n-1)
ir1 = ir1 + 1;
xs1(ir1,:) = xs(jr,:);
end
end

% Calculate PDF.
% Sum over solutions up to order 'nn' for Diophantine Equation.
sum1 = 0;
sr = size(xs1,1);
sc = size(xs1,2);
% Summation loop over common (intersecting) solutions.
for jr = 1:sr
x(1:sc) = xs1(jr,1:sc);
%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

```