Code covered by the BSD License

# Generation of Random Variates

### James Huntley (view profile)

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;
elseif(xj == 1)
pdf(jx) = Bratio;
elseif(xj > 1 && xj <= 2)
x1 = 2 - xj;
x2 = x1 / (1-xj);
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```