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

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