Code covered by the BSD License  

Highlights from
Generation of Random Variates

image thumbnail

Generation of Random Variates

by

 

generates random variates from over 870 univariate distributions

multinom_pdf(n, sp1, nn, p)
% multinom_pdf.m - evaluates a Multinomial Distribution.
%   See "Univariate Discrete Distributions", Johnson, Kemp & Kotz, 
%   J. Wiley, 2005, p.523.
%
%   Created by: J. Huntley, 11/20/08
%
%   Loads 'partition0-N.dat'
%

function[pdf] = multinom_pdf(n, sp1, nn, p)

%persistent ss gamnnp1

%if(isempty(ss))
    ss = sp1 - 1;
    gamnnp1 = gammaln(nn+1);
%end

pdf = p(1)^nn;
if(n > eps)
    % Fetch pre-stored partitions of 'n'. Frequencies returned in array, "d".
    pname = ['partition' num2str(n)];
    load(pname, 'd');

    % Select rows of 'd' with only min(ss,spc) columns populated and store in array, 'xs'. 
    % AND that sum to 'nn'.
    xlim = min(ss,size(d,2));
    xs = int8(d(:,1:xlim));

    % Find subset of rows of 'xs' of length 'xlim' that still satisfy the Diophantine Equation and
    sc = size(xs,2);
    sd = size(d,2);
    sr = size(xs,1);
    ir1 = 0;
    xs1 = [];
    for jr = 1:sr           
        if(d(jr,sc+1:end) == 0 & sum(d(jr,1:sc)) <= nn)
            ir1 = ir1 + 1;
            xs1(ir1,:) = int8(xs(jr,:));
        end
        if(sc >= sd  & sum(d(jr,1:sc)) <= nn)
            ir1 = ir1 + 1;
            xs1(ir1,:) = xs(jr,:);
        end             
    end
    
    % Estimate 'r0'.
    r0 = [];
    for jr = 1:ir1
        sumd = sum(xs1(jr,1:xlim));
        r0(jr) = max(0,nn-sumd);
    end

    %xs1n = xs1
    %r0

    % Calculate PDF.
    % Sum over the subset of solutions to Diaphantine Eq. such that
    % the sum of all the 'r's (including 'r0') = 'nn'. 
    sum1 = 0;
    sr = size(xs1,1);
    sc = size(xs1,2);                         % Summation loop over common (intersecting) solutions.
    for jr = 1:sr
        r(1:sc) = xs1(jr,1:sc);
        %prodpdr = p(1)^r0(jr) / gamma(r0(jr)+1);
        prodpdr = (r0(jr)*log(p(1)) - gammaln(r0(jr)+1));
        for jc = 1:sc
            %pdr = p(jc+1)*r(jc) / gamma(r(jc)+1);
            pdr = (r(jc)*log(p(jc+1)) - gammaln(r(jc)+1));
            prodpdr = prodpdr + pdr;
        end
        sum1 = sum1 + exp(prodpdr);
        end % summation loop over common solutions.
    pdf = exp(log(sum1) + gamnnp1);

end % if n(jn) > eps. 

return

Contact us