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

bin2ordk_pdf(n, k, p, nn)
% bin2ordk_pdf.m - evaluates a Type II Binomial Order K Probability denisity.
%   See "Univariate Discrete Distributions", Johnson, Kemp, and Kotz,
%   Wiley, p.462, 2005.  See also "http://mathworld.wolfram.com/Partition.html".
%
%   Note:  Be sure that 'nn' > 2*'k' for normalization.
%
%   For nn = 40 (& p=0.5), 'distribs20.m' requires ~45 minutes to generate 1000 points!!!
%
%   Created by  J. Huntley,  10/17/07.
%
%   Loads 'partition1-N.dat'
%

function [pdf] = bin2ordk_pdf(n, k, p, nn)

%persistent q qdp logqdp pn 

%Initializations.
%if(isempty(q))
    q = 1 - p;
    qdp = q / p;
    logqdp = log(qdp);
    pn = p^nn;
%end

% Evaluate PDF.
if(n == 0)
    pdf = pn;

elseif(n > 0)

    % Main summantion loop over 'i'.
    sumi = 0;
    for ii = 1:nn
        i = ii - 1;
        n0 = nn - i;
       
        % Fetch pre-stored partitions of 'n'. Frequencies returned in array, "d".
        pname = ['partition' num2str(n0)];
        load(pname, 'd');
        spc = size(d,2);
        d = double(d);

        % Select rows of 'd' with only min(nn,spc) columns populated and store in array, 'xs'.
        xlim = min(nn,spc);
        xs = d(:,1:xlim);

        % Find subsets of rows of 'xs' such that the sum of the last n-k elements equals
        % 'n' (call this subset 'xs1') and 'n-1' (call this subset 'xs2').
        sc = size(xs,2);
        sr = size(xs,1);
        ir1 = 0;
        ir2 = 0;
        xs1 = [];
        xs2 = [];
        for jr = 1:sr
            if(sum(xs(jr,k+1:end)) == n-1)
                ir1 = ir1 + 1;
                xs1(ir1,:) = xs(jr,:);
            end
            if(sum(xs(jr,k+1:end)) == n-2)
                ir2 = ir2 + 1;
                xs2(ir2,:) = xs(jr,:);
            end
        end

        % Calculate PDF.
        % Sum over solutions up to order 'nn' for Diophantine Equation. 
        sum1 = 0;
        if(ii <= k)
            sr = size(xs1,1);
            sc = size(xs1,2);                
        elseif(ii > k)
            sr = size(xs2,1);
            sc = size(xs2,2);
        end
        % Summation loop over common (intersecting) solutions.
        for jr = 1:sr                
            if(ii <= k)
                x(1:sc) = xs1(jr,1:sc);
            elseif(ii > k)
                x(1:sc) = xs2(jr,1:sc);
            end
            %fx(1:sc) = factorial(x(1:sc));
            fx(1:sc) = gammaln(x(1:sc)+1);                
            sumx = sum(x);
            %pfx = prod(fx);
            pfx = sum(fx);
            %sum1 = sum1 + factorial(sumx)*(qdp)^sumx / pfx;
            sum1 = sum1 + exp(gammaln(sumx+1)+ sumx*logqdp - pfx);
        end % summation loop over common solutions.
        sumi = sumi + sum1;      
    end % main summation loop over i.
    pdf = pn * sumi; 

end % conditional for n > 0.

return


    

Contact us