Code covered by the BSD License

# Generation of Random Variates

### James Huntley (view profile)

generates random variates from over 870 univariate distributions

stutgenwar_pdf(n, k, a, c, m)
% stutgenwar_pdf.m - evaluates a Stuttering Generalized Waring Probability denisity.
%   See "The Stuttering Generalized Waring Distribution", J.Panaretos & E.
%   Xekalaki, Statistics and Prob. Lett. V4, No.6, 1986, 313.
%   http://mpra.ub.uni-muenchen.de/6250/1/MPRA_paper_6250.pdf.
%
%  Created by Jim Huntley,  12/29/09
%
%

function [pdf] = stutgenwar_pdf(n, k, a, c, m)

%persistent M cM apcM gamlnacM gamlna pn

%Initializations.
%if(isempty(pn))
M = sum(m);
cM = exp(gammaln(c+M) - gammaln(c));
apcM = exp(gammaln(a+c+M) - gammaln(a+c));
gamlnacM = gammaln(a+c+M);
gamlna = gammaln(a);
pn = cM / apcM;
%end

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

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

% Reduce array "dd" to only those rows whose first k columns contain solutions
% to the Diophantine equations.
ir = 0;
minkspc = min(k,spc);
for jr = 1:spr
if(spc <= k)
ir = ir + 1;
dd(ir,1:minkspc) = d(jr,1:minkspc);
elseif(spc > k)
sumgtk = sum(d(jr,k+1:spc));
if(sumgtk <= 0)
ir = ir + 1;
dd(ir,1:k) = d(jr,1:k);
end
end
end

% Calculate PDF.
% Sum over all partition frequencies to get solution to Diophantine Equation.
sum1 = 0;
dd1 = size(dd,1);
for jr = 1:dd1
for jc = 1:minkspc
x(jc) = dd(jr,jc);
fx(jc) = gammaln(x(jc)+1);
end
sumx = sum(x);
asumx = exp(gammaln(a+sumx) - gamlna);
apcpMsumx = exp(gammaln(a+c+M+sumx) - gamlnacM);
prod1 = 1;
for jp = 1:min(k,spc)
msumx = (gammaln(m(jp)+x(jp)) - gammaln(m(jp)));
prod1 = prod1 * exp(msumx - fx(jp));
end
sum1 = sum1 + asumx * prod1 / apcpMsumx;
end
pdf = pn * sum1;

end % conditional for n > 0.

return