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

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