Code covered by the BSD License  

Highlights from
Generation of Random Variates

image thumbnail

Generation of Random Variates

by

James Huntley (view profile)

 

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