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

negmulti_pdf(n, sp1, k, nnmax, p)
% negmulti_pdf.m - evaluates a Negative Multinomial Distribution.
%   See "Univariate Discrete Distributions", Johnson, Kemp & Kotz, 
%   J. Wiley, 2005, p.552.
%
%   Created by: J. Huntley
%
%   NOTE: Text appears to give incomplete description of this pdf.  In
%   order to normalize this PDF, must also sum over all nn ('r' in text).
%   As a result, parameter range is quite limited for reasonable running
%   times.  In general, (sp1-1)*nnmax+1 ~< 35, even though partitions are 
%   already pre-stored up to a value of 44. It also helps to have the p's
%   in descending order.
%   Also note that Pr(X=0) = p(1)^k.
%
%   Requires: 'partitionXX.mat', pre-stored partitions created by
%   'create_partition.m'. Note that max(XX) is currently 44.
%

function[sumpdf] = negmulti_pdf(n, sp1, k, nnmax, p)

% Initializations.

nn = 1:nnmax;
pdf0 = p(1)^k;
sumpdf = pdf0;
       
if(n > eps)
    sumpdf = 0;
    % 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(sp1-1,size(d,2));
    xs = int8(d(:,1:xlim));

    for jnn = 1:size(nn,2)
        zeroterm = exp(gammaln(k+nn(jnn)) + k*log(p(1)) - gammaln(k));
        % 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(jnn))
                ir1 = ir1 + 1;
                xs1(ir1,:) = int8(xs(jr,:));
            end
            if(sc >= sd  & sum(d(jr,1:sc)) == nn(jnn))
                ir1 = ir1 + 1;
                xs1(ir1,:) = int8(xs(jr,:));
            end            
        end

        % Calculate PDF.
        % Sum over the subset of solutions to Diaphantine Eq. such that
        % the sum of all the 'r's = '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 = 1;
            prodpdr = 0;
            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;
                prodpdr = prodpdr + pdr;
            end
            %sum1 = sum1 + prodpdr
            sum1 = sum1 + exp(prodpdr);
        end % summation loop over common solutions.
        pdf = sum1 * zeroterm;
        sumpdf = sumpdf + pdf;

    end % Main Loop over 'nn'.

end % if n(jn) > eps. 

return

Contact us