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

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