Code covered by the BSD License

# Generation of Random Variates

### James Huntley (view profile)

generates random variates from over 870 univariate distributions

bin1ordk_pdf(n, k, p, nn)
```% bin1ordk_pdf.m - evaluates a Type I Binomial Order K Probability denisity.
%   See "Univariate Discrete Distributions", Johnson, Kemp, and Kotz,
%
%   Created by  J. Huntley,  10/12/07.
%
%

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

%persistent q qdp logqdp pn nnn

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

sumi = 0;
for ii = 1:k
i = ii - 1;
n0 = nnn-i-n*k;
if(n0 == 0)
pdf = pn;
elseif(n0 > 0)
% Fetch pre-stored partitions of 'n'. Frequencies returned in array, "d".
pname = ['partition' num2str(n0)];
spc = size(d,2);
spr = size(d,1);
d = double(d);

% Select rows of 'd' with only min(k,spc) columns populated and store in array, 'xs'.
indx = 0;
for jr = 1:spr
if(size(find(d(jr,k+1:spc)),2)== 0)
indx = indx + 1;
dd(indx,:) = d(jr,:);
end
end
xlim = min(k,spc);
xs = dd(1:indx,1:xlim);

% Calculate PDF.
% Sum over solutions up to order 'k' for Diophantine Equation.
sum1 = 0;
sx1 = size(xs,1);
for jr = 1:sx1
for jc = 1:xlim
x(jc) = xs(jr,jc);
%fx(jc) = factorial(x(jc));
fx(jc) = gammaln(x(jc)+1);
end
sumx = sum(x);
%pfx = prod(fx)*n;
pfx = sum(fx)+ gammaln(n+1);
%sum1 = sum1 + factorial(sumx+n)*qdp^sumx / pfx;
sum1 = sum1 + exp(gammaln(sumx+n+1)+ sumx*logqdp - pfx);
end
sumi = sumi + sum1;
pdf = pn * sumi;
clear L pp d dd xs
end  % (n-k) > 0
end

return

```