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

betasum_pdf(x, p1, q1, p2, q2)
% betasum_pdf.m - evaluates a Beta Sum Probability Density.
%   See "Continuous Univariate Distributions", Johnson, Kotz, and
%   Balakrishnan, V2, J. Wiley, , p.263, 1996. 
%
%               Vector Form of PDF !!!
%
%  Created by Jim Huntley,  11/10/06
%

function [pdf] = betasum_pdf(x, p1, q1, p2, q2)

%persistent ap b1p b2p cp am b1m b2m cm
%persistent Bstarm Bstarp Bratio
%persistent gamap gamcp gamcmap gamam gamcm gamcmam
%persistent tmin tmax nt delt t

% Initializations.
%if(isempty(ap))
    ap = q1;
    b1p = 1-p1;
    b2p = 1-p2;
    cp = q1+q2;
    am = p2;
    b1m = 1-q1;
    b2m = 1-q2;
    cm = p1+p2;
    Bstarm = gammaln(p1+q1) + gammaln(p2+q2) - (gammaln(q1)+gammaln(q2)+gammaln(p1+p2));
    Bstarp = gammaln(p2+q2) + gammaln(p1+q1) - (gammaln(p2)+gammaln(p1)+gammaln(q2+q1));
    Bratio = beta(p1+q2-1+eps,p2+q1-1+eps) / (beta(p1,q1) * beta(p2,q2));
    gamap = gammaln(ap);
    gamcp = gammaln(cp);
    gamcmap = gammaln(cp-ap);
    gamam = gammaln(am);
    gamcm = gammaln(cm);
    gamcmam = gammaln(cm-am);
    sx = size(x,2);
    tol = 1e-8;
    trace = [];
    %warning off all;
    tmin = 0;
    tmax = 1;
    nt = 1000;
    delt = 1 / (nt-1);
    t = 0:delt:1;
%end

sx = size(x,2);
tol = 1e-8;
trace = [];
%warning off all;
% Evaluate and plot Beta Difference PDF.
for jx = 1:sx
    xj = x(jx);
    if(xj >= 0 && xj < 1)
        x1 = xj / (xj-1);
        x2 = xj;
        pdf(jx) = quadl(@kfun1,tmin,tmax,tol,trace,am,b1m,b2m,cm,x1,x2,p1,q1,p2,q2,Bstarm,gamam,gamcm,gamcmam);
    elseif(xj == 1)
        pdf(jx) = Bratio;
    elseif(xj > 1 && xj <= 2)
        x1 = 2 - xj;
        x2 = x1 / (1-xj);
        pdf(jx) = quadl(@kfun2,tmin,tmax,tol,trace,ap,b1p,b2p,cp,x1,x2,p1,q1,p2,q2,Bstarp,gamap,gamcp,gamcmap);
    end
end

return


function [k] = kfun1(t,a,b1,b2,c,x1,x2,p1,q1,p2,q2,Bstar,gama,gamc,gamcma)

%k = Bstar .* gamc .* x2.^(p1+p2-1) .* (1-x2).^(q1-1) .* t.^(a-1) .* (1-t).^(c-a-1) ...
%    .* (1-t.*x1).^(-b1) .* (1-t.*x2).^(-b2) ./ (gama .* gamcma);
t1 = x2.^(p1+p2-1);
t2 = (1-x2).^(q1-1);
t3 = t.^(a-1);
t4 = (1-t).^(c-a-1);
t5 = (1-t.*x1).^(-b1);
t6 = (1-t.*x2).^(-b2);
tp = max(eps, t1.*t2.*t3.*t4.*t5.*t6);
k = exp(Bstar + gamc +log(tp) - (gama + gamcma));
                  
return

function [k] = kfun2(t,a,b1,b2,c,x1,x2,p1,q1,p2,q2,Bstar,gama,gamc,gamcma)

%k =  Bstar .* gamc .* (1-x1).^(p2-1) .* x1.^(q1+q2-1) .* t.^(a-1) .* (1-t).^(c-a-1) ...
%    .* (1-t.*x1).^(-b1) .* (1-t.*x2).^(-b2) ./ (gama .* gamcma);
t1 = (1-x1).^(p2-1);
t2 = x1.^(q1+q2-1);
t3 = t.^(a-1);
t4 = (1-t).^(c-a-1);
t5 = (1-t.*x1).^(-b1);
t6 = (1-t.*x2).^(-b2);
tp = max(eps, t1.*t2.*t3.*t4.*t5.*t6);
k = exp(Bstar + gamc +log(tp) - (gama + gamcma));
                  
return

Contact us