Streamlined gamrnd algorithms.
parent
8ac4860150
commit
4206cabe52
|
@ -1,6 +1,6 @@
|
||||||
function g = ahrens_dieter(a, b)
|
function g = ahrens_dieter(a, b)
|
||||||
|
|
||||||
% Returns gamma variates, see Devroye (1986) page 410.
|
% Returns gamma variates, see Devroye (1986) page 425.
|
||||||
%
|
%
|
||||||
% INPUTS
|
% INPUTS
|
||||||
% - a [double] n*1 vector, first hyperparameter.
|
% - a [double] n*1 vector, first hyperparameter.
|
||||||
|
@ -30,34 +30,24 @@ nn = length(a);
|
||||||
mm = nn;
|
mm = nn;
|
||||||
bb = (exp(1)+a)/exp(1);
|
bb = (exp(1)+a)/exp(1);
|
||||||
cc = 1./a;
|
cc = 1./a;
|
||||||
INDEX = 1:mm;
|
index = 1:nn;
|
||||||
index = INDEX;
|
U = NaN(nn,1);
|
||||||
UW = NaN(nn,2);
|
W = NaN(nn,1);
|
||||||
V = NaN(nn,1);
|
V = NaN(nn,1);
|
||||||
X = NaN(nn,1);
|
X = NaN(nn,1);
|
||||||
|
|
||||||
while mm
|
while mm
|
||||||
UW(index,:) = rand(mm,2);
|
U(index) = rand(mm,1);
|
||||||
V(index) = UW(index,1).*bb(index);
|
W(index) = rand(mm,1);
|
||||||
state1 = find(V(index)<=1);
|
V(index) = U(index).*bb(index);
|
||||||
state2 = find(V(index)>1);
|
id1 = index(V(index)<=1);
|
||||||
ID = [];
|
id2 = setdiff(index, id1);
|
||||||
if ~isempty(state1)
|
X(id1) = V(id1).^cc(id1);
|
||||||
X(index(state1)) = V(index(state1)).^cc(index(state1));
|
id3 = id1(W(id1)>exp(-X(id1)));
|
||||||
ID = INDEX(index(state1(UW(index(state1),2)>exp(-X(index(state1))))));
|
X(id2) = -log(cc(id2).*(bb(id2)-V(id2)));
|
||||||
end
|
id4 = id2(W(id2)>X(id2).^(a(id2)-1));
|
||||||
if ~isempty(state2)
|
index = [id3, id4];
|
||||||
X(index(state2)) = -log(cc(index(state2)).*(bb(index(state2))-V(index(state2))));
|
mm = length(index);
|
||||||
if isempty(ID)
|
|
||||||
ID = INDEX(index(state2(UW(index(state2),2)>X(index(state2)).^(a(index(state2))-1))));
|
|
||||||
else
|
|
||||||
ID = [ID, INDEX(index(state2(UW(index(state2),2)>X(index(state2)).^(a(index(state2))-1))))];
|
|
||||||
end
|
|
||||||
end
|
|
||||||
mm = length(ID);
|
|
||||||
if mm
|
|
||||||
index = ID;
|
|
||||||
end
|
|
||||||
end
|
end
|
||||||
|
|
||||||
g = X.*b ;
|
g = X.*b ;
|
||||||
|
|
|
@ -40,13 +40,8 @@ while mm
|
||||||
UV(index,:) = rand(mm,2);
|
UV(index,:) = rand(mm,2);
|
||||||
X(index) = UV(index,1).^aa(index);
|
X(index) = UV(index,1).^aa(index);
|
||||||
Y(index) = UV(index,2).^cc(index);
|
Y(index) = UV(index,2).^cc(index);
|
||||||
id = find(X+Y>1);
|
index = index(X(index)+Y(index)>1);
|
||||||
if isempty(id)
|
mm = length(index);
|
||||||
mm = 0;
|
|
||||||
else
|
|
||||||
index = INDEX(id);
|
|
||||||
mm = length(index);
|
|
||||||
end
|
|
||||||
end
|
end
|
||||||
|
|
||||||
Z = gamrnd(2*ones(nn,1), ones(nn,1));
|
Z = gamrnd(2*ones(nn,1), ones(nn,1));
|
||||||
|
|
|
@ -1,6 +1,6 @@
|
||||||
function g = best_1978(a ,b)
|
function g = best_1978(a ,b)
|
||||||
|
|
||||||
% Returns gamma variates, see Devroye (1986) page 410.
|
% Returns gamma variates, see Devroye (1986) page 410 and Best (1978).
|
||||||
%
|
%
|
||||||
% INPUTS
|
% INPUTS
|
||||||
% - a [double] n*1 vector, first hyperparameter.
|
% - a [double] n*1 vector, first hyperparameter.
|
||||||
|
@ -30,34 +30,24 @@ nn = length(a);
|
||||||
mm = nn;
|
mm = nn;
|
||||||
bb = a-1;
|
bb = a-1;
|
||||||
cc = 3*a-.75;
|
cc = 3*a-.75;
|
||||||
UV = NaN(nn,2);
|
U = NaN(nn,1);
|
||||||
Y = NaN(nn,1);
|
Y = NaN(nn,1);
|
||||||
X = NaN(nn,1);
|
X = NaN(nn,1);
|
||||||
Z = NaN(nn,1);
|
Z = NaN(nn,1);
|
||||||
W = NaN(nn,1);
|
W = NaN(nn,1);
|
||||||
index = 1:nn;
|
index = 1:nn;
|
||||||
INDEX = index;
|
|
||||||
|
|
||||||
while mm
|
while mm
|
||||||
UV(index,:) = rand(mm,2);
|
U(index) = rand(mm,1);
|
||||||
W(index) = UV(index,1).*(1-UV(index,1));
|
W(index) = U(index).*(1-U(index)); % e
|
||||||
Y(index) = sqrt(cc(index)./W(index)).*(UV(index,1)-.5);
|
Y(index) = sqrt(cc(index)./W(index)).*(U(index)-.5); % f
|
||||||
X(index) = bb(index)+Y(index);
|
X(index) = bb(index)+Y(index); % x
|
||||||
jndex = index(X(index)>=0);
|
id1 = index(X(index)<0); % Reject.
|
||||||
Jndex = setdiff(index,jndex);
|
id2 = setdiff(index, id1);
|
||||||
if ~isempty(jndex)
|
Z(id2) = 64.0*(W(id2).^3).*(rand(length(id2),1).^2); % d
|
||||||
Z(jndex) = 64*W(jndex).*W(jndex).*W(jndex).*UV(jndex,2).*UV(jndex,2);
|
id3 = id2(Z(id2)>1.0-2.0*Y(id2).*Y(id2)./X(id2)); % Reject.
|
||||||
kndex = jndex(Z(jndex)<=1-2*Y(jndex).*Y(jndex)./X(jndex));
|
id4 = id3(log(Z(id3))>2.0*(bb(id3).*log(X(id3)./bb(id3))-Y(id3))); % Reject.
|
||||||
Kndex = setdiff(jndex, kndex);
|
index = [id1, id4];
|
||||||
if ~isempty(Kndex)
|
|
||||||
lndex = Kndex(log(Z(Kndex))<=2*(bb(Kndex).*log(X(Kndex)./bb(Kndex))-Y(Kndex)));
|
|
||||||
Lndex = setdiff(Kndex, lndex);
|
|
||||||
else
|
|
||||||
Lndex = [];
|
|
||||||
end
|
|
||||||
new_index = INDEX(Lndex);
|
|
||||||
end
|
|
||||||
index = union(new_index, INDEX(Jndex));
|
|
||||||
mm = length(index);
|
mm = length(index);
|
||||||
end
|
end
|
||||||
|
|
||||||
|
|
|
@ -1,6 +1,6 @@
|
||||||
function g = best_1983(a, b)
|
function g = best_1983(a, b)
|
||||||
|
|
||||||
% Returns gamma variates, see Devroye (1986) page 426.
|
% Returns gamma variates, see Devroye (1986) page 426 and Best (1983) page 187.
|
||||||
%
|
%
|
||||||
% INPUTS
|
% INPUTS
|
||||||
% - a [double] n*1 vector, first hyperparameter.
|
% - a [double] n*1 vector, first hyperparameter.
|
||||||
|
@ -28,55 +28,31 @@ function g = best_1983(a, b)
|
||||||
|
|
||||||
nn = length(a);
|
nn = length(a);
|
||||||
mm = nn;
|
mm = nn;
|
||||||
tt = .07 + .75*sqrt(1-a);
|
index = 1:nn;
|
||||||
bb = 1 + exp(-tt).*a./tt;
|
U = NaN(nn,1);
|
||||||
|
Ustar = NaN(nn, 1);
|
||||||
|
P = NaN(nn,1);
|
||||||
|
X = NaN(nn,1);
|
||||||
|
Y = NaN(nn,1);
|
||||||
|
zz = .07 + .75*sqrt(1-a);
|
||||||
|
bb = 1 + exp(-zz).*a./zz;
|
||||||
cc = 1./a;
|
cc = 1./a;
|
||||||
INDEX = 1:mm;
|
|
||||||
index = INDEX;
|
|
||||||
UW = NaN(nn,2);
|
|
||||||
V = NaN(nn,1);
|
|
||||||
X = NaN(nn,1);
|
|
||||||
Y = NaN(nn,1);
|
|
||||||
|
|
||||||
while mm
|
while mm
|
||||||
UW(index,:) = rand(mm,2);
|
U(index) = rand(mm,1);
|
||||||
V(index) = UW(index,1).*bb(index);
|
Ustar(index) = rand(mm, 1);
|
||||||
state1 = find(V(index)<=1);
|
P(index) = U(index).*bb(index);
|
||||||
state2 = find(V(index)>1);
|
id1 = index(P(index)<=1);
|
||||||
ID = [];
|
id2 = setdiff(index, id1); % Goto 4.
|
||||||
if ~isempty(state1)
|
X(id1) = zz(id1).*(P(id1).^cc(id1));
|
||||||
X(index(state1)) = tt(index(state1)).*V(index(state1)).^cc(index(state1));
|
id3 = id1(Ustar(id1)>((2-X(id1))./(2+X(id1))));
|
||||||
test11 = UW(index(state1),2) <= (2-X(index(state1)))./(2+X(index(state1))) ;
|
id5 = id3(Ustar(id3)>exp(-X(id3)));
|
||||||
id11 = find(~test11);
|
X(id2) = -log(cc(id2).*zz(id2).*(bb(id2)-P(id2))); % This is 4.
|
||||||
if ~isempty(id11)
|
Y(id2) = X(id2)./zz(id2);
|
||||||
test12 = UW(index(state1(id11)),2) <= exp(-X(index(state1(id11)))) ;
|
id4 = id2(Ustar(id2).*(a(id2)+(1-a(id2)).*Y(id2))>1);
|
||||||
id12 = find(~test12);
|
id6 = id4(Ustar(id4)>Y(id4).^(a(id4)-1));
|
||||||
else
|
index = [id5, id6];
|
||||||
id12 = [];
|
mm = length(index);
|
||||||
end
|
|
||||||
ID = INDEX(index(state1(id11(id12))));
|
|
||||||
end
|
|
||||||
if ~isempty(state2)
|
|
||||||
X(index(state2)) = -log(cc(index(state2)).*tt(index(state2)).*(bb(index(state2))-V(index(state2)))) ;
|
|
||||||
Y(index(state2)) = X(index(state2))./tt(index(state2)) ;
|
|
||||||
test21 = UW(index(state2),2).*(a(index(state2)) + Y(index(state2)) - a(index(state2)).*Y(index(state2)) ) <= 1 ;
|
|
||||||
id21 = find(~test21);
|
|
||||||
if ~isempty(id21)
|
|
||||||
test22 = UW(index(state2(id21)),2) <= Y(index(state2(id21))).^(a(index(state2(id21)))-1) ;
|
|
||||||
id22 = find(~test22);
|
|
||||||
else
|
|
||||||
id22 = [];
|
|
||||||
end
|
|
||||||
if isempty(ID)
|
|
||||||
ID = INDEX(index(state2(id21(id22))));
|
|
||||||
else
|
|
||||||
ID = [ID,INDEX(index(state2(id21(id22))))];
|
|
||||||
end
|
|
||||||
end
|
|
||||||
mm = length(ID);
|
|
||||||
if mm
|
|
||||||
index = ID;
|
|
||||||
end
|
|
||||||
end
|
end
|
||||||
|
|
||||||
g = X.*b;
|
g = X.*b;
|
||||||
|
|
|
@ -40,13 +40,8 @@ while mm
|
||||||
UV(index,:) = rand(mm,2);
|
UV(index,:) = rand(mm,2);
|
||||||
X(index) = UV(index,1).^aa(index);
|
X(index) = UV(index,1).^aa(index);
|
||||||
Y(index) = UV(index,2).^bb(index);
|
Y(index) = UV(index,2).^bb(index);
|
||||||
id = find(X+Y>1);
|
index = index(X(index)+Y(index)>1);
|
||||||
if isempty(id)
|
mm = length(index);
|
||||||
mm = 0;
|
|
||||||
else
|
|
||||||
index = INDEX(id);
|
|
||||||
mm = length(index);
|
|
||||||
end
|
|
||||||
end
|
end
|
||||||
|
|
||||||
g = (exprnd(ones(nn,1)).*(X./(X+Y))).*b;
|
g = (exprnd(ones(nn,1)).*(X./(X+Y))).*b;
|
|
@ -1,6 +1,6 @@
|
||||||
function g = knuth(a, b)
|
function g = knuth(a, b)
|
||||||
|
|
||||||
% Returns gamma variates, see Bauwens, Lubrano & Richard (1999) page 316.
|
% Returns gamma variates, see Knuth (1981) page 129.
|
||||||
%
|
%
|
||||||
% INPUTS
|
% INPUTS
|
||||||
% - a [double] n*1 vector, first hyperparameter.
|
% - a [double] n*1 vector, first hyperparameter.
|
||||||
|
@ -29,29 +29,18 @@ function g = knuth(a, b)
|
||||||
nn = length(a);
|
nn = length(a);
|
||||||
mm = nn;
|
mm = nn;
|
||||||
bb = sqrt(2*a-1);
|
bb = sqrt(2*a-1);
|
||||||
dd = 1./(a-1);
|
|
||||||
Y = NaN(nn,1);
|
Y = NaN(nn,1);
|
||||||
X = NaN(nn,1);
|
X = NaN(nn,1);
|
||||||
INDEX = 1:mm;
|
index = 1:mm;
|
||||||
index = INDEX;
|
|
||||||
|
|
||||||
while mm
|
while mm
|
||||||
Y(index) = tan(pi*rand(mm,1));
|
Y(index) = tan(pi*rand(mm,1));
|
||||||
X(index) = Y(index).*bb(index) + a(index) - 1 ;
|
X(index) = Y(
|
||||||
idy1 = find(X(index)>=0);
|
index).*bb(index) + a(index) - 1;
|
||||||
idn1 = setdiff(index,index(idy1));
|
id1 = index(X(index)<=0); % Rejected draws.
|
||||||
if ~isempty(idy1)
|
id2 = setdiff(index, id1);
|
||||||
test = log(rand(length(idy1),1)) <= ...
|
id3 = id2(rand(length(id2), 1)>(1+Y(id2).*Y(id2)).*exp((a(id2)-1).*(log(X(id2))-log(a(id2)-1))-bb(id2).*Y(id2))); % Rejected draws.
|
||||||
log(1+Y(index(idy1)).*Y(index(idy1))) + ...
|
index = [id1, id3];
|
||||||
(a(index(idy1))-1).*log(X(index(idy1)).*dd(index(idy1))) - ...
|
|
||||||
Y(index(idy1)).*bb(index(idy1)) ;
|
|
||||||
idy2 = find(test);
|
|
||||||
idn2 = setdiff(idy1, idy1(idy2));
|
|
||||||
else
|
|
||||||
idy2 = [];
|
|
||||||
idn2 = [];
|
|
||||||
end
|
|
||||||
index = [ INDEX(idn1) , INDEX(index(idn2)) ] ;
|
|
||||||
mm = length(index);
|
mm = length(index);
|
||||||
end
|
end
|
||||||
|
|
||||||
|
|
|
@ -37,7 +37,6 @@ X = NaN(nn, 1);
|
||||||
index = 1:nn;
|
index = 1:nn;
|
||||||
|
|
||||||
while mm
|
while mm
|
||||||
% Generate Weibull
|
|
||||||
Z(index) = -log(rand(mm, 1));
|
Z(index) = -log(rand(mm, 1));
|
||||||
Y(index) = Z(index).^cc(index);
|
Y(index) = Z(index).^cc(index);
|
||||||
INDEX = index(rand(mm,1)>aa(index).*exp(Z(index)-Y(index)));
|
INDEX = index(rand(mm,1)>aa(index).*exp(Z(index)-Y(index)));
|
||||||
|
|
Loading…
Reference in New Issue