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

distribs31(userin,iseed)
% distribs31.m - program to generate random numbers from any of a wide variety of 
%                distributions via numerical inversion of the CDF or other techniques.
%
% Created by:  Jim Huntley,  12/06/2011.
%
% Calls: 'gen_distrib.m', 'cdf4d_m', 'pdf4d_m', + '*_cdf' (See "distrib_full" for 
%         list), '*_pdf' (See "distrib_full" for list) + *_rnd. 
%
% Called by:  GUI_02.m.
%
% Inputs:  User-Supplied Distribution Selection and Parameters -
%           as specified in spreadsheet, 'noise_distrib22.xls', and 
%           taken from GUI input.
%
% Outputs: Plots of "time series", CDF, PDF, and sample distribution (histogram).
%          Diagnostic output in Matlab workspace.   ASCII time series output can be
%          saved via GUI to 'user_named.dat'
%

function [xhold,fhold,ghold,jhold,chold,xsort,gsort,fsort,isint,xmax_value,xmin_value] = distribs31(userin,iseed)

global hObjectE11 hObjectE12 hObjectE16 hObjectE20
global c_value

%Get Initial Input Data from GUI.

% Unpack Input Data Structure.
d_indx = userin.d_indx;
distrib = userin.distrib;
p_value = userin.p_value;
xmax_value = userin.xmax_value;
xmin_value = userin.xmin_value;
ftype = userin.ftype;
fnorm = userin.fnorm;
nsamples = userin.nsamples;
name = userin.name;
fxmin = userin.fxmin;
fxmax = userin.fxmax;
fwarn = userin.fwarn;

ifdone = 0;

% Compute additional Input Parameters Required for Non-Central Chi-Squared Distribution.
% Note that these Additional Parameters are not Saved - only SQRT(SUM(Squares)) is Kept.
if(strcmp(distrib,'nc_chisq') == 1)							
    %fprintf('\n The %s Distribution requires additional value(s) for %d input parameter(s).\n', name, p_value(1));
    out_string = char(['INPUT  ' num2str(p_value(1)) ' ADDITIONAL PARAMETER(S) BELOW. Then hit "Enter" to Continue']);
    set(hObjectE16,'String',out_string);
    set(hObjectE20,'Visible','on');
    set(hObjectE20,'BackgroundColor',[1.0 0.969 0.922]);
    pause
    waitfor(hObjectE20,'BackgroundColor','white');
    sum = 0;
    for jp = 1:p_value(1)
            pv = c_value(jp);
            sum = sum + pv^2;                
    end
    p_value(3) = sqrt(sum);
end 
c_value = [];

% Overwrite User's Choices for Support Limits because Max and Min are directly Specified
% by the Input Parameters.
if(fxmin > 0)
    xmin_value = sign(fxmin)*p_value(abs(fxmin))+1e-7*sign(fxmin);
end
if(fxmax > 0)
    xmax_value = sign(fxmax)*p_value(abs(fxmax))-1e-7;
end

% Overwrite User's Choices for Support Limits because Max and Min are
% Specified by the Following Functions of the Input Parameters.
if(strcmp(distrib,'absorp') == 1)
    xmax_value = min(p_value(2),p_value(3));
end
if(strcmp(distrib,'amoroso') == 1)
    if(p_value(4) > 0)
        xmin_value = p_value(1) + 1e-4;
        xmax_value = p_value(5);
    elseif(p_value(4) < 0)
        xmax_value = p_value(1) - 1e-4;
        xmin_value = -p_value(5);
    end
end
if(strcmp(distrib,'beta1') == 1)
    xmin_value = p_value(1) + 1e-4; 
    xmax_value = p_value(2) - 1e-4; 
end
if(strcmp(distrib,'bin1ordk') == 1)
    xmax_value = p_value(3)/p_value(1);
end
if(strcmp(distrib,'bin2ordk') == 1)
    xmax_value = p_value(3) - p_value(1);
end
if(strcmp(distrib,'bin3ordk') == 1)
    xmax_value = p_value(3) - p_value(1) + 1;
