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

ipnegbin_pdf(n, ppi, rho, r)
% ipnegbin_pdf.m - evaluates an Inflated Parameter Negative 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/27/09
%
%   Loads files 'partition1-N.dat'
%

function [pdf] = ipnegbin_pdf(n, ppi, rho, r)

%persistent log1mpr logrho pn

%Initializations.
%if(isempty(pn))
    pn = ppi^r;
    log1mpr = log((1-ppi)*(1-rho));
    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);
        %sumxprm1 = sumx + r - 1;
        sumxpr = sumx + r;
        sumdx = sum(dx); 
        %pfx = prod(fx)*factorial(r-1);
        pfx = sum(fx)+ gammaln(r);
        %sum1 = sum1 + factorial(sumxprm1)*((1-ppi)*(1-rho))^sumx*rho^sumdx / pfx;
        sum1 = sum1 + exp(gammaln(sumxpr) + sumx*log1mpr + sumdx*logrho - pfx);
    end
    pdf = pn * sum1; 

end % conditional for n > 0.

return


    

Contact us