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

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