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

poisordk_pdf(n, k, lambda)
% poisordk_pdf.m - evaluates a Poisson Order K Probability denisity.
%   See "Univariate Discrete Distributions", Johnson, Kemp, and Kotz,
%   Wiley, p.460, 2005.  See also "http://mathworld.wolfram.com/Partition.html".
%
%   Created by  J. Huntley,  09/05/07.
%
%   Loads files 'partition1-N.dat'
%

function [pdf] = poisordk_pdf(n, k, lambda)

%Initializations.
coef = exp(-k*lambda);
ll = log(lambda);

if(n == 0)
    pdf = coef;
elseif(n > 0)
    
    % Fetch pre-stored partitions of 'n'. Frequencies returned in array, "d".
    pname = ['partition' num2str(n)];
    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);
        pfx = sum(fx);
        %sum1 = sum1 + factorial(sumx)*(q/p)^sumx / pfx;
        sum1 = sum1 + exp(sumx*ll - pfx);
    end
    pdf = coef * sum1;
end  % n > 0

return


    

Contact us