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

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