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

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