Code covered by the BSD License

# Generation of Random Variates

### 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

```