gamrnd

 This function produces independent random variates from the Gamma distribution.

  INPUTS 
    a       [double]    n*1 vector of positive parameters.
    b       [double]    n*1 vector of positive parameters.    
    method  [string]    'BawensLubranoRichard' or anything else (see below).

  OUTPUT 
    rnd     [double]    n*1 vector of independent variates from the gamma(a,b) distribution.
                        rnd(i) is gamma distributed with mean a(i)b(i) and variance a(i)b(i)^2.
  
  ALGORITHMS     
    Described in Bauwens, Lubrano and Richard (1999, page 316) and Devroye (1986, chapter 9).

0001 function rnd = gamrnd(a,b,method)
0002 % This function produces independent random variates from the Gamma distribution.
0003 %
0004 %  INPUTS
0005 %    a       [double]    n*1 vector of positive parameters.
0006 %    b       [double]    n*1 vector of positive parameters.
0007 %    method  [string]    'BawensLubranoRichard' or anything else (see below).
0008 %
0009 %  OUTPUT
0010 %    rnd     [double]    n*1 vector of independent variates from the gamma(a,b) distribution.
0011 %                        rnd(i) is gamma distributed with mean a(i)b(i) and variance a(i)b(i)^2.
0012 %
0013 %  ALGORITHMS
0014 %    Described in Bauwens, Lubrano and Richard (1999, page 316) and Devroye (1986, chapter 9).
0015 
0016 % Copyright (C) 2006-2011 Dynare Team
0017 %
0018 % This file is part of Dynare.
0019 %
0020 % Dynare is free software: you can redistribute it and/or modify
0021 % it under the terms of the GNU General Public License as published by
0022 % the Free Software Foundation, either version 3 of the License, or
0023 % (at your option) any later version.
0024 %
0025 % Dynare is distributed in the hope that it will be useful,
0026 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0027 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0028 % GNU General Public License for more details.
0029 %
0030 % You should have received a copy of the GNU General Public License
0031 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0032 
0033 if (nargin < 2)
0034     error('gamrnd:: Two input arguments are needed!');
0035 end
0036 
0037 if nargin==2
0038     method= 'BauwensLubranoRichard';
0039     if any(a<1)
0040         method = 'Devroye';
0041         Devroye.small = 'Best'; % 'Weibull' , 'Johnk' , 'Berman' , 'GS' , 'Best'
0042                                 % REMARK: The first algorithm (Weibull) is producing too much extreme values.
0043     end
0044     if ~strcmpi(method,'BauwensLubranoRichard')
0045         Devroye.big = 'Best'; % 'Cheng' , 'Best'
0046                               % REMARK 1: The first algorithm (Cheng) is still producing obviously wrong simulations.
0047                               % REMARK 2: The second algorithm seems slightly slower than the algorithm advocated by Bauwens,
0048                               %           Lubrano and Richard, but the comparison depends on the value of a (this should be
0049                               %           investigated further).
0050     end
0051 else
0052     error('gamrnd:: Selection of method not yet implemented')
0053 end
0054 
0055 [ma,na] = size(a);
0056 [mb,nb] = size(b);
0057 
0058 if ma~=mb || na~=nb
0059     error('gamrnd:: Input arguments must have the same size!');
0060 end
0061 
0062 if na~=1
0063     error('gamrnd:: Input arguments must be column vectors');
0064 end
0065 
0066 if (any(a<0)) || (any(b<0)) || (any(a==Inf)) || (any(b==Inf))
0067     error('gamrnd:: Input arguments must be finite and positive!');
0068 end
0069 
0070 [tests,integer_idx,double_idx] = isint(a);
0071 
0072 number_of_integer_a = length(integer_idx);
0073 number_of_double_a = length(double_idx);
0074 
0075 rnd = NaN(ma,1);
0076 
0077 if number_of_integer_a
0078     small_idx = find(a(integer_idx)<30);
0079     big_idx = find(a(integer_idx)>=30);
0080     number_of_small_a = length(small_idx);
0081     number_of_big_a = length(big_idx);
0082     if number_of_small_a
0083         % Exact sampling.
0084         for i=1:number_of_small_a
0085             rnd(integer_idx(small_idx(i))) = sum(exprnd(ones(a(integer_idx(small_idx(i))),1)))*b(integer_idx(small_idx(i)));
0086         end
0087     end
0088     if number_of_big_a
0089         % Gaussian approximation.
0090         rnd(integer_idx(big_idx)) = sqrt(a(integer_idx(big_idx))).* b(integer_idx(big_idx)) .* randn(number_of_big_a, 1) + a(integer_idx(big_idx)) .* b(integer_idx(big_idx));
0091     end
0092 end
0093 
0094 
0095 if number_of_double_a
0096     if strcmpi(method,'BauwensLubranoRichard')
0097         % Algorithm given in Bauwens, Lubrano & Richard (1999) page 316.
0098         rnd(double_idx) = knuth_algorithm(a(double_idx),b(double_idx));
0099     else% Algorithm given in  Devroye (1986, chapter 9)
0100         small_idx = find(a(double_idx)<1);
0101         big_idx = find(a(double_idx)>1);
0102         number_of_small_a = length(small_idx);
0103         number_of_big_a = length(big_idx);
0104         if number_of_small_a
0105             if strcmpi(Devroye.small,'Weibull')
0106                 % Algorithm given in Devroye (1986, page 415) [Rejection from the Weibull density]
0107                 rnd(double_idx(small_idx)) = weibull_rejection_algorithm(a(double_idx(small_idx)),b(double_idx(small_idx)));
0108             elseif strcmpi(Devroye.small,'Johnk')
0109                 % Algorithm given in Devroye (1986, page 418) [Johnk's gamma generator]
0110                 rnd(double_idx(small_idx)) = johnk_algorithm(a(double_idx(small_idx)),b(double_idx(small_idx)));
0111             elseif strcmpi(Devroye.small,'Berman')
0112                 % Algorithm given in Devroye (1986, page 418) [Berman's gamma generator]
0113                 rnd(double_idx(small_idx)) = berman_algorithm(a(double_idx(small_idx)),b(double_idx(small_idx)));
0114             elseif strcmpi(Devroye.small,'GS')
0115                 % Algorithm given in Devroye (1986, page 425) [Ahrens and Dieter, 1974]
0116                 rnd(double_idx(small_idx)) = ahrens_dieter_algorithm(a(double_idx(small_idx)),b(double_idx(small_idx)));
0117             elseif strcmpi(Devroye.small,'Best')
0118                 % Algorithm given in Devroye (1986, page 426) [Best, 1983]
0119                 rnd(double_idx(small_idx)) = best_1983_algorithm(a(double_idx(small_idx)),b(double_idx(small_idx)));
0120             end
0121         end
0122         if number_of_big_a
0123             if strcmpi(Devroye.big,'Cheng')
0124                 % Algorithm given in Devroye (1986, page 413) [Cheng's rejection algorithm GB]
0125                 rnd(double_idx(big_idx)) = cheng_algorithm(a(double_idx(big_idx)),b(double_idx(big_idx)));
0126             elseif strcmpi(Devroye.big,'Best')
0127                 % Algorithm given in Devroye (1986, page 410) [Best's rejection algorithm XG]
0128                 rnd(double_idx(big_idx)) = best_1978_algorithm(a(double_idx(big_idx)),b(double_idx(big_idx)));
0129             end
0130         end
0131     end
0132 end
0133 
0134 
0135 function gamma_variates = weibull_rejection_algorithm(a,b)
0136 nn = length(a);
0137 mm = nn;
0138 cc = 1./a ;
0139 dd = a.^(a./(1-a)).*(1-a);
0140 ZE = NaN(nn,2);
0141 X  = NaN(nn,1);
0142 INDEX = 1:mm;
0143 index = INDEX;
0144 while mm
0145     ZE(index,:) = exprnd(ones(mm,2));
0146     X(index) = ZE(index,1).^cc(index);
0147     id = find( (ZE(:,1)+ZE(:,2) > dd + X) );
0148     if isempty(id)
0149         mm = 0;
0150     else
0151         index = INDEX(id);
0152         mm = length(index);
0153     end
0154 end
0155 gamma_variates = X.*b;
0156 
0157 
0158 function gamma_variates = johnk_algorithm(a,b)
0159 nn = length(a);
0160 mm = nn;
0161 aa = 1./a ;
0162 bb = 1./b ;
0163 INDEX = 1:mm;
0164 index = INDEX;
0165 UV = NaN(nn,2);
0166 X  = NaN(nn,1);
0167 Y  = NaN(nn,1);
0168 while mm
0169     UV(index,:) = rand(mm,2);
0170     X(index) = UV(index,1).^aa(index);
0171     Y(index) = UV(index,2).^bb(index);
0172     id = find(X+Y>1);
0173     if isempty(id)
0174         mm = 0;
0175     else
0176         index = INDEX(id);
0177         mm = length(index);
0178     end
0179 end
0180 gamma_variates = exprnd(ones(nn,1)).*X./(X+Y);
0181 
0182 
0183 function gamma_variates = berman_algorithm(a,b)    
0184 nn = length(a);
0185 mm = nn;
0186 aa = 1./a ;
0187 cc = 1./(1-a) ;
0188 INDEX = 1:mm;
0189 index = INDEX;
0190 UV = NaN(nn,2);
0191 X  = NaN(nn,1);
0192 Y  = NaN(nn,1);
0193 while mm
0194     UV(index,:) = rand(mm,2);
0195     X(index) = UV(index,1).^aa(index);
0196     Y(index) = UV(index,2).^cc(index);
0197     id = find(X+Y>1);
0198     if isempty(id)
0199         mm = 0;
0200     else
0201         index = INDEX(id);
0202         mm = length(index);
0203     end
0204 end
0205 Z = gamrnd(2*ones(nn,1),ones(nn,1));
0206 gamma_variates = Z.*X.*b ;
0207 
0208 
0209 function gamma_variates = ahrens_dieter_algorithm(a,b)    
0210 nn = length(a);
0211 mm = nn;
0212 bb = (exp(1)+a)/exp(1);
0213 cc = 1./a;
0214 INDEX = 1:mm;
0215 index = INDEX;
0216 UW = NaN(nn,2);
0217 V  = NaN(nn,1);
0218 X  = NaN(nn,1);
0219 while mm
0220     UW(index,:) = rand(mm,2);
0221     V(index) = UW(index,1).*bb(index);
0222     state1 = find(V(index)<=1);
0223     state2 = find(V(index)>1);%setdiff(index,index(state1));
0224     ID = [];
0225     if ~isempty(state1)
0226         X(index(state1)) = V(index(state1)).^cc(index(state1));
0227         test1 = UW(index(state1),2) <= exp(-X(index(state1))) ;
0228         id1 = find(~test1);
0229         ID = INDEX(index(state1(id1)));
0230     end
0231     if ~isempty(state2)
0232         X(index(state2)) = -log(cc(index(state2)).*(bb(index(state2))-V(index(state2))));
0233         test2 = UW(index(state2),2) <= X(index(state2)).^(a(index(state2))-1);
0234         id2 = find(~test2);
0235         if isempty(ID)
0236             ID = INDEX(index(state2(id2)));
0237         else
0238             ID = [ID,INDEX(index(state2(id2)))];
0239         end
0240     end
0241     mm = length(ID);
0242     if mm
0243         index = ID;
0244     end
0245 end
0246 gamma_variates = X.*b ;
0247 
0248 
0249 function gamma_variates = best_1983_algorithm(a,b)    
0250 nn = length(a);
0251 mm = nn;
0252 tt = .07 + .75*sqrt(1-a);
0253 bb = 1 + exp(-tt).*a./tt;
0254 cc = 1./a;
0255 INDEX = 1:mm;
0256 index = INDEX;
0257 UW = NaN(nn,2);
0258 V  = NaN(nn,1);
0259 X  = NaN(nn,1);
0260 Y  = NaN(nn,1);
0261 while mm
0262     UW(index,:) = rand(mm,2);
0263     V(index) = UW(index,1).*bb(index);
0264     state1 = find(V(index)<=1);
0265     state2 = find(V(index)>1);%setdiff(index,index(state1));
0266     ID = [];
0267     if ~isempty(state1)
0268         X(index(state1)) = tt(index(state1)).*V(index(state1)).^cc(index(state1));
0269         test11 = UW(index(state1),2) <= (2-X(index(state1)))./(2+X(index(state1))) ;
0270         id11 = find(~test11);
0271         if ~isempty(id11)
0272             test12 = UW(index(state1(id11)),2) <= exp(-X(index(state1(id11)))) ;
0273             id12 = find(~test12);
0274         else
0275             id12 = [];
0276         end
0277         ID = INDEX(index(state1(id11(id12))));
0278     end
0279     if ~isempty(state2)
0280         X(index(state2)) = -log(cc(index(state2)).*tt(index(state2)).*(bb(index(state2))-V(index(state2)))) ;
0281         Y(index(state2)) = X(index(state2))./tt(index(state2)) ;
0282         test21 = UW(index(state2),2).*(a(index(state2)) + Y(index(state2)) - a(index(state2)).*Y(index(state2)) ) <= 1 ;
0283         id21 = find(~test21);
0284         if ~isempty(id21)
0285             test22 = UW(index(state2(id21)),2) <= Y(index(state2(id21))).^(a(index(state2(id21)))-1) ;
0286             id22 = find(~test22);
0287         else
0288             id22 = [];
0289         end
0290         if isempty(ID)
0291             ID = INDEX(index(state2(id21(id22))));
0292         else
0293             ID = [ID,INDEX(index(state2(id21(id22))))];
0294         end
0295     end
0296     mm = length(ID);
0297     if mm
0298         index = ID;
0299     end
0300 end
0301 gamma_variates = X.*b ;
0302 
0303 
0304 function  gamma_variates = knuth_algorithm(a,b)
0305 nn = length(a);
0306 mm = nn;
0307 bb = sqrt(2*a-1);
0308 dd = 1./(a-1);
0309 Y = NaN(nn,1);
0310 X = NaN(nn,1);
0311 INDEX = 1:mm;
0312 index = INDEX;
0313 while mm
0314     Y(index) = tan(pi*rand(mm,1));
0315     X(index) = Y(index).*bb(index) + a(index) - 1 ;
0316     idy1 = find(X(index)>=0);
0317     idn1 = setdiff(index,index(idy1));
0318     if ~isempty(idy1)
0319         test = log(rand(length(idy1),1)) <= ...
0320                log(1+Y(index(idy1)).*Y(index(idy1))) + ...
0321                (a(index(idy1))-1).*log(X(index(idy1)).*dd(index(idy1))) - ...
0322                Y(index(idy1)).*bb(index(idy1)) ;
0323         idy2 = find(test);
0324         idn2 = setdiff(idy1,idy1(idy2));
0325     else
0326         idy2 = [];
0327         idn2 = [];
0328     end
0329     index = [ INDEX(idn1) , INDEX(index(idn2)) ] ;
0330     mm = length(index);
0331 end
0332 gamma_variates = X.*b;
0333 
0334 
0335 function  gamma_variates = cheng_algorithm(a,b)
0336 nn = length(a);
0337 mm = nn;
0338 bb = a-log(4);
0339 cc = a+sqrt(2*a-1);
0340 UV = NaN(nn,2);
0341 Y  = NaN(nn,1);
0342 X  = NaN(nn,1);
0343 Z  = NaN(nn,1);
0344 R  = NaN(nn,1);
0345 index = 1:nn;
0346 INDEX = index;
0347 while mm
0348     UV(index,:) = rand(mm,2);
0349     Y(index) = a(index).*log(UV(index,2)./(1-UV(index,2)));
0350     X(index) = a(index).*exp(UV(index,2));
0351     Z(index) = UV(index,1).*UV(index,2).*UV(index,2);
0352     R(index) = bb(index) + cc(index).*Y(index)-X(index);
0353     test1 = (R(index) >= 4.5*Z(index)-(1+log(4.5)));
0354     jndex = index(find(test1));
0355     Jndex = setdiff(index,jndex);
0356     if ~isempty(Jndex)
0357         test2 = (R(Jndex) >= log(Z(Jndex)));
0358         lndex = Jndex(find(test2));
0359         Lndex = setdiff(Jndex,lndex);
0360     else
0361         Lndex = [];
0362     end
0363     index = INDEX(Lndex);
0364     mm = length(index);
0365 end
0366 gamma_variates = X.*b;
0367 
0368 
0369 function  gamma_variates = best_1978_algorithm(a,b)
0370 nn = length(a);
0371 mm = nn;
0372 bb = a-1;
0373 cc = 3*a-.75;
0374 UV = NaN(nn,2);
0375 Y  = NaN(nn,1);
0376 X  = NaN(nn,1);
0377 Z  = NaN(nn,1);
0378 W  = NaN(nn,1);
0379 index = 1:nn;
0380 INDEX = index;
0381 while mm
0382     UV(index,:) = rand(mm,2);
0383     W(index) = UV(index,1).*(1-UV(index,1));
0384     Y(index) = sqrt(cc(index)./W(index)).*(UV(index,1)-.5);
0385     X(index) = bb(index)+Y(index);
0386     jndex = index(find(X(index)>=0));
0387     Jndex = setdiff(index,jndex);
0388     if ~isempty(jndex)
0389         Z(jndex) = 64*W(jndex).*W(jndex).*W(jndex).*UV(jndex,2).*UV(jndex,2);
0390         kndex = jndex(find(Z(jndex)<=1-2*Y(jndex).*Y(jndex)./X(jndex)));
0391         Kndex = setdiff(jndex,kndex);
0392         if ~isempty(Kndex)
0393             lndex = Kndex(find(log(Z(Kndex))<=2*(bb(Kndex).*log(X(Kndex)./bb(Kndex))-Y(Kndex))));
0394             Lndex = setdiff(Kndex,lndex);
0395         else
0396             Lndex = [];
0397         end
0398         new_index = INDEX(Lndex);
0399         %mm = length(index);
0400     end
0401     index = union(new_index,INDEX(Jndex));
0402     mm = length(index);
0403 end
0404 gamma_variates = X.*b;