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

...
% disades_rnd.m - generates 'nsamples' of a Discrete Ades random variable.
%   from "Univariate Discrete Distributions", Johnson, Kemp, & Kotz,
%   Wiley, p.495, 2005.
%
%   Created by: J. Huntley,  12/08/06
%

function [xhold,fhold,ghold,jhold,xsort,fsort,gsort,chold,isint] = ...
          disades_rnd(userin,p_value,nsamples)

% Initializations.
hbar = waitbar(0,'Please wait...');
a = p_value(1);
b = 1 / p_value(2);
c = p_value(3);
xmax = userin.xmax_value;
xmin = userin.xmin_value;
userin.distrib = 'gamma'; 
userin.ftype = 4;
userin.p_value = [a 0 b];


% Sample Gamma Distribution.
[xhold,fhold,ghold,jhold,xsort,fsort,gsort,chold,isint] =  gen_distrib2(userin,0);

% Generate samples from Ades Distribution.
for jx = 1:nsamples
    xxx = 0;
    if(xhold(jx) > 1)
        xxx = log(xhold(jx))^c;
    end
    xhold(jx) = xxx;
%    xhold0(jx) = jx;
end
%pdf0 = xhold;
%figure(40)
%plot(xhold0,pdf0);
%drawnow;

% Generate Discrete Ades PDF.
Pr(1:xmax+1) = 0;
%nhold = 0:xmax;
for js = 1:nsamples
    if(xhold(js) < 0.5)
        Pr(1) = Pr(1) + 1;
    elseif(xhold(js) >= 0.5)
        indx = fix(xhold(js)+0.5) + 1;
        Pr(indx) = Pr(indx) + 1;
    end
end
pdf = Pr;
%figure(41)
%plot(nhold,pdf);
%drawnow;

% Generate Discrete Ades CDF.
cdf(1:xmax+1) = 0;
for jn = 1:xmax+1;
    for jp = 1:jn
        cdf(jn) = cdf(jn) + pdf(jp);
    end
end
cdf = cdf ./ nsamples;
%figure(42)
%plot(nhold,cdf);
%drawnow;

% Sample CDF for Discrete Ades Distribution.
for js = 1:nsamples 
    waitbar(js / nsamples);
    urd1 = rand(1);   			% uniform random deviate (0 : 1)
    %Find closest cdf to urd1.    
    xm = xmin;
    for jy = 1:xmax+1
        if(urd1 < cdf(jy))
            xm = jy - 1;
            break;
        end
    end
    xhold(js) = xm;
end % Loop over Samples for Discrete Ades Distribution.

 close(hbar); 

return

Contact us