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

waring_pdf(n, c, a)
% waring_pdf.m - evaluates a Waring Probability Density.
%   See "Dataplot Reference Manual, WARPDF", AUX-361, NIST, 3/26/97..
%
%  Created by Jim Huntley,  11/06/03
%

function [pdf] = waring_pdf(n, c, a)

%persistent gamlna logc gamlncp1 logcma

%if(isempty(logc))
    gamlna = gammaln(a);
    logc = log(c);
    gamlncp1 = gammaln(c+1);
    logcma = log(c-a);
%end

if(c+n+1 > 160)					% Use Stirling's Approx. to reduce ratio of Gamma functions.
    p1 = 1 + 1/(12*(a+n)) + 1/(288*(a+n)^2) - 139/(51840*(a+n)^3) - 571/(2488320*(a+n)^4);
    p2 = 1 + 1/(12*(c+1+n)) + 1/(288*(c+1+n)^2) - 139/(51840*(c+1+n)^3) - 571/(2488320*(c+1+n)^4);
    pdf = exp(c-a+1) * ((a+n)/(c+1+n))^(n) * (a+n)^(a-0.5) * p1 / ((c+n+1)^(c+0.5) * p2);
elseif(c+n+1 <= 160)            
    %pdf = ((c-a) * gamma(a+n) * gamma(c+1)) / (c * gamma(c+n+1) * gamma(a));
    pdf = exp(logcma + gammaln(a+n) + gamlncp1 - (logc + gammaln(c+n+1) + gamlna));
end

return

Contact us