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

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