stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Created namespace for gamrnd routines. 

Also fixed implementation of Jonk's algorithm.
time-shift
Stéphane Adjemia (Scylla) 2018-12-12 15:32:24 +01:00 committed by Stéphane Adjemian (Hermes)
parent 71667846e5
commit 426c63a735
Signed by untrusted user who does not match committer: stepan
GPG Key ID: A6D44CB9C64CE77B
9 changed files with 493 additions and 283 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

63
matlab/distributions/+gamrnd/ahrens_dieter.m Normal file
View File

@ -0,0 +1,63 @@
function g = ahrens_dieter(a, b)
% Returns gamma variates, see Devroye (1986) page 410.
%
% INPUTS
% - a [double] n*1 vector, first hyperparameter.
% - b [double] n*1 vector, second hyperparameter.
%
% OUTPUTS
% - g [double] n*1 vector, gamma variates.
% Copyright (C) 2006-2018 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
nn = length(a);
mm = nn;
bb = (exp(1)+a)/exp(1);
cc = 1./a;
INDEX = 1:mm;
index = INDEX;
UW = NaN(nn,2);
V = NaN(nn,1);
X = NaN(nn,1);
while mm
UW(index,:) = rand(mm,2);
V(index) = UW(index,1).*bb(index);
state1 = find(V(index)<=1);
state2 = find(V(index)>1);
ID = [];
if ~isempty(state1)
X(index(state1)) = V(index(state1)).^cc(index(state1));
ID = INDEX(index(state1(UW(index(state1),2)>exp(-X(index(state1))))));
end
if ~isempty(state2)
X(index(state2)) = -log(cc(index(state2)).*(bb(index(state2))-V(index(state2))));
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
g = X.*b ;

53
matlab/distributions/+gamrnd/berman.m Normal file
View File

@ -0,0 +1,53 @@
function g = berman(a,b)
% Returns gamma variates, see Devroye (1986) page 418.
%
% INPUTS
% - a [double] n*1 vector, first hyperparameter.
% - b [double] n*1 vector, second hyperparameter.
%
% OUTPUTS
% - g [double] n*1 vector, gamma variates.
% Copyright (C) 2006-2018 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
nn = length(a);
mm = nn;
aa = 1./a ;
cc = 1./(1-a) ;
INDEX = 1:mm;
index = INDEX;
UV = NaN(nn,2);
X = NaN(nn,1);
Y = NaN(nn,1);
while mm
UV(index,:) = rand(mm,2);
X(index) = UV(index,1).^aa(index);
Y(index) = UV(index,2).^cc(index);
id = find(X+Y>1);
if isempty(id)
mm = 0;
else
index = INDEX(id);
mm = length(index);
end
end
Z = gamrnd(2*ones(nn,1), ones(nn,1));
g = (Z.*X).*b ;

64
matlab/distributions/+gamrnd/best_1978.m Normal file
View File

@ -0,0 +1,64 @@
function g = best_1978(a ,b)
% Returns gamma variates, see Devroye (1986) page 410.
%
% INPUTS
% - a [double] n*1 vector, first hyperparameter.
% - b [double] n*1 vector, second hyperparameter.
%
% OUTPUTS
% - g [double] n*1 vector, gamma variates.
% Copyright (C) 2006-2018 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
nn = length(a);
mm = nn;
bb = a-1;
cc = 3*a-.75;
UV = NaN(nn,2);
Y = NaN(nn,1);
X = NaN(nn,1);
Z = NaN(nn,1);
W = NaN(nn,1);
index = 1:nn;
INDEX = index;
while mm
UV(index,:) = rand(mm,2);
W(index) = UV(index,1).*(1-UV(index,1));
Y(index) = sqrt(cc(index)./W(index)).*(UV(index,1)-.5);
X(index) = bb(index)+Y(index);
jndex = index(X(index)>=0);
Jndex = setdiff(index,jndex);
if ~isempty(jndex)
Z(jndex) = 64*W(jndex).*W(jndex).*W(jndex).*UV(jndex,2).*UV(jndex,2);
kndex = jndex(Z(jndex)<=1-2*Y(jndex).*Y(jndex)./X(jndex));
Kndex = setdiff(jndex, kndex);
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);
end
g = X.*b;

82
matlab/distributions/+gamrnd/best_1983.m Normal file
View File

