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

findcmax(d_indx,distrib,p_value,c_indx)
%findcmax.m - find maximum allowed value for 'c' such that all values of
%the pdf >= 0.  Used with Compund Binomial Distributions (Takacs).
%
%   Author: J. Huntley, 06/08/11
%

function [cmax] = findcmax(d_indx,distrib,p_value,c_indx)

global hObjectE16 hObjectE20
global c_value

chold = p_value(c_indx);
cc = 0;              %1st guess for maximum value of 'c'.
nmin = 0;
nmax = p_value(1);
nh = nmin:1:nmax;
mp = 0;
delc = 0.001;
delcmin = 1e-10;
cchold = cc;
cdig = 6;
c_value = 0;
%Find upper limit on 'c' by marching outward from 'cmin = 0' 
%and checking for negative values of 'pdf'.
while (mp >= 0)
    p_value(c_indx) = cc;
    pdfv = pdf4d(nh, d_indx, p_value, distrib, c_value);
    mp= min(pdfv);
    delc_old = delc;
    while(mp < 0)
        delc_old = delc;
        delc = 0.5 * delc;                 
        p_value(c_indx) = cc;
        pdfv = pdf4d(nh, d_indx, p_value, distrib, c_value);
        mp= min(pdfv);
        if(delc < delcmin | mp > 0)
            delc;
            cchold = cc;
            break
        end 
        cc = cc - delc_old + delc;
    end
    cc = cc + delc;
end
cmax = fix(10^cdig * cchold) / 10^cdig;     %Truncate 'cmax' to 'cdig' digits to ensure 'pdf' > 0.

%If "c" > cmax, enter a New Value for "c".
if(chold > cmax)
    %fprintf('\n MAXIMUM ALLOWED VALUE FOR "c" is %d',cmax);
    out_string = char(['MAXIMUM ALLOWED VALUE FOR "c" is  ' num2str(cmax) ' Hit "Enter" to Continue.']);
    set(hObjectE16,'String',out_string);
    pause
    pv = [];  
    while(isempty(pv) | pv > cmax)
        %pv = input('\n Enter a New Value for "c" > 0 & <= cmax ');
        out_string = 'Enter a New Value for "c" > 0 & <= cmax ';
        set(hObjectE16,'String',out_string);
        set(hObjectE20,'Visible','on');
        set(hObjectE20,'BackgroundColor',[1.0 0.969 0.922]);
        pause
        waitfor(hObjectE20,'BackgroundColor','white');
        pv = c_value;
    end
    cmax = pv;
else 
    cmax = chold;
    p_value(c_indx) = chold;
end    

return

Contact us