end
if(strcmp(distrib,'binbinmix') == 1)
    xmax_value = p_value(1)*p_value(3);
end
if(strcmp(distrib,'binchi') == 1 | strcmp(distrib,'binerlang') == 1 | ...
   strcmp(distrib,'binexp') == 1 | strcmp(distrib,'binmax') == 1 | ...
   strcmp(distrib,'binray') == 1 | strcmp(distrib,'binweib') == 1)
    cmax = findcmax(d_indx,distrib,p_value,3);
    p_value(3) = cmax;
end
if(strcmp(distrib,'binchisq') == 1 | strcmp(distrib,'bintrnorm') == 1)
    cmax = findcmax(d_indx,distrib,p_value,2);
    p_value(2) = cmax;
end
if(strcmp(distrib,'bingengam') == 1)
    cmax = findcmax(d_indx,distrib,p_value,4);
    p_value(4) = cmax;
end
if(strcmp(distrib,'binomhyp') == 1)
    xmax_value = min(p_value(2),p_value(3));
end
if(strcmp(distrib,'binomsum') == 1)
    xmax_value = p_value(1) + p_value(3);
end
if(strcmp(distrib,'cannibal') == 1)
    xmax_value = p_value(1) + p_value(2);
end
if(strcmp(distrib,'cantor') == 1)
    xmin_value = 10^(-p_value(1));
    xmax_value = 1 - 10^(-p_value(1));
end
if(strcmp(distrib,'circline') == 1)
    xmax_value = 2*p_value(1);
end
if(strcmp(distrib,'dirmultnom') == 1)
    xmax_value = (exp(gammaln(p_value(2)+p_value(1))-gammaln(p_value(2))-gammaln(p_value(1)+1)));
    xmin_value = 1;
end
if(strcmp(distrib,'diskline') == 1)
    xmax_value = 2*p_value(1);
end
if(strcmp(distrib,'dtrbinom') == 1)
    xmax_value = p_value(1)-p_value(4);
end
if(strcmp(distrib,'dbtrlapl') == 1)
    xmin_value = log(2*p_value(2)); 
    xmax_value = -log(2*p_value(1)); 
end
if(strcmp(distrib,'dbtrlogis') == 1)
    xmin_value = (1-((1-p_value(2))/p_value(2))^p_value(3))/p_value(3); 
    xmax_value = (1-((1-p_value(1))/p_value(1))^p_value(3))/p_value(3); 
end
if(strcmp(distrib,'exthyper') == 1)
    xmin_value = max(0,p_value(5)-p_value(2)); 
    xmax_value = min(p_value(1),p_value(5));
end
if(strcmp(distrib,'extstir2f') == 1)
    xmax_value = min(p_value(1),p_value(2));
end
if(strcmp(distrib,'extstir2c') == 1)
    xmax_value = p_value(1)+p_value(2);
end
if(strcmp(distrib,'fatlife') == 1)
    xmin_value = p_value(2) + 1e-7;
end
if(strcmp(distrib,'fishtip') == 1)
    xmin_value = p_value(1) + 1e-4;
end    
if(strcmp(distrib,'genbeta') == 1)
    xmax_value = p_value(2)/(1-p_value(3))^(1/p_value(1));
end
if(strcmp(distrib,'genfrech') == 1)
    xmin_value = p_value(1) + 1e-7;
end
if(strcmp(distrib,'genfishtip') == 1)
    xmin_value = p_value(1) + 1e-7;
end
if(strcmp(distrib,'genlam') == 1)
    xmin_value = p_value(1)-1/(p_value(2)*p_value(3));
    xmax_value = p_value(1)+1/(p_value(2)*p_value(4));
end 
if(strcmp(distrib,'genlog5') == 1)
    if(p_value(1)>0)
        xmax_value = 1/p_value(1);
    end
    if(p_value(1)<0)
        xmin_value = 1/p_value(1);
    end        
