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

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