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

negfacmulti_pdf(n, sp1, k, N)
% negfacmulti_pdf.m - evaluates a Negative Factorial Multinomial Distribution.
%   See "Univariate Discrete Distributions", Johnson, Kemp & Kotz, 
%   J. Wiley, 2005, p.522.
%
%   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] = negfacmulti_pdf(n, sp1, k, N)

%persistent NN nnmax nn nmax zeroterm P0

%if(isempty(NN))
    NN = sum(N(1:sp1));
    nnmax = NN-k;
    nn = 0:nnmax;
    sumn = 0;
    for js = 2:sp1
        sumn = sumn + (js-1)*N(js);
    end
    nmax = sumn;
    zeroterm = (N(1) - k + 1) * binomial_coef(N(1),k-1);
    P0 = (N(1)-k+1)*N(1) / ((NN-k+1)*(k-1)*binomial_coef(NN,k-1));
%end

sumpdf = 0;
if(n <= nmax)
    pdf = P0;   
    if(n > eps)
        sumpdf = P0;
        % 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)
            % Find subset of rows of 'xs' of length 'xlim' that still satisfy
            % the Diophantine Equation
            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);
                prodbcNr = 1;
                    for jc = 1:sc
                        bcNr = binomial_coef(N(jc+1),r(jc));
                        prodbcNr = prodbcNr * bcNr;
                    end                
                    sum1 = sum1 + prodbcNr;
                end % summation loop over common solutions.
                pdf = sum1 * zeroterm / (binomial_coef(NN,k+nn(jnn)-1)*(NN-k-nn(jnn)+1));
                sumpdf = sumpdf + pdf;

        end % Main Loop over 'nn'.

    end % if n > eps. 
    
end % if(n <= nmax)

return

Contact us