From 426c63a735cefa7bc6f8a9de3c9e3a2590df9104 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemia=20=28Scylla=29?= Date: Wed, 12 Dec 2018 15:32:24 +0100 Subject: [PATCH] Created namespace for gamrnd routines. Also fixed implementation of Jonk's algorithm. --- matlab/distributions/+gamrnd/ahrens_dieter.m | 63 ++++ matlab/distributions/+gamrnd/berman.m | 53 ++++ matlab/distributions/+gamrnd/best_1978.m | 64 ++++ matlab/distributions/+gamrnd/best_1983.m | 82 +++++ matlab/distributions/+gamrnd/cheng.m | 59 ++++ matlab/distributions/+gamrnd/johnk.m | 52 +++ matlab/distributions/+gamrnd/knuth.m | 58 ++++ .../distributions/+gamrnd/weibull_rejection.m | 50 +++ matlab/missing/stats/gamrnd.m | 295 +----------------- 9 files changed, 493 insertions(+), 283 deletions(-) create mode 100644 matlab/distributions/+gamrnd/ahrens_dieter.m create mode 100644 matlab/distributions/+gamrnd/berman.m create mode 100644 matlab/distributions/+gamrnd/best_1978.m create mode 100644 matlab/distributions/+gamrnd/best_1983.m create mode 100644 matlab/distributions/+gamrnd/cheng.m create mode 100644 matlab/distributions/+gamrnd/johnk.m create mode 100644 matlab/distributions/+gamrnd/knuth.m create mode 100644 matlab/distributions/+gamrnd/weibull_rejection.m diff --git a/matlab/distributions/+gamrnd/ahrens_dieter.m b/matlab/distributions/+gamrnd/ahrens_dieter.m new file mode 100644 index 000000000..a4d5fac9e --- /dev/null +++ b/matlab/distributions/+gamrnd/ahrens_dieter.m @@ -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 . + +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 ; diff --git a/matlab/distributions/+gamrnd/berman.m b/matlab/distributions/+gamrnd/berman.m new file mode 100644 index 000000000..a3f330d20 --- /dev/null +++ b/matlab/distributions/+gamrnd/berman.m @@ -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 . + +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 ; diff --git a/matlab/distributions/+gamrnd/best_1978.m b/matlab/distributions/+gamrnd/best_1978.m new file mode 100644 index 000000000..8c2d9b9a2 --- /dev/null +++ b/matlab/distributions/+gamrnd/best_1978.m @@ -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 . + +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; diff --git a/matlab/distributions/+gamrnd/best_1983.m b/matlab/distributions/+gamrnd/best_1983.m new file mode 100644 index 000000000..0be8736a0 --- /dev/null +++ b/matlab/distributions/+gamrnd/best_1983.m @@ -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 . + +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; diff --git a/matlab/distributions/+gamrnd/cheng.m b/matlab/distributions/+gamrnd/cheng.m new file mode 100644 index 000000000..8ed154db7 --- /dev/null +++ b/matlab/distributions/+gamrnd/cheng.m @@ -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 . + +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; \ No newline at end of file diff --git a/matlab/distributions/+gamrnd/johnk.m b/matlab/distributions/+gamrnd/johnk.m new file mode 100644 index 000000000..0939446ec --- /dev/null +++ b/matlab/distributions/+gamrnd/johnk.m @@ -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 . + +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; \ No newline at end of file diff --git a/matlab/distributions/+gamrnd/knuth.m b/matlab/distributions/+gamrnd/knuth.m new file mode 100644 index 000000000..4099f5049 --- /dev/null +++ b/matlab/distributions/+gamrnd/knuth.m @@ -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 . + +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; \ No newline at end of file diff --git a/matlab/distributions/+gamrnd/weibull_rejection.m b/matlab/distributions/+gamrnd/weibull_rejection.m new file mode 100644 index 000000000..fdf343fbb --- /dev/null +++ b/matlab/distributions/+gamrnd/weibull_rejection.m @@ -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 . + +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; \ No newline at end of file diff --git a/matlab/missing/stats/gamrnd.m b/matlab/missing/stats/gamrnd.m index af2124753..b80b2e83f 100644 --- a/matlab/missing/stats/gamrnd.m +++ b/matlab/missing/stats/gamrnd.m @@ -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 \ No newline at end of file