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_cdf(n, c, a)
% waring_cdf.m - evaluates a Cumulative Waring Distribution.
%   See "Dataplot Reference Manual, WARCDF", AUX-361, NIST, 3/26/97..
%
%  Created by Jim Huntley,  11/06/03
%
%

function [cdf] = waring_cdf(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

for js = 1:size(n,2)
    nlm = n(js);
    p = zeros(1,nlm);    
	for jn = 1:nlm					% Sum terms from 1 to limit of n.
        if(c+jn+1 > 160)			% Use Stirling's Approx. to reduce ratio of Gamma functions.
            p1 = 1 + 1/(12*(a+jn)) + 1/(288*(a+jn)^2) - 139/(51840*(a+jn)^3) - 571/(2488320*(a+jn)^4);
            p2 = 1 + 1/(12*(c+1+jn)) + 1/(288*(c+1+jn)^2) - 139/(51840*(c+1+jn)^3) - 571/(2488320*(c+1+jn)^4);
            p(jn) = exp(c-a+1) * ((a+jn)/(c+1+jn))^(jn) * (a+jn)^(a-0.5) * p1 / ((c+jn+1)^(c+0.5) * p2);
            %p(jn) = p(jn) * (c-a) * gamma(c+1) / (c * gamma(a));
            p(jn) = p(jn) * (c-a) * exp(gamlncp1 - (logc + gamlna));
        elseif(c+jn+1 <= 160)            
            %p(jn) = ((c-a) * gamma(a+jn) * gamma(c+1)) / (c * gamma(c+jn+1) * gamma(a));
            p(jn) = exp(logcma + gammaln(a+jn) + gamlncp1 - (logc + gammaln(c+jn+1) + gamlna));            
        end
    end
    
    cdf(js) = sum(p);
end
cdf = cdf + (c-a) / c;			% Add n=0 term.

return

Contact us