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

nbinordk_pdf(n, k, p, r, b)
% nbinordk_pdf.m - evaluates a Negative Binomial Order K Probability denisity.
%   See "Univariate Discrete Distributions", Johnson, Kemp, and Kotz,
%   Wiley, p.458, 2005.  See also "http://mathworld.wolfram.com/Partition.html".
%
%   For nn = 40 (& p=0.5), 'distribs20.m' requires ~45 minutes to generate 1000 points!!!
%
%   Created by  J. Huntley,  7/19/07.
%
%   Loads files 'partition1-N.dat'
%

function [pdf] = nbinordk_pdf(n, k, p, r, b)

%persistent q logqdp

%Initializations.
%if(isempty(q))
    q = 1 - p;
    logqdp = log(q/p);
%end

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

return


    

Contact us