end
if(strcmp(distrib,'genpar2') == 1)
    if(p_value(1) >= 0)
        xmin_value = p_value(2);
    elseif(p_value(1) < 0)
        xmin_value = 0;
        xmax_value = p_value(2) - 1/p_value(1) - eps;
    end
end
if(strcmp(distrib,'genunif') == 1)
    xmax_value = p_value(4)+p_value(1)^(-1/p_value(2));
end
if(strcmp(distrib,'genurn') == 1)
    xmax_value = fix(-p_value(2)/p_value(4))+1;
end
if(strcmp(distrib,'genweib') == 1)
    xmin_value = p_value(1) + 1e-7;
end
if(strcmp(distrib,'genxval') == 1)
    if(p_value(1) > 0)
        xmax_value = p_value(2) + p_value(3)/p_value(1);
    end
    if(p_value(1) < 0)
        xmin_value = p_value(2) + p_value(3)/p_value(1);
    end
end
if(strcmp(distrib,'modhans') == 1)
    u2 = p_value(1);
    u4 = p_value(2);
    u6 = p_value(3);
    u8 = p_value(4);
    del = 9*u8*(5*u4-9*u2^2) - 49*u6^2 + 210*u2*u4*u6 - 125*u4^3;
    c0 = 2*(14*u2*u6^2 - 9*u2*u4*u8 - 5*u4^2*u6) / del;
    c1 = (9*u8*(3*u2^2-u4) + 7*u6^2 - 50*u2*u4*u6 + 25*u4^3) / del;
    c2 = 2*(5*u2*u4^2 -u6*(6*u2^2-u4)) / del;
    c = [c2 c1 c0];
    r = roots(c);
    a2 = r(2);
    a = sqrt(a2)
    xmin_value = -a + eps
    xmax_value = a - eps        
end
if(strcmp(distrib,'harris') == 1)
    xmax_value = xmax_value*(p_value(1) + 1);
end
if(strcmp(distrib,'hoskwal') == 1)
    if(p_value(1) > 0)
        xmax_value = p_value(2)/p_value(1);
    end
end
if(strcmp(distrib,'invgaus3p') == 1)
    xmin_value = p_value(3)+eps;
end
if(strcmp(distrib,'invhyper') == 1)
    if(p_value(1) > 0)
        xmax_value = p_value(3) + p_value(1) - p_value(2);
    end
end
if(strcmp(distrib,'irwinhall') == 1)
    xmin_value = p_value(1) * p_value(2);
    xmax_value = p_value(1) * (p_value(2) + 1);
end
if(strcmp(distrib,'kappa4p') == 1)
    h = p_value(1);
    k = p_value(2);
    tsi = p_value(3);
    alpha = p_value(4);
    eps0 = 1e-4;
    if(h > 0 && k > 0)
        xmin_value = tsi + alpha*(1-h^(-k))/k + eps0;
        xmax_value = tsi + alpha/k - eps0;
    elseif(h > 0 && k == 0)
        xmin_value = tsi + alpha*log(h) + eps0;
    elseif(h > 0 && k < 0)
        xmin_value = tsi + alpha*(1-h^(-k))/k + eps0;
    elseif(h <= 0 && k > 0)
        xmax_value = tsi + alpha/k - eps0;
    elseif(h <= 0 && k < 0)
        xmin_value = tsi + alpha/k + eps0;
    end   
end
if(strcmp(distrib,'kies') == 1)
    xmax_value = p_value(2)-0.0001;
end
if(strcmp(distrib,'lagbeta') == 1)
    xmax_value = 1 / p_value(2);
end
if(strcmp(distrib,'levy') == 1)
    xmin_value = p_value(2)+0.0001;
end
if(strcmp(distrib,'levy2') == 1)
    xmin_value = p_value(1)+1e-8;
end
if(strcmp(distrib,'lmstat') == 1)
    xmax_value = p_value(2) - eps;
end
if(strcmp(distrib,'malik1') == 1)
    xmin_value = p_value(2)*log(p_value(3)) + eps;
