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

stutgenwar_pdf(n, k, a, c, m)
% stutgenwar_pdf.m - evaluates a Stuttering Generalized Waring Probability denisity.
%   See "The Stuttering Generalized Waring Distribution", J.Panaretos & E.
%   Xekalaki, Statistics and Prob. Lett. V4, No.6, 1986, 313.
%   http://mpra.ub.uni-muenchen.de/6250/1/MPRA_paper_6250.pdf.
%
%  Created by Jim Huntley,  12/29/09
%
%   Loads files 'partition1-N.dat'
%

function [pdf] = stutgenwar_pdf(n, k, a, c, m)

%persistent M cM apcM gamlnacM gamlna pn

%Initializations.
%if(isempty(pn))
    M = sum(m);
    cM = exp(gammaln(c+M) - gammaln(c));
    apcM = exp(gammaln(a+c+M) - gammaln(a+c));
    gamlnacM = gammaln(a+c+M);
    gamlna = gammaln(a);
    pn = cM / apcM;
%end

% Evaluate PDF.
if(n == 0)
    pdf = pn;

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);
    
    % Reduce array "dd" to only those rows whose first k columns contain solutions
    % to the Diophantine equations.
    ir = 0;
    minkspc = min(k,spc);
    for jr = 1:spr
        if(spc <= k)
            ir = ir + 1;
            dd(ir,1:minkspc) = d(jr,1:minkspc);
        elseif(spc > k)
            sumgtk = sum(d(jr,k+1:spc));
            if(sumgtk <= 0)
                ir = ir + 1;
                dd(ir,1:k) = d(jr,1:k);
            end
        end
    end             

    % Calculate PDF.
    % Sum over all partition frequencies to get solution to Diophantine Equation. 
        sum1 = 0;
        dd1 = size(dd,1);
        for jr = 1:dd1
            for jc = 1:minkspc
                x(jc) = dd(jr,jc);
                fx(jc) = gammaln(x(jc)+1);
            end
            sumx = sum(x);
            asumx = exp(gammaln(a+sumx) - gamlna);
            apcpMsumx = exp(gammaln(a+c+M+sumx) - gamlnacM);
            prod1 = 1;
            for jp = 1:min(k,spc)
                msumx = (gammaln(m(jp)+x(jp)) - gammaln(m(jp)));
                prod1 = prod1 * exp(msumx - fx(jp));
            end
            sum1 = sum1 + asumx * prod1 / apcpMsumx;
        end
        pdf = pn * sum1;

end % conditional for n > 0.

return


    

Contact us