@ -0,0 +1,82 @@
function g = best_1983(a, b)
% Returns gamma variates, see Devroye (1986) page 426.
%
% INPUTS
% - a [double] n*1 vector, first hyperparameter.
% - b [double] n*1 vector, second hyperparameter.
%
% OUTPUTS
% - g [double] n*1 vector, gamma variates.
% Copyright (C) 2006-2018 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
nn = length(a);
mm = nn;
tt = .07 + .75*sqrt(1-a);
bb = 1 + exp(-tt).*a./tt;
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
UW(index,:) = rand(mm,2);
V(index) = UW(index,1).*bb(index);
state1 = find(V(index)<=1);
state2 = find(V(index)>1);
ID = [];
if ~isempty(state1)
X(index(state1)) = tt(index(state1)).*V(index(state1)).^cc(index(state1));
test11 = UW(index(state1),2) <= (2-X(index(state1)))./(2+X(index(state1))) ;
id11 = find(~test11);
if ~isempty(id11)
test12 = UW(index(state1(id11)),2) <= exp(-X(index(state1(id11)))) ;
id12 = find(~test12);
else
id12 = [];
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
g = X.*b;

59
matlab/distributions/+gamrnd/cheng.m Normal file
View File

@ -0,0 +1,59 @@
function g = cheng(a, b)
% Returns gamma variates, see Devroye (1986) page 413.
%
% INPUTS
% - a [double] n*1 vector, first hyperparameter.
% - b [double] n*1 vector, second hyperparameter.
%
% OUTPUTS
% - g [double] n*1 vector, gamma variates.
% Copyright (C) 2006-2018 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
nn = length(a);
mm = nn;
bb = a-log(4);
cc = a+sqrt(2*a-1);
UV = NaN(nn,2);
Y = NaN(nn,1);
X = NaN(nn,1);
Z = NaN(nn,1);
R = NaN(nn,1);
index = 1:nn;
INDEX = index;
while mm
UV(index,:) = rand(mm,2);
Y(index) = a(index).*log(UV(index,2)./(1-UV(index,2)));
X(index) = a(index).*exp(UV(index,2));
Z(index) = UV(index,1).*UV(index,2).*UV(index,2);
R(index) = bb(index) + cc(index).*Y(index)-X(index);
jndex = index(R(index) >= 4.5*Z(index)-(1+log(4.5)));
Jndex = setdiff(index, jndex);
if ~isempty(Jndex)
lndex = Jndex(R(Jndex) >= log(Z(Jndex)));
Lndex = setdiff(Jndex, lndex);
else
Lndex = [];
end
index = INDEX(Lndex);
mm = length(index);
end
g = X.*b;

52
matlab/distributions/+gamrnd/johnk.m Normal file
View File

@ -0,0 +1,52 @@
function g = johnk(a, b)
% Returns gamma variates, see Devroye (1986) page 418.
%
% INPUTS
% - a [double] n*1 vector, first hyperparameter.
% - b [double] n*1 vector, second hyperparameter.
%
% OUTPUTS
% - g [double] n*1 vector, gamma variates.
% Copyright (C) 2006-2018 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
nn = length(a);
mm = nn;
aa = 1./a;
bb = 1./(1-a);
INDEX = 1:mm;
index = INDEX;
UV = NaN(nn,2);
X = NaN(nn,1);
Y = NaN(nn,1);
while mm
UV(index,:) = rand(mm,2);
X(index) = UV(index,1).^aa(index);
Y(index) = UV(index,2).^bb(index);
id = find(X+Y>1);
if isempty(id)
mm = 0;
else
index = INDEX(id);
mm = length(index);
end
end
g = (exprnd(ones(nn,1)).*(X./(X+Y))).*b;

58
matlab/distributions/+gamrnd/knuth.m Normal file
View File

@ -0,0 +1,58 @@
function g = knuth(a, b)
% Returns gamma variates, see Bauwens, Lubrano & Richard (1999) page 316.
%
% INPUTS
% - a [double] n*1 vector, first hyperparameter.
% - b [double] n*1 vector, second hyperparameter.
%
% OUTPUTS
% - g [double] n*1 vector, gamma variates.
% Copyright (C) 2006-2018 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
nn = length(a);
mm = nn;
bb = sqrt(2*a-1);
dd = 1./(a-1);
Y = NaN(nn,1);
X = NaN(nn,1);
INDEX = 1:mm;
index = INDEX;
while mm
Y(index) = tan(pi*rand(mm,1));
X(index) = Y(index).*bb(index) + a(index) - 1 ;
idy1 = find(X(index)>=0);
idn1 = setdiff(index,index(idy1));
if ~isempty(idy1)
test = log(rand(length(idy1),1)) <= ...
log(1+Y(index(idy1)).*Y(index(idy1))) + ...
(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);
end
g = X.*b;

