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

bin1ordk_pdf(n, k, p, nn)
% bin1ordk_pdf.m - evaluates a Type I 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".
%
%   Created by  J. Huntley,  10/12/07.
%
%   Loads 'partition1-N.dat'
%

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

%persistent q qdp logqdp pn nnn

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

sumi = 0;
for ii = 1:k
    i = ii - 1;
    n0 = nnn-i-n*k;
    if(n0 == 0)
        pdf = pn;
    elseif(n0 > 0) 
        % Fetch pre-stored partitions of 'n'. Frequencies returned in array, "d".
        pname = ['partition' num2str(n0)];
        load(pname, 'd');
        spc = size(d,2);
        spr = size(d,1);
        d = double(d);

        % Select rows of 'd' with only min(k,spc) columns populated and store in array, 'xs'.
        indx = 0;
        for jr = 1:spr
            if(size(find(d(jr,k+1:spc)),2)== 0)
                indx = indx + 1;
                dd(indx,:) = d(jr,:);
            end        
        end
        xlim = min(k,spc);
        xs = dd(1:indx,1:xlim);

        % Calculate PDF.
        % Sum over solutions up to order 'k' for Diophantine Equation. 
        sum1 = 0;
        sx1 = size(xs,1);
        for jr = 1:sx1
            for jc = 1:xlim
                x(jc) = xs(jr,jc);
                %fx(jc) = factorial(x(jc));
                fx(jc) = gammaln(x(jc)+1);
            end
            sumx = sum(x);
            %pfx = prod(fx)*n;
            pfx = sum(fx)+ gammaln(n+1);
            %sum1 = sum1 + factorial(sumx+n)*qdp^sumx / pfx;
            sum1 = sum1 + exp(gammaln(sumx+n+1)+ sumx*logqdp - pfx);
        end
        sumi = sumi + sum1;    
        pdf = pn * sumi;
        clear L pp d dd xs
    end  % (n-k) > 0
end

return


    

Contact us