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

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