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

betadiff_pdf(x, p1, q1, p2, q2)
% betadiff_pdf.m - evaluates a Beta Difference Probability Density.
%   See "Continuous Univariate Distributions", Johnson, Kotz, and
%   Balakrishnan, V2, J. Wiley, , p.263, 1996 + "Bayesian Analysis of the Difference of Two Proportions",
%   T. Pham-Gia + N. Turkkan, Comm. Stat. - Theory Meth., 22(6), p.1755, 1993. 
%
%               Vector Form of PDF !!!
%
%  Created by Jim Huntley,  10/30/06
%
%  WARNING:  PDF may develop discontinuities at x = 0 for large values of p1+q1+p2+q2!!!
%


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

%persistent ap b1p b2p cp am b1m b2m cm B21 B12 A logA
%persistent gamap gamcp gamcmap gamam gamcm gamcmam
%persistent tmin tmax nt delt t

% Initializations.

%if(isempty(ap))
    ap = q1;
    b1p = p1+p2+q1+q2-2;
    b2p = 1-p1;
    cp = p2+q1;
    am = q2;
    b1m = 1-p2;
    b2m = p1+p2+q1+q2-2;
    cm = p1+q2;
    B21 = log(beta(p2,q1));
    B12 = log(beta(p1,q2));
    A = beta(p1,q1) * beta(p2,q2);
    logA = log(A);
    gamap = gammaln(ap);
    gamcp = gammaln(cp);
    gamcmap = gammaln(cp-ap);
    gamam = gammaln(am);
    gamcm = gammaln(cm);
    gamcmam = gammaln(cm-am);
    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)
        x1 = 1 - xj;
        x2 = 1 - xj^2;
        pdf(jx) = quadl(@kfun1,tmin,tmax,tol,trace,ap,b1p,b2p,cp,x1,x2,p1,q1,p2,q2,gamap,gamcp,gamcmap,B21,logA);
    elseif(xj < 0)
        x1 = 1 - xj^2;
        x2 = 1 + xj;
        pdf(jx) = quadl(@kfun2,tmin,tmax,tol,trace,am,b1m,b2m,cm,x1,x2,p1,q1,p2,q2,gamam,gamcm,gamcmam,B12,logA);
    end
end

% Fudge to interpolate across discontinuity at x = 0.
dpdf = diff(pdf);
index = find(abs(dpdf) > 0.25);
if(~isempty(index))
    di = index(2) - index(1); 
    y(1) = pdf(index(1)-1);
    y(2) = pdf(index(1));
    y(3) = beta(p1+p2-1,q1+q2-1) / A;
    y(4) = pdf(index(2)+1);
    y(5) = pdf(index(2)+2);
    xx(1) = x(index(1)-1);
    xx(2) = x(index(1));
    xx(3) = 0;
    xx(4) = x(index(2)+1);
    xx(5) = x(index(2)+2);
    for ji = 1:di
        pdf(index(1)+ji) = interp1(xx,y,x(index(1)+ji),'cubic');
    end
end

return


function [k] = kfun1(t,a,b1,b2,c,x1,x2,p1,q1,p2,q2,gamap,gamcp,gamcmap,B21,logA)

%k = B21 .* gamcp .* (1-x1).^(q1+q2-1) .* x1.^(p2+q1-1) .* t.^(a-1) .* (1-t).^(c-a-1) ...
%    .* (1-t.*x1).^(-b1) .* (1-t.*x2).^(-b2) ./ (gamap .* gamcmap .* A);
t1 = (1-x1).^(q1+q2-1);
t2 = x1.^(p2+q1-1);
t3 = t.^(a-1);
t4 = (1-t).^(c-a-1);
t5 = (1-t.*x1).^(-b1);
t6 = (1-t.*x2).^(-b2);
pt = max(eps,t1 .* t2 .* t3 .* t4 .* t5 .* t6);
k = exp(B21 + gamcp + log(pt) - (gamap + gamcmap + logA));
                  
return

function [k] = kfun2(t,a,b1,b2,c,x1,x2,p1,q1,p2,q2,gamam,gamcm,gamcmam,B12,logA)

%k =  B12 .* gamcm .* (1-x2).^(q1+q2-1) .* x2.^(p1+q2-1) .* t.^(a-1) .* (1-t).^(c-a-1) ...
%    .* (1-t.*x1).^(-b1) .* (1-t.*x2).^(-b2) ./ (gamam .* gamcmam .* A);
t1 = (1-x2).^(q1+q2-1);
t2 = x2.^(p2+q1-1);
t3 = t.^(a-1);
t4 = (1-t).^(c-a-1);
t5 = (1-t.*x1).^(-b1);
t6 = (1-t.*x2).^(-b2);
pt = max(eps,t1 .* t2 .* t3 .* t4 .* t5 .* t6);
k =  exp(B12 + gamcm + log(pt) - (gamam + gamcmam + logA));
                  
return

Contact us