50
matlab/distributions/+gamrnd/weibull_rejection.m Normal file
View File

@ -0,0 +1,50 @@
function g = weibull_rejection(a, b)
% Returns gamma variates, see Devroye (1986) page 415.
%
% INPUTS
% - a [double] n*1 vector, first hyperparameter.
% - b [double] n*1 vector, second hyperparameter.
%
% OUTPUTS
% - g [double] n*1 vector, gamma variates.
% Copyright (C) 2006-2018 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
nn = length(a);
mm = nn;
cc = 1./a ;
dd = a.^(a./(1-a)).*(1-a);
ZE = NaN(nn,2);
X = NaN(nn,1);
INDEX = 1:mm;
index = INDEX;
while mm
ZE(index,:) = exprnd(ones(mm,2));
X(index) = ZE(index,1).^cc(index);
id = find( (ZE(:,1)+ZE(:,2) > dd + X) );
if isempty(id)
mm = 0;
else
index = INDEX(id);
mm = length(index);
end
end
g = X.*b;

295
matlab/missing/stats/gamrnd.m
View File

@ -1,4 +1,5 @@
function rnd = gamrnd(a,b,method)
function rnd = gamrnd(a, b, method)
% This function produces independent random variates from the Gamma distribution.
%
% INPUTS
@ -67,7 +68,7 @@ if (any(a<0)) || (any(b<0)) || (any(a==Inf)) || (any(b==Inf))
error('gamrnd:: Input arguments must be finite and positive!');
end
[tests,integer_idx,double_idx] = isint(a);
[~,integer_idx,double_idx] = isint(a);
number_of_integer_a = length(integer_idx);
number_of_double_a = length(double_idx);
@ -95,7 +96,7 @@ end
if number_of_double_a
if strcmpi(method,'BauwensLubranoRichard')
% Algorithm given in Bauwens, Lubrano & Richard (1999) page 316.
rnd(double_idx) = knuth_algorithm(a(double_idx),b(double_idx));
rnd(double_idx) = gamrnd.knuth(a(double_idx),b(double_idx));
else% Algorithm given in Devroye (1986, chapter 9)
small_idx = find(a(double_idx)<1);
big_idx = find(a(double_idx)>1);
@ -104,301 +105,29 @@ if number_of_double_a
if number_of_small_a
if strcmpi(Devroye.small,'Weibull')
% Algorithm given in Devroye (1986, page 415) [Rejection from the Weibull density]
rnd(double_idx(small_idx)) = weibull_rejection_algorithm(a(double_idx(small_idx)),b(double_idx(small_idx)));
rnd(double_idx(small_idx)) = gamrnd.weibull_rejection(a(double_idx(small_idx)),b(double_idx(small_idx)));
elseif strcmpi(Devroye.small,'Johnk')
% Algorithm given in Devroye (1986, page 418) [Johnk's gamma generator]
rnd(double_idx(small_idx)) = johnk_algorithm(a(double_idx(small_idx)),b(double_idx(small_idx)));
rnd(double_idx(small_idx)) = gamrnd.johnk(a(double_idx(small_idx)),b(double_idx(small_idx)));
elseif strcmpi(Devroye.small,'Berman')
% Algorithm given in Devroye (1986, page 418) [Berman's gamma generator]
rnd(double_idx(small_idx)) = berman_algorithm(a(double_idx(small_idx)),b(double_idx(small_idx)));
rnd(double_idx(small_idx)) = gamrnd.berman(a(double_idx(small_idx)),b(double_idx(small_idx)));
elseif strcmpi(Devroye.small,'GS')
% Algorithm given in Devroye (1986, page 425) [Ahrens and Dieter, 1974]
rnd(double_idx(small_idx)) = ahrens_dieter_algorithm(a(double_idx(small_idx)),b(double_idx(small_idx)));
rnd(double_idx(small_idx)) = gamrnd.ahrens_dieter(a(double_idx(small_idx)),b(double_idx(small_idx)));
elseif strcmpi(Devroye.small,'Best')
% Algorithm given in Devroye (1986, page 426) [Best, 1983]
rnd(double_idx(small_idx)) = best_1983_algorithm(a(double_idx(small_idx)),b(double_idx(small_idx)));
rnd(double_idx(small_idx)) = gamrnd.best_1983(a(double_idx(small_idx)),b(double_idx(small_idx)));
end
end
if number_of_big_a
if strcmpi(Devroye.big,'Cheng')
% Algorithm given in Devroye (1986, page 413) [Cheng's rejection algorithm GB]
rnd(double_idx(big_idx)) = cheng_algorithm(a(double_idx(big_idx)),b(double_idx(big_idx)));
rnd(double_idx(big_idx)) = gamrnd.cheng(a(double_idx(big_idx)),b(double_idx(big_idx)));
elseif strcmpi(Devroye.big,'Best')
% Algorithm given in Devroye (1986, page 410) [Best's rejection algorithm XG]
rnd(double_idx(big_idx)) = best_1978_algorithm(a(double_idx(big_idx)),b(double_idx(big_idx)));
rnd(double_idx(big_idx)) = gamrnd.best_1978(a(double_idx(big_idx)),b(double_idx(big_idx)));
end
end
end
end
function gamma_variates = weibull_rejection_algorithm(a,b)
nn = length(a);
mm = nn;
cc = 1./a ;
dd = a.^(a./(1-a)).*(1-a);
ZE = NaN(nn,2);
X = NaN(nn,1);
INDEX = 1:mm;
index = INDEX;
while mm
ZE(index,:) = exprnd(ones(mm,2));
X(index) = ZE(index,1).^cc(index);
id = find( (ZE(:,1)+ZE(:,2) > dd + X) );
if isempty(id)
mm = 0;
else
index = INDEX(id);
mm = length(index);
end
end
gamma_variates = X.*b;
function gamma_variates = johnk_algorithm(a,b)
nn = length(a);
mm = nn;
aa = 1./a ;
bb = 1./b ;
INDEX = 1:mm;
index = INDEX;
UV = NaN(nn,2);
X = NaN(nn,1);
Y = NaN(nn,1);
while mm
UV(index,:) = rand(mm,2);
X(index) = UV(index,1).^aa(index);
Y(index) = UV(index,2).^bb(index);
id = find(X+Y>1);
if isempty(id)
mm = 0;
else
index = INDEX(id);
mm = length(index);
end
end
gamma_variates = exprnd(ones(nn,1)).*X./(X+Y);
function gamma_variates = berman_algorithm(a,b)
nn = length(a);
mm = nn;
aa = 1./a ;
cc = 1./(1-a) ;
INDEX = 1:mm;
index = INDEX;
UV = NaN(nn,2);
X = NaN(nn,1);
Y = NaN(nn,1);
while mm
UV(index,:) = rand(mm,2);
X(index) = UV(index,1).^aa(index);
Y(index) = UV(index,2).^cc(index);
id = find(X+Y>1);
if isempty(id)
mm = 0;
else
index = INDEX(id);
mm = length(index);
end
end
Z = gamrnd(2*ones(nn,1),ones(nn,1));
gamma_variates = Z.*X.*b ;
function gamma_variates = ahrens_dieter_algorithm(a,b)
nn = length(a);
mm = nn;
bb = (exp(1)+a)/exp(1);
cc = 1./a;
INDEX = 1:mm;
index = INDEX;
UW = NaN(nn,2);
V = NaN(nn,1);
X = NaN(nn,1);
while mm
UW(index,:) = rand(mm,2);
V(index) = UW(index,1).*bb(index);
state1 = find(V(index)<=1);
state2 = find(V(index)>1);%setdiff(index,index(state1));
ID = [];
if ~isempty(state1)
X(index(state1)) = V(index(state1)).^cc(index(state1));
test1 = UW(index(state1),2) <= exp(-X(index(state1))) ;
id1 = find(~test1);
ID = INDEX(index(state1(id1)));
end
if ~isempty(state2)
X(index(state2)) = -log(cc(index(state2)).*(bb(index(state2))-V(index(state2))));
test2 = UW(index(state2),2) <= X(index(state2)).^(a(index(state2))-1);
id2 = find(~test2);
if isempty(ID)
ID = INDEX(index(state2(id2)));
else
ID = [ID,INDEX(index(state2(id2)))];
end
end
mm = length(ID);
if mm
index = ID;
end
end
gamma_variates = X.*b ;
function gamma_variates = best_1983_algorithm(a,b)
nn = length(a);
mm = nn;
tt = .07 + .75*sqrt(1-a);
bb = 1 + exp(-tt).*a./tt;
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
UW(index,:) = rand(mm,2);
V(index) = UW(index,1).*bb(index);
state1 = find(V(index)<=1);
state2 = find(V(index)>1);%setdiff(index,index(state1));
ID = [];
if ~isempty(state1)
X(index(state1)) = tt(index(state1)).*V(index(state1)).^cc(index(state1));
test11 = UW(index(state1),2) <= (2-X(index(state1)))./(2+X(index(state1))) ;
id11 = find(~test11);
if ~isempty(id11)
test12 = UW(index(state1(id11)),2) <= exp(-X(index(state1(id11)))) ;
id12 = find(~test12);
else
id12 = [];
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
gamma_variates = X.*b ;
function gamma_variates = knuth_algorithm(a,b)
nn = length(a);
mm = nn;
bb = sqrt(2*a-1);
dd = 1./(a-1);
Y = NaN(nn,1);
X = NaN(nn,1);
INDEX = 1:mm;
index = INDEX;
while mm
Y(index) = tan(pi*rand(mm,1));
X(index) = Y(index).*bb(index) + a(index) - 1 ;
idy1 = find(X(index)>=0);
idn1 = setdiff(index,index(idy1));
if ~isempty(idy1)
test = log(rand(length(idy1),1)) <= ...
log(1+Y(index(idy1)).*Y(index(idy1))) + ...
(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);
end
gamma_variates = X.*b;
function gamma_variates = cheng_algorithm(a,b)
nn = length(a);
mm = nn;
bb = a-log(4);
cc = a+sqrt(2*a-1);
UV = NaN(nn,2);
Y = NaN(nn,1);
X = NaN(nn,1);
Z = NaN(nn,1);
R = NaN(nn,1);
index = 1:nn;
INDEX = index;
while mm
UV(index,:) = rand(mm,2);
Y(index) = a(index).*log(UV(index,2)./(1-UV(index,2)));
X(index) = a(index).*exp(UV(index,2));
Z(index) = UV(index,1).*UV(index,2).*UV(index,2);
R(index) = bb(index) + cc(index).*Y(index)-X(index);
test1 = (R(index) >= 4.5*Z(index)-(1+log(4.5)));
jndex = index(find(test1));
Jndex = setdiff(index,jndex);
if ~isempty(Jndex)
test2 = (R(Jndex) >= log(Z(Jndex)));
lndex = Jndex(find(test2));
Lndex = setdiff(Jndex,lndex);
else
Lndex = [];
end
index = INDEX(Lndex);
mm = length(index);
end
gamma_variates = X.*b;
function gamma_variates = best_1978_algorithm(a,b)
nn = length(a);
mm = nn;
bb = a-1;
cc = 3*a-.75;
UV = NaN(nn,2);
Y = NaN(nn,1);
X = NaN(nn,1);
Z = NaN(nn,1);
W = NaN(nn,1);
index = 1:nn;
INDEX = index;
while mm
UV(index,:) = rand(mm,2);
W(index) = UV(index,1).*(1-UV(index,1));
Y(index) = sqrt(cc(index)./W(index)).*(UV(index,1)-.5);
X(index) = bb(index)+Y(index);
jndex = index(find(X(index)>=0));
Jndex = setdiff(index,jndex);
if ~isempty(jndex)
Z(jndex) = 64*W(jndex).*W(jndex).*W(jndex).*UV(jndex,2).*UV(jndex,2);
kndex = jndex(find(Z(jndex)<=1-2*Y(jndex).*Y(jndex)./X(jndex)));
Kndex = setdiff(jndex,kndex);
if ~isempty(Kndex)
lndex = Kndex(find(log(Z(Kndex))<=2*(bb(Kndex).*log(X(Kndex)./bb(Kndex))-Y(Kndex))));
Lndex = setdiff(Kndex,lndex);
else
Lndex = [];
end
new_index = INDEX(Lndex);
%mm = length(index);
end
index = union(new_index,INDEX(Jndex));
mm = length(index);
end
gamma_variates = X.*b;
end