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

dblpois_pdf(n, mu, psi, nmax)
% dblpois_pdf.m - evaluates a Double Poisson Probability Density.
%   See "Univariate Discrete Distributions", Johnson, Kemp & Kotz,
%   J. Wiley, p.489, 2005.
%
%   NOTE:  Efron Normalization is good only for mu AND psi >~ 3 !!! 
%   If mu AND psi <= 3, then program will fit normalizing constant to a 2-D
%   CDF surface with < 1% normalization error. If mu OR psi <= 3, then program uses 
%   Efron approximation and normalization error depends strongly on psi/mu.  
%   If this ratio is high, then so is normalization error.  In this case,
%   the user should turn on the auto-normalization in the Excel database!!!
%
%  Created by Jim Huntley,  12/20/06
%

function [pdf] = dblpois_pdf(n, mu, psi, nmax)

%persistent coef

%if(isempty(coef))
    mulim = 3.0;
    psilim = 3.0;
    dmu = 0.1;
    dpsi = 0.1;
    if(mu > mulim+eps || psi > psilim+eps)
        c = 1 / (1 + (1-psi) * (1+1/(mu*psi)) / (12*mu*psi));     % Efron approximation.
    elseif(mu <= mulim+eps && psi <= psilim+eps)
        mumin = dmu;
        mumax = mulim;        
        mu0 = mumin:dmu:mumax;
        psimin = dpsi;
        psimax = psilim;        
        psi0 = psimin:dpsi:psimax;
        randfrac = 0.0;                 % If program warns about singularity in surface fit,
                                        % then make randfrac > 0 until warning goes away.
        pdf0(1:nmax) = 0;
        
        % Make CDF(mu0,psi0) surface.
        for jm = 1:size(mu0,2)
            for jp = 1:size(psi0,2)
                c = 1;
                coef = log(c * sqrt(psi0(jp)));
                for jn = 1:nmax
                    pdf0(jn) = exp(coef + (-mu0(jm)*psi0(jp)-jn) + jn*log(jn) + ...
                              (jn*psi0(jp))*log(mu0(jm)*exp(1)/jn)- gammaln(jn+1));
                end
                ipdf0(jm,jp) = sum(pdf0);
                clear pdf0
            end
        end

        % Plot CDF(mu0,psi0) surface.
        %[X,Y] = meshgrid(mu0,psi0,ipdf0');
        %figure(41)
        %surf(X,Y,ipdf0');
        %drawnow;

        % Add noise to stabilize fit.
        ipdf0 = ipdf0 + (rand(size(ipdf0,2),size(ipdf0,1))-0.5)*randfrac;
        %figure(42)
        %surf(X,Y,ipdf0');

        % Weighted least-squares fit to polynomial surface.
        orderx = 6;
        ordery = 6;
        w = ones(size(ipdf0',1),size(ipdf0',2));
        [p, S] = polyfit2w(X,Y,ipdf0',w',orderx,ordery);

        % Computer residuals (errors).
        %for je = 1:size(ipdf0,2)
        %    for jr = 1:size(ipdf0,1)      
        %        bherr(jr,je) = polyval2d(p,jr*dmu,je*dpsi) - ipdf0(jr,je);            
        %    end
        %end

        % Plot Residuals.
        %[XX,YY] = meshgrid(mu0,psi0,ipdf0');
        %figure(43)
        %surf(XX,YY,bherr);
        %drawnow;
        
        % Compute Normalization constant.
        c = 1 / polyval2d(p,mu,psi);
    end
    coef = log(c * sqrt(psi));
%end

pdf = exp(coef + (-mu*psi-n) + n*log(n) + (n*psi)*log(mu*exp(1)/n)...
      - gammaln(n+1));

return


Contact us