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

distribs()
% distribs.m - calling 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,  09/25/06.
%
% Calls: 'gen_distrib.m', 'cdf2d_m', 'cdf3d_m', 'pdf2d_m', + '*_cdf' (See "distrib_full" for 
%         list), '*_pdf' (See "distrib_full" for list) + *_rnd. 
%
% Called by:  None.
%
% Inputs:  User-Supplied Distribution Selection and Parameters -
%           as specified in spreadsheet, 'noise_distrib21.xls'. 
%           Interactive Program with ASCII interface for User Selections.
%
% Outputs: Plots of "time series", CDF, PDF, and sample distribution (histogram).
%          Diagnostic output in Matlab workspace.   ASCII time series output can be
%          saved to 'user_named.dat'
%

function [] = distribs()

%************************************************************************************************************
% Initializations.
keyin = '';

% Main Loop over User Options.
while(isempty(keyin))
    
    % Get Distribution Database from MS-Excel Spreadsheet.
    [AA, BB] = xlsread('noise_distribs21.xls');

    distrib = BB(2:end,1);
    ndistrib = size(distrib,1)-6;
    distrib_full = BB(2:end,2);

    % Distribution Parameters.
    params1 = BB(2:end,3);
    params2 = BB(2:end,4);
    params3 = BB(2:end,5);
    params4 = BB(2:end,6);
    params5 = BB(2:end,7);
    params = [params1 params2 params3 params4 params5];         % Distributions currently limited to 5 parameters max.
    nparams = size(params,2);

    % Distribution Support Limits
    fmin = BB(2:end,8);
    fmax = BB(2:end,9);
    xlim = [fmin fmax];
    xmin = AA(1:end,1);
    xmax = AA(1:end,2);

    % Solution Types, Normalization, (Simple) Special Limits, & Warnings.
    ftype = AA(1:end,3);
    fnorm = AA(1:end,4);
    fxmin = AA(1:end,5);
    fxmax = AA(1:end,6);
    iflong = AA(1:end,7);

