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

ipbinom_pdf(n, ppi, rho, N)
% ipbinom_pdf.m - evaluates an Inflated Parameter Binomial Probability denisity.
%   See "Inflated-Parameter Family of Generalized Power Series Distributions
%   and their Application in Analysis of Overdispersed Insurance Data", 
%   N. Kolev et al. ARCH 00V211.PDF
%
%  Created by Jim Huntley,  10/23/09
%
%   Loads files 'partition1-N.dat'
%

function [pdf] = ipbinom_pdf(n, ppi, rho, N)

%persistent glnNp1 logp1mr log1mp logrho pn 

%Initializations.
%if(isempty(pn))
    glnNp1 = gammaln(N+1);
    pn = (1-ppi)^N;
    logp1mr = log(ppi*(1-rho));
    log1mp = log(1-ppi);
    logrho = log(rho);
%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);
    d = double(d);

    % Calculate PDF.
    % Sum over all partition frequencies to get solution to Diophantine Equation. 
    sum1 = 0;
    d1 = size(d,1);
    for jr = 1:d1
         for jc = 1:spc
            x(jc) = d(jr,jc);
            dx(jc) = (jc-1) * x(jc);
            %fx(jc) = factorial(x(jc));
            fx(jc) = gammaln(x(jc)+1);
        end
        sumx = sum(x);
        N0 = N - sumx;
        if(N0 >= 0)
            glnN0p1 = gammaln(N0+1);
            sumdx = sum(dx); 
            %pfx = prod(fx) * factorial(N0);            
            pfx = sum(fx) + glnN0p1;
            %sum1 = sum1 + factorial(N) * (1-ppi)^N0 * (ppi*(1-rho))^sumx * rho^sumdx / pfx;
            sum1 = sum1 + exp(glnNp1 + N0*log1mp + sumx*logp1mr + sumdx*logrho - pfx);
        end
    end
    pdf = sum1; 

end % conditional for n > 0.

return


    

Contact us