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

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