end
if(strcmp(distrib,'malik2') == 1)
    xmin_value = p_value(3)^p_value(2) + eps;
end
if(strcmp(distrib,'minnegbin') == 1)
    xmax_value = p_value(1) - 1;
end
if(strcmp(distrib,'modbinom1') == 1)
    xmax_value = fix(p_value(1)/p_value(2));
end
if(strcmp(distrib,'modpois1') == 1)
    xmax_value = fix(1/p_value(2));
end
if(strcmp(distrib,'modstev') == 1)
    xmin_value = fix(p_value(1)/p_value(3)+eps+1);
    xmax_value = min(p_value(1),p_value(2));
end
if(strcmp(distrib,'multinom') == 1)
    xmax_value = (p_value(1)-1) * p_value(2) + 1;
end
if(strcmp(distrib,'mxgenbet') == 1)
    xmax_value = p_value(4) / (1-p_value(5))^(1/p_value(3));
end
if(strcmp(distrib,'mrgenbeta') == 1)
    xmax_value = p_value(4) - 0.0002;
end
if(strcmp(distrib,'nandi') == 1)
    xmax_value = fix(30/(1-p_value(3))^2);                 % heuristc
end
if(strcmp(distrib,'naor') == 1)
    xmax_value = p_value(1)-1;
end
if(strcmp(distrib,'nbinordk') == 1)
    xmin_value = p_value(1)*p_value(3);
end
if(strcmp(distrib,'negmulti') == 1)
    xmax_value = (p_value(1)-1) * p_value(3) + 1;
end
if(strcmp(distrib,'nsuccess') == 1)
    xmax_value = fix((p_value(3)+1)/(p_value(2)+1));
end
if(strcmp(distrib,'occupy') == 1)
    xmin_value = p_value(1)-p_value(2);
end
if(strcmp(distrib,'parabolic') == 1)
    xmin_value = p_value(1)-p_value(2);
    xmax_value = p_value(1)+p_value(2);
end
if(strcmp(distrib,'pareto4') == 1)
    xmin_value = p_value(1)+0.005;
end
if(strcmp(distrib,'pearson5') == 1)
    xmin_value = p_value(3)+0.01;
end
if(strcmp(distrib,'phani') == 1)
    xmax_value = p_value(2)-0.0001;
end
if(strcmp(distrib,'poshyper') == 1)
    xmax_value = min(p_value(1),p_value(3));
end
if(strcmp(distrib,'powlognorm') == 1)
    xmin_value = p_value(3)+0.01;
end
if(strcmp(distrib,'raisdcos') == 1)
    xmin_value = p_value(1)-p_value(2);
    xmax_value = p_value(1)+p_value(2);
end
if(strcmp(distrib,'recip') == 1)
    xmin_value = 1/p_value(1);
    xmax_value = 1;
end
if(strcmp(distrib,'rectline') == 1)
    xmax_value = sqrt(p_value(1)^2+p_value(2)^2);
end
if(strcmp(distrib,'rgp') == 1)
    xmax_value = fix(10/(1-p_value(1)*p_value(2))^2);      % heuristc
end
if(strcmp(distrib,'runs') == 1)
    xmax_value = min(p_value(1)+1,p_value(2));
end
if(strcmp(distrib,'schmeis') == 1)
    xmin_value = sign(-p_value(1))*(abs(-p_value(1)))^p_value(2);
    xmax_value = sign(1-p_value(1))*(abs(1-p_value(1)))^p_value(2);
end
if(strcmp(distrib,'stacy') == 1)
    if(p_value(3) > 0)
        xmin_value = 1e-4;
        xmax_value = p_value(4);
    elseif(p_value(3) < 0)
        xmax_value = -1e-4;
        xmin_value = -p_value(4);
    end
end
if(strcmp(distrib,'stevens') == 1)
    xmax_value = min(p_value(1),p_value(2));
