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

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

function[pdf] = factmult_pdf(n, sp1, nn, N)

%persistent ss NN denom

%if(isempty(ss))
    ss = sp1 - 1;
    NN = sum(N(1:ss+1));
    denom = 1 / binomial_coef(NN,nn);
%end

pdf = 0;
if(n > eps)
    % Fetch pre-stored partitions of 'n'. Frequencies returned in array, "d".
    pname = ['partition' num2str(n)];
    load(pname, 'd');            
    spc = size(d,2)
    spr = size(d,1);
    xlim = min(ss,spc);
    d = double(d)

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

    % 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,:) = 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);
        prodbcNr = binomial_coef(N(1),r0(jr));
        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 * denom;

end % if n(jn) > eps. 

return

Contact us