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

stutpois_pdf(n, k, lambda)
% stutpois_pdf.m - evaluates a Stuttering Poisson 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,  12/07/06.
%
%   Loads files 'partition1-N.dat'
%

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

%Initializations.
coef = exp(-k*lambda);
loglam = 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 = pro(fx);
        %pfx = exp(sum(fx));
        pfx = sum(fx);
        %sum1 = sum1 + lambda^sumx / pfx;
        sum1 = sum1 + exp(sumx*loglam - pfx);
    end
    pdf = coef * sum1;
end  % n > 0

return


    

Contact us