%**************************************************************************************************************    
    
    % Solicit User Inputs.
    % First, Print a List of the Distributions to the Screen.
    fprintf('    ************************************************************************************************************************\n');
    fprintf('    * WELCOME TO THE RANDOM NUMBER GENERATOR!!!  Select by number from the following Distributions.                        *\n');
    fprintf('    ************************************************************************************************************************\n');
    for jd = 1:4:ndistrib                                   % Set up user selection screen with 4 columns of distributions.
        name1 = char(distrib_full(jd));
        nblanks = 30 - length(name1);                       % Max length of distribution name is 29 characters.
        if(jd < 10)
            nblanks = nblanks + 1;
        end
        if(jd + 1 <= ndistrib)
            name2 = char(distrib_full(jd+1));
            if(jd+1 < 10)
                nblanks = nblanks + 1;
            end
            if(jd+1 > 99)
                nblanks = nblanks - 1;
            end
        end
        if(jd + 2 <= ndistrib)
            nblanks2 = 30 - length(name2);
            if(jd+2 < 10)
                nblanks2 = nblanks2 + 1;
            end 
            if(jd+2 > 99)
                nblanks2 = nblanks2 - 1;
            end
            name3 = char(distrib_full(jd+2));
        end
        if(jd+3 <= ndistrib)
            nblanks3 = 30 - length(name3);
            if(jd+3 < 10)
                nblanks3 = nblanks3 + 1;
            end 
            if(jd+3 > 99 && jd+2 > 99)
                nblanks3 = nblanks3 - 1;
            end
            name4 = char(distrib_full(jd+3));            
            fprintf(' %d  %s %s %d  %s %s %d  %s %s %d  %s\n',jd,name1,blanks(nblanks),jd+1,name2,blanks(nblanks2),jd+2,name3,blanks(nblanks3),jd+3,name4);
        elseif(jd+1 > ndistrib)
            fprintf(' %d  %s\n',jd, name1);
        elseif(jd+2 > ndistrib)
            fprintf(' %d  %s %s %d  %s\n',jd, name1,blanks(nblanks),jd+1,name2);
        elseif(jd+3 > ndistrib)
            fprintf(' %d  %s %s %d  %s %s %d  %s\n',jd,name1,blanks(nblanks),jd+1,name2,blanks(nblanks2),jd+2,name3);    
        end
    end
    fprintf(' ****************************************************************************************************************************\n');
    
    % Next, Solicit User's Choice of Distribution
    d_indx = [];
    while(isempty(d_indx))
        d_indx = input('\n Enter the Number of the desired Distribution > ');
        if(d_indx < 1 || d_indx > ndistrib)
            d_indx = [];
        end
    end
    name = char(distrib_full(d_indx));
    fprintf('\n You chose %d, the %s Distribution. \n', d_indx, name);
    fprintf('\n ****************************************************************************\n\n');
    pcount = 0;
    for jp = 1:nparams
        if(~isempty(params{d_indx,jp}))            
            name_p = char(params(d_indx,jp));
            pcount = pcount + 1;
            name_v{jp} = name_p;
        end
    end
    
    % Next, Input Bounds on the Support of the Distribution.
    % Print the formal Bounds for this Distribution.
    fprintf('\n\n Each distribution requires Maximum & Minimum values for its support region.\n');
    name_l1 = char(xlim(d_indx,1));
    name_l2 = char(xlim(d_indx,2));
    fprintf('\n The Formal Support for the %s Distribution is:  Xmax = %s and Xmin = %s.\n',name, name_l2, name_l1);
    fprintf('\n *************************************************************************************************');
    fprintf('\n * The User MAY specify their own support as long as the Formal Limits are not exceeded.   *');
    fprintf('\n * The program will run faster if the MAGNITUDES of both support Limits are as small as possible. *');
    fprintf('\n * For example, if the Formal Limit of Xmax is "inf", the the user might try                     *');
    fprintf('\n * Xmax = 10, 100 or even 1000 - whichever support limit will contain "all" the output samples.   *');
    fprintf('\n * NOTE: User inputs for support Limits will be ignored if the Distribution Limits are already    *');
    fprintf('\n *       defined by its input parameters.                                                        *');
    fprintf('\n *************************************************************************************************');
    
    % Solicit User's Choice of Max and Min.
    fprintf('\n\n Default value for Xmax is %d', xmax(d_indx));
    xmax_value = input('\n Enter a Value for Xmax > ');
    if(isempty(xmax_value))
        xmax_value = xmax(d_indx);
    end
    fprintf('\n Default Value for Xmin is %d',xmin(d_indx));
    xmin_value = input('\n Enter a Value for Xmin > ');
    if(isempty(xmin_value))
        xmin_value = xmin(d_indx);
    end
    
    % Next, User Inputs the Distribution Parameters.
    fprintf('\n The %s Distribution requires value(s) for %d input parameter(s)', name, pcount);
    pv = [];
    if(pcount > 0)
        p_value(1:pcount) = 0;
        fprintf('\n Parameter Name(s) is (are):\n');
        for jp = 1:pcount
            pv = [];  
            while(isempty(pv))
                fprintf('\n %s', name_v{jp});
                pv = input('\n Enter a Value for this parameter > ');
            end
            p_value(jp) = pv;
        end
    end
    % 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{d_indx},'nc_chisq') == 1)							
        fprintf('\n The %s Distribution requires additional value(s) for %d input parameter(s).\n', name, p_value(1));
        sum = 0;
        for jp = 1:p_value(1)
            pv = [];  
            while(isempty(pv))
                name_m = ['m(',num2str(jp),')'];
                fprintf('\n %s', name_m);
                pv = input('\n Enter a Value for this additional parameter > ');
                sum = sum + pv^2;                
            end
        end
        p_value(3) = sqrt(sum);
    end   
    if(isempty(pv))
        p_value = [];
    end
    
    % Overwrite User's Choices for Support Limits because Max and Min are directly Specified
    % by the Input Parameters.
    if(fxmin(d_indx) > 0)
        xmin_value = sign(fxmin(d_indx))*p_value(abs(fxmin(d_indx)))+1e-7*sign(fxmin(d_indx));
    end
    if(fxmax(d_indx) > 0)
        xmax_value = sign(fxmax(d_indx))*p_value(abs(fxmax(d_indx)))-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{d_indx},'absorp') == 1)
        xmax_value = min(p_value(2),p_value(3));
    end
    if(strcmp(distrib{d_indx},'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{d_indx},'beta1') == 1)
        xmin_value = p_value(1) + 1e-4; 
        xmax_value = p_value(2) - 1e-4; 
    end
    if(strcmp(distrib{d_indx},'bin1ordk') == 1)
        xmax_value = p_value(3)/p_value(1);
    end
    if(strcmp(distrib{d_indx},'bin2ordk') == 1)
        xmax_value = p_value(3) - p_value(1);
    end
    if(strcmp(distrib{d_indx},'bin3ordk') == 1)
        xmax_value = p_value(3) - p_value(1) + 1;
    end
    if(strcmp(distrib{d_indx},'binbinmix') == 1)
        xmax_value = p_value(1)*p_value(3);
    end
    if(strcmp(distrib{d_indx},'binomhyp') == 1)
        xmax_value = min(p_value(2),p_value(3));
    end
    if(strcmp(distrib{d_indx},'binomsum') == 1)
        xmax_value = p_value(1) + p_value(3);
    end
    if(strcmp(distrib{d_indx},'cannibal') == 1)
        xmax_value = p_value(1) + p_value(2);
    end
    if(strcmp(distrib{d_indx},'circline') == 1)
        xmax_value = 2*p_value(1);
    end
    if(strcmp(distrib{d_indx},'diskline') == 1)
        xmax_value = 2*p_value(1);
    end
    if(strcmp(distrib{d_indx},'dtrbinom') == 1)
        xmax_value = p_value(1)-p_value(4);
    end
    if(strcmp(distrib{d_indx},'dbtrlapl') == 1)
        xmin_value = log(2*p_value(2)); 
        xmax_value = -log(2*p_value(1)); 
    end
    if(strcmp(distrib{d_indx},'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{d_indx},'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{d_indx},'extstir2f') == 1)
        xmax_value = min(p_value(1),p_value(2));
    end
    if(strcmp(distrib{d_indx},'extstir2c') == 1)
        xmax_value = p_value(1)+p_value(2);
    end
    if(strcmp(distrib{d_indx},'fatlife') == 1)
        xmin_value = p_value(2) + 1e-7;
    end
    if(strcmp(distrib{d_indx},'fishtip') == 1)
        xmin_value = p_value(1) + 1e-4;
    end    
    if(strcmp(distrib{d_indx},'genbeta') == 1)
        xmax_value = p_value(2)/(1-p_value(3))^(1/p_value(1));
    end
    if(strcmp(distrib{d_indx},'genfishtip') == 1)
        xmin_value = p_value(1) + 1e-7;
    end
    if(strcmp(distrib{d_indx},'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{d_indx},'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{d_indx},'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{d_indx},'genunif') == 1)
        xmax_value = p_value(4)+p_value(1)^(-1/p_value(2));
    end
    if(strcmp(distrib{d_indx},'genurn') == 1)
        xmax_value = fix(-p_value(2)/p_value(4))+1;
    end
    if(strcmp(distrib{d_indx},'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{d_indx},'gurland') == 1)
        xmax_value = p_value(2) + p_value(3) - 1;
    end
    if(strcmp(distrib{d_indx},'hansmann') == 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);
        b2 = r(1);
        b = sqrt(b2);
        xmin_value = -b + eps;
        xmax_value = b - eps;        
    end
    if(strcmp(distrib{d_indx},'harris') == 1)
        xmax_value = xmax_value*(p_value(1) + 1);
    end
    if(strcmp(distrib{d_indx},'hoskwal') == 1)
        if(p_value(1) > 0)
            xmax_value = p_value(2)/p_value(1);
        end
    end
    if(strcmp(distrib{d_indx},'invgaus3p') == 1)
        xmin_value = p_value(3)+eps;
    end
    if(strcmp(distrib{d_indx},'invhyper') == 1)
        if(p_value(1) > 0)
            xmax_value = p_value(3) + p_value(1) - p_value(2);
        end
    end
    if(strcmp(distrib{d_indx},'irwinhall') == 1)
        xmin_value = p_value(1) * p_value(2);
        xmax_value = p_value(1) * (p_value(2) + 1);
    end
    if(strcmp(distrib{d_indx},'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{d_indx},'kies') == 1)
        xmax_value = p_value(2)-0.0001;
    end
    if(strcmp(distrib{d_indx},'lagbeta') == 1)
        xmax_value = 1 / p_value(2);
    end
    if(strcmp(distrib{d_indx},'levy') == 1)
        xmin_value = p_value(2)+0.0001;
    end
    if(strcmp(distrib{d_indx},'malik1') == 1)
        xmin_value = p_value(2)*log(p_value(3)) + eps;
    end
    if(strcmp(distrib{d_indx},'malik2') == 1)
        xmin_value = p_value(3)^p_value(2) + eps;
    end
    if(strcmp(distrib{d_indx},'minnegbin') == 1)
        xmax_value = p_value(1) - 1;
    end
    if(strcmp(distrib{d_indx},'modbinom1') == 1)
        xmax_value = fix(p_value(1)/p_value(2));
    end
    if(strcmp(distrib{d_indx},'modpois1') == 1)
        xmax_value = fix(1/p_value(2));
    end
    if(strcmp(distrib{d_indx},'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{d_indx},'multinom') == 1)
        xmax_value = (p_value(1)-1) * p_value(2);
    end
    if(strcmp(distrib{d_indx},'nandi') == 1)
        xmax_value = fix(30/(1-p_value(3))^2);                 % heuristc
    end
    if(strcmp(distrib{d_indx},'naor') == 1)
        xmax_value = p_value(1)-1;
    end
    if(strcmp(distrib{d_indx},'nbinordk') == 1)
        xmin_value = p_value(1)*p_value(3);
    end
    if(strcmp(distrib{d_indx},'negmulti') == 1)
        xmax_value = (p_value(1)-1) * p_value(3) + 1;
    end
    if(strcmp(distrib{d_indx},'nsuccess') == 1)
        xmax_value = fix((p_value(3)+1)/(p_value(2)+1));
    end
    if(strcmp(distrib{d_indx},'occupy') == 1)
        xmin_value = p_value(1)-p_value(2);
    end
    if(strcmp(distrib{d_indx},'parabolic') == 1)
        xmin_value = p_value(1)-p_value(2);
        xmax_value = p_value(1)+p_value(2);
    end
    if(strcmp(distrib{d_indx},'pareto4') == 1)
        xmin_value = p_value(1)+0.005;
    end
    if(strcmp(distrib{d_indx},'pearson5') == 1)
        xmin_value = p_value(3)+0.01;
    end
    if(strcmp(distrib{d_indx},'phani') == 1)
        xmax_value = p_value(2)-0.0001;
    end
    if(strcmp(distrib{d_indx},'poshyper') == 1)
        xmax_value = min(p_value(1),p_value(3));
    end
    if(strcmp(distrib{d_indx},'powlognorm') == 1)
        xmin_value = p_value(3)+0.01;
    end
    if(strcmp(distrib{d_indx},'raisdcos') == 1)
        xmin_value = p_value(1)-p_value(2);
        xmax_value = p_value(1)+p_value(2);
    end
    if(strcmp(distrib{d_indx},'recip') == 1)
        xmin_value = 1/p_value(1);
        xmax_value = 1;
    end
    if(strcmp(distrib{d_indx},'rectline') == 1)
        xmax_value = sqrt(p_value(1)^2+p_value(2)^2);
    end
    if(strcmp(distrib{d_indx},'rgp') == 1)
        xmax_value = fix(10/(1-p_value(1)*p_value(2))^2);      % heuristc
    end
    if(strcmp(distrib{d_indx},'runs') == 1)
        xmax_value = min(p_value(1)+1,p_value(2));
    end
    if(strcmp(distrib{d_indx},'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{d_indx},'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{d_indx},'stevens') == 1)
        xmax_value = min(p_value(1),p_value(2));
    end
    if(strcmp(distrib{d_indx},'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{d_indx},'stutgenwar') == 1)
        p_value(4) = xmax_value ;
    end
    if(strcmp(distrib{d_indx},'tadikab') == 1)
        xmax_value = p_value(1)+p_value(4);
    end
    if(strcmp(distrib{d_indx},'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{d_indx},'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{d_indx},'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{d_indx},'trgenhlog') == 1)
        xmax_value = (1-((1-p_value(2))/(1+p_value(2)))^p_value(1))/p_value(1);
    end
    if(strcmp(distrib{d_indx},'tsallis') == 1)
        xmax_value = 1/(p_value(1)*(1-p_value(2)));
    end
    if(strcmp(distrib{d_indx},'vonmises') == 1)
        xmin_value = -pi+p_value(2)+eps;
        xmax_value =  pi+p_value(2)-eps;
    end
    if(strcmp(distrib{d_indx},'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{d_indx},'wcauchy') == 1)
        xmax_value =  2*pi+p_value(2)-eps;
    end
    if(strcmp(distrib{d_indx},'wirelength') == 1)
        xmax_value =  4*p_value(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{d_indx},'crommelin') + strcmp(distrib{d_indx},'factmult') + ...
       strcmp(distrib{d_indx},'multinom') + strcmp(distrib{d_indx},'negmulti') + ...
       strcmp(distrib{d_indx},'negfacmulti') + strcmp(distrib{d_indx},'parnprod1') + ...
       strcmp(distrib{d_indx},'parnprod2') + strcmp(distrib{d_indx},'stutgenwar') == 1) % Special cases that take more inputs.
        if(strcmp(distrib{d_indx},'crommelin') + strcmp(distrib{d_indx},'multinom') + ...
            strcmp(distrib{d_indx},'negmulti') == 1)			      % Special cases that sum to 1.
            fprintf('\n For this Distribution, THE SUM OF ALL ADDITIONAL PARAMETERS MUST = 1 !!!\n');
        end
        if(strcmp(distrib{d_indx},'crommelin') == 1)
            nterms = fix((xmax_value-eps)/p_value(2))+1;
        else
            nterms = p_value(1);
        end
        fprintf('\n The %s Distribution requires additional value(s) for %d input parameter(s).\n', name, nterms);       
        c_value(1:nterms) = 0;
        for jp = 1:nterms
            if(strcmp(distrib{d_indx},'negfacmulti') & jp == 1)
                fprintf('\n For this Distribution, THE 1st ADDITIONAL PARAMETER MUST >= k !!!\n');
            end
            pv = [];  
            while(isempty(pv))
                name_m = ['p(',num2str(jp),')'];
                fprintf('\n %s', name_m);
                pv = input('\n Enter a Value for this additional parameter > ');
                c_value(jp) = pv;
            end
            if(strcmp(distrib{d_indx},'parnprod2') == 1)
                xmin_value = prod(c_value);
            end
        end
    end 
    
    fprintf('\n\n Final value for Xmax is %d', xmax_value);
    fprintf('\n\n Final value for Xmin is %d \n', xmin_value);
    
    % User Inputs the Number of Samples to be Generated.
    % For Some Distributions, Warn if Calculations Could Take a Long Time.
    fprintf('\n **************************************************************************************\n\n');
    if(iflong(d_indx) == 1)    
        fprintf('\n WARNING: The %s Distribution may take a long time to generate a large number of samples!!!\n', name);
    end
    nsamples = [];
    while(isempty(nsamples))
        nsamples = input('\n Enter the Number of Samples to Generate > ');
    end
    fprintf('\n');
  
    % Finally, User Inputs the Seed for the Uniform Random Number Generator.
    fprintf('\n **************************************************************************************\n\n');
    iseed0 = 0;
    fprintf('\n Default value for Seed is %d', iseed0);   
    iseed = input('\n Enter a Value for Random Number Generator Seed > ');
    if(isempty(iseed))
        iseed = iseed0;
    end
    fprintf('\n\n\n');
    %  End of User Inputs.
    
    %***************************************************************************************************
    % Now, Prepare for Calculation of Selected Random Distribution.
    
    % Create Input Data-Structure.
    userin.distrib = distrib{d_indx};
    userin.p_value = p_value;
    userin.d_indx = d_indx;
    userin.xmax_value = xmax_value;
    userin.xmin_value = xmin_value;
    userin.ftype = ftype(d_indx);
    userin.fnorm = fnorm(d_indx);
    userin.nsamples = nsamples;
    userin.name = name;
    userin.c_value = c_value;
    
    % 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(d_indx) <= 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(d_indx) >= 50)
        tic
        rand('state',0);
        isint = 0;
        name5 = char(strcat(distrib{d_indx},'_rnd'));
        if(ftype(d_indx) == 51)
            [xhold,fhold,ghold,jhold,xsort,fsort,gsort,chold,isint] = ...
                feval(name5,userin,p_value,nsamples);
        end
        if(ftype(d_indx) == 53)
            [xhold,xhold2,xhold3,fhold,ghold,jhold,xsort,fsort,gsort,chold,isint] = ...
                feval(name5,userin,p_value,nsamples);
        end
        toc
    end
  
    %***************************************************************************************************************    
    % End of Random Distribution Generation.
    
    % Now, Prepare Outputs for Printing, Plotting and Writing to Disk (user
    % option).
    
    % Print Information for User.
    fprintf('\n\n ****************************************************************');
    fprintf('\n *            RANDOM NUMBER GENERATION COMPLETED!               *');
    fprintf('\n *                                                              *');
    fprintf('\n * Four plots will now be generated for most Distributions:     *');
    fprintf('\n * 1) Raw Samples                                               *');
    fprintf('\n * 2) CDF                                                       *');
    fprintf('\n * 3) PDF                                                       *');
    fprintf('\n * 4) Ordered Samples (Histogram)                               *');
    fprintf('\n *                                                              *');
    fprintf('\n *                       PLEASE WAIT !!!                        *');
    fprintf('\n ****************************************************************');
    
    % Plot Outputs.
    if(ftype(d_indx) == 2)
        [xsort] = sort(fhold);
        xmin_value = xsort(1);
        xmax_value = xsort(end);
        xsort = xmin_value:(xmax_value-xmin_value)/(nsamples-1):xmax_value;
    end
    xdel = (xmax_value - xmin_value)/100;
    if(isint == 1)
        xdel = 1;
    end
    
    % Plot Time series.
    figure(1)
    xx = xhold;
    if(ftype(d_indx) == 2)                  % Don't sort outputs of Inverse PPF functions until histogram plot.
        xx = fhold;
    end

    plot(jhold, xx)                         % Plot time series of noise.
    title('Sample Values')
    xlabel('Sample');
    ylabel('Value');
    if(ftype(d_indx) ~= 2)                  % Don't sort outputs of Inverse PPF functions until histogram plot.
        [xsort, index] = sort(xhold);
    end
    drawnow;
    
    
    % Plot CDF and PDF.
    if(ftype(d_indx) <= 4)                  % Don't plot CDF or PDF for Special Distributions, Recursive Distributions.
        for js = 1:nsamples            
            if(ftype(d_indx) ~= 2)          % Don't sort outputs of Inverse PPF functions until histogram plot.
                fsort(js) = fhold(index(js));   
                if(ftype(d_indx) == 3)      % Don't sort PDF of PDF-only function.
                    fsort(js) = fhold(js);      
                end
                gsort(js) = ghold(index(js));
            end          
        end
        
        figure(2)
        plot(xsort, gsort)				    % Plot CDF of noise.
        title([name ' CDF']);
        xlabel('X');
        ylabel('Cumulative Distribution');
        axis([xmin_value xmax_value 0 1]);
        drawnow;
     
        figure(3)
        xv = xsort;
        if(ftype(d_indx) == 3)
            xv = chold;                     % Use equispaced 'x' for PDF only functions.    
        end
    
        plot(xv, fsort)                     % Plot PDF of noise.
        title([name ' PDF']);
        xlabel('X');
        ylabel('Probability Density');
        ymax = max(fsort);
        if(strcmp(distrib{d_indx},'genxval') == 1)
            ymax = 1;
        end
        axis([xmin_value xmax_value 0 ymax]);
        drawnow;
    end
    
    % Plot Histogram.
    figure(4)
    if(ftype(d_indx) == 1)
        xmin_value = xmin_value - 1.1e-7;
    end
    x = xmin_value:xdel:xmax_value;
    if(ftype(d_indx) == 2)
        [xsort, index] = sort(fhold);
        xdel = (xsort(end)-xsort(1))/100;
        x = xsort(1):xdel:xsort(end);
    end
    if(ftype(d_indx) >= 50)                 % Scale histogram to output data.
        xdel = (xsort(end)-xsort(1))/100;
        if(strcmp(distrib{d_indx},'disades') == 1)
            xdel = 1; 
        end
        x = xsort(1):xdel:xsort(end);
    end
    
    hist(xsort, x);  				        % Plot ordered distribution (histogram).
    title([name ' Sample Distribution']);
    xlabel('X');
    ylabel('Number of Samples');
    drawnow;    

%****************************************************************************************************************************    
    
    % Re-Start or Save Data?
    keyin = input('\n\n\n STOP(Type "1", then "Enter"),  STOP & SAVE(Type "2", then "Enter"),  TRY AGAIN(Type "Enter")? > ');
    
end % Main Loop over User Options.

% Save Output File.
if(keyin == 2)
    filein = input('\n Enter the Name of the File to Write (in single quotes with no extension) > ');
    filename = strcat(filein,'.dat');
    fid = fopen(filename,'w');
    for js = 1:nsamples
        fprintf(fid,'%12.7f\n',xhold(js));
    end
    status = fclose(fid);
    count = fprintf('\n\n ASCII File  %s  written to current directory. \n\n', filename);
end

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

return



   
         
         
           
               




           
           

           

        

                

Contact us