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

bin3ordk_pdf(n, k, p, nn)
% bin3ordk_pdf.m - evaluates a Type III Binary order K Probability denisity.
%   See "Univariate Discrete Distributions", Johnson, Kemp, and Kotz,
%   J. Wiley, p.526, 2005.  
%
%   Note sign error in text.
%
%   Created by  J. Huntley,  10/18/07.
%
%   Loads 'partition1-N.dat'
%

function [pdf] = bin3ordk_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;
        xs = [];

       % 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;            
        xs1 = [];           
        for jr = 1:sr  
            sumj = 0;
            for jc = k+1:sc
                sumj = sumj + (jc-k)*xs(jr,jc);
            end              
            %if(sumj == n-max(0,ii-k)-1)
            if(sumj == n-1)
                ir1 = ir1 + 1;
                xs1(ir1,:) = xs(jr,:);
            end
        end    

        % Calculate PDF.
        % Sum over solutions up to order 'nn' for Diophantine Equation. 
        sum1 = 0;
        sr = size(xs1,1);
        sc = size(xs1,2);             
        % Summation loop over common (intersecting) solutions.
        for jr = 1:sr                
            x(1:sc) = xs1(jr,1:sc);  
            %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