end
if(strcmp(distrib,'strchexp') == 1)
    if(p_value(2) > 0)
        xmin_value = 1e-4;
        xmax_value = p_value(3);
    elseif(p_value(2) < 0)
        xmax_value = -1e-4;
        xmin_value = -p_value(3);
    end
end
if(strcmp(distrib,'tadikab') == 1)
    xmax_value = p_value(1)+p_value(4);
end
if(strcmp(distrib,'trunbur') == 1)
    xmin_value = (((1-p_value(5))^(-1/p_value(2))-1) / p_value(3))^(1/p_value(1));
    xmax_value = (((1-p_value(4))^(-1/p_value(2))-1) / p_value(3))^(1/p_value(1));
end
if(strcmp(distrib,'truncauchy') == 1)
    xmin_value = p_value(3)-p_value(1)*p_value(4);
    xmax_value = p_value(3)+p_value(2)*p_value(4);
end
if(strcmp(distrib,'trundexp') == 1)
    xmin_value = p_value(1)-p_value(3)*p_value(2)*log(2);
    xmax_value = p_value(1)+p_value(3)*p_value(2)*log(2);
end
if(strcmp(distrib,'trgenhlog') == 1)
    xmax_value = (1-((1-p_value(2))/(1+p_value(2)))^p_value(1))/p_value(1);
end
if(strcmp(distrib,'tsallis') == 1)
    xmax_value = 1/(p_value(1)*(1-p_value(2)));
end
if(strcmp(distrib,'vonmises') == 1)
    xmin_value = -pi+p_value(2)+eps;
    xmax_value =  pi+p_value(2)-eps;
end
if(strcmp(distrib,'wrgp') == 1)                                       % heuristic
    if(p_value(3) <= 0.35) 
        xmax_value = max(50,fix(1e4*(p_value(3)/.35)^6));                     % heuristic
    elseif(p_value(3) > 0.35 && p_value(3) <= 1)
        xmax_value = max(50,fix(1e4*(.35/p_value(3))^9));
    elseif(p_value(3) > 1)
        xmax_value = fix(50*(.35/p_value(3)));   
    end
end
if(strcmp(distrib,'wcauchy') == 1)
    xmax_value =  2*pi+p_value(2)-eps;
end
if(strcmp(distrib,'wirelength') == 1)
    xmax_value =  4*p_value(1);
end    
if(strcmp(distrib,'ztpoissum') == 1)
    xmin_value = p_value(2) - 1;
    xmax_value = p_value(1) + 1;
end

% User Inputs Additional Parameters for Special Distributions.  Note
% that these Additional Parameters are all saved in 'c_value'.
c_value = [];
if(strcmp(distrib,'crommelin') + strcmp(distrib,'factmult') + ...
    strcmp(distrib,'multinom') + strcmp(distrib,'negmulti') + ...
    strcmp(distrib,'negfacmulti') + strcmp(distrib,'parnprod1') + ...
    strcmp(distrib,'parnprod2') + strcmp(distrib,'stutgenwar') + ...
    strcmp(distrib,'dirmultnom') == 1)                     % Special cases that take more inputs.
    if(strcmp(distrib,'crommelin') + strcmp(distrib,'multinom') + ...
        strcmp(distrib,'negmulti') == 1)			       % Special cases that sum to 1.
        set(hObjectE16,'String','THE SUM OF ADDITIONAL PARAMETERS MUST = 1.  Hit "Enter" to Continue');
        pause
    end
    if(strcmp(distrib,'parnprod2')== 1)
        set(hObjectE16,'String','THE PRODUCT OF ADDITIONAL PARAMETERS MUST < Xmax.  Hit "Enter" to Continue');
        pause
    end
    if(strcmp(distrib,'crommelin') == 1)
        nterms = p_value(5);
    elseif(strcmp(distrib,'dirmultnom') == 1)
        nterms = p_value(2);
    else
        nterms = p_value(1);
    end
    out_string = char(['INPUT  ' num2str(nterms) ' ADDITIONAL PARAMETER(S) BELOW. Then hit "Enter" to Continue']);
    if(strcmp(distrib,'negfacmulti'))
        out_string = strcat(outstring,' 1st ADDITIONAL PARAMETER MUST >= k !!!  Hit "Enter" to Continue.');
    end
    set(hObjectE16,'String',out_string);
    set(hObjectE20,'Visible','on');
    set(hObjectE20,'BackgroundColor',[1.0 0.969 0.922]);
    pause
    waitfor(hObjectE20,'BackgroundColor','white');
     
    if(strcmp(distrib,'parnprod2') == 1)
        xmin_value = prod(c_value);
    end
