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

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