end 

set(hObjectE11,'String',num2str(xmax_value));
set(hObjectE12,'String',num2str(xmin_value));

%***************************************************************************************************
% Now, Prepare for Calculation of Selected Random Distribution.

% Re-create Input Data-Structure.
userin.distrib = distrib;
userin.p_value = p_value;
userin.d_indx = d_indx;
userin.xmax_value = xmax_value;
userin.xmin_value = xmin_value;
userin.ftype = ftype;
userin.fnorm = fnorm;
userin.nsamples = nsamples;
userin.name = name;
userin.c_value = c_value;
%userin.iseed = iseed;

% Pre-Allocate.
xhold(1:nsamples) = 0;
xhold2(1:nsamples) = 0;
xhold3(1:nsamples) = 0;
fhold(1:nsamples) = 0;
ghold(1:nsamples) = 0;
jhold(1:nsamples) = 0;
xsort(1:nsamples) = 0;
fsort(1:nsamples) = 0;
gsort(1:nsamples) = 0;
chold(1:nsamples) = 0;


% Call Random Distribution Generator.
%*******************************************************************************************************************    
%
% Generate samples of selected distribution.  Five separate techniques are used here:
%
% 1. Integer (Discrete) Distributions - find CDF value closest to URD & use corresponding
%                                       domain value as output sample for this distribution.
%                                       Discrete distributions are further divided into those
%                                       having closed-form CDFs and those without (need normalization).
% 2. Real Distributions (PPF Inv.)    - use PPF(CDF^(-1)) to find samples of
%                                       Tukey's Generalized Lambda, "G and H", Wakeby, or Davies
%                                       Distributions.
% 3. Real Distributions (Integration) - integrate PDF vector to estimate CDF vector, interpolate 
%                                       between CDF values closest to URD & use interpolated domain value
%                                       as output sample for this distribution. "Graphical Inversion".
% 4. Real Distributions (Inversion)   - invert CDF - URD = 0 to find corresponding domain value
%                                       for use as output sample for this distribution. This
%                                       is the standard technique for sampling most distributions.
%
%*******************************************************************************************************************        
%
if(ftype <= 4)
    [xhold,fhold,ghold,jhold,xsort,fsort,gsort,chold,isint] =  gen_distrib2(userin,iseed);
%
%***********************************************************************************************
%
% 5. Combination Distributions        - generate random numbers from distributions having neither
%                                       CDF, PDF, nor PPF in closed form by generating them as 
%                                       combination of simple distributions.
%                                       '50' - indicates function of uniform random distrbution only.
%                                       '51' - indicates function of 1 non-uniform distribution.
%                                       '52' - indicates function of 2 non-uniform distributions.
%                                       '53' - indicates function of 3 non-uniform distributions.                                       
%
%***********************************************************************************************

elseif(ftype >= 50)
    tic
    %rand('state',iseed);
    isint = 0;
    name5 = char(strcat(distrib,'_rnd'));
    hfunc = str2func(name5);
    if(ftype == 51)
        [xhold,fhold,ghold,jhold,xsort,fsort,gsort,chold,isint] = ...
            hfunc(userin,p_value,nsamples);
    end
    if(ftype == 53)
        [xhold,xhold2,xhold3,fhold,ghold,jhold,xsort,fsort,gsort,chold,isint] = ...
            hfunc(userin,p_value,nsamples);
    end
    toc
end

%****************************************************************************************************************************

return



   
         
         
           
               




           
           

           

        

                

Contact us