From 01ba142329b69d274d08479cb8b9160cb7a5bd93 Mon Sep 17 00:00:00 2001
From: stepan <stepan@ac1d8469-bf42-47a9-8791-bf33cf982152>
Date: Tue, 17 Feb 2009 16:00:40 +0000
Subject: [PATCH] v4.1:: Changes related to the gamma variate generator. +
 Changed headers. + Allows vector of coefficients when calling the Beta random
 generator.

git-svn-id: https://www.dynare.org/svn/dynare/trunk@2422 ac1d8469-bf42-47a9-8791-bf33cf982152
---
 matlab/distributions/toolbox/betarnd.m | 41 ++++++++++++++++++--------
 matlab/distributions/toolbox/gamrnd.m  | 24 ++++++++++-----
 2 files changed, 45 insertions(+), 20 deletions(-)

diff --git a/matlab/distributions/toolbox/betarnd.m b/matlab/distributions/toolbox/betarnd.m
index 281902898..e3a906524 100644
--- a/matlab/distributions/toolbox/betarnd.m
+++ b/matlab/distributions/toolbox/betarnd.m
@@ -1,11 +1,19 @@
 function rnd = betarnd(a, b)
-% BETARND  Random samples from the Beta distribution
-%  RND = betarnd(A,B) returns a random sample from the
-%  Beta distribution with parameters A and B (i.e. mean of
-%  the distribution is A/(A+B) and variance is
-%  A*B/(A+B)^2/(A+B+1) ).
+% This function produces independent random variates from the Beta distribution.
+%
+%  INPUTS 
+%    a       [double]    n*1 vector of positive parameters.
+%    b       [double]    n*1 vector of positive parameters.
+%
+%  OUTPUT 
+%    rnd     [double]    n*1 vector of independent variates from the beta(a,b) distribution.
+%                        rnd(i) is beta distributed with variance a(i)/(a(i)+b(i)) and 
+%                        variance a(i)b(i)/(a(i)+b(i))^2/(a(i)+b(i)+1).    
+%  
+%  ALGORITHMS     
+%    Described and Devroye (1986, chapter 9).
 
-% Copyright (C) 2008 Dynare Team
+% Copyright (C) 2008-2009 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -26,13 +34,20 @@ if (nargin ~= 2)
   error('betarnd: you must give two arguments');
 end
 
-if (~isscalar(a) || ~isscalar(b))
-  error('betarnd: A and B should be scalar parameters');
+if (any(a<0)) || (any(b<0)) || (any(a==Inf)) || (any(b==Inf))
+    error('betarnd:: Input arguments must be finite and positive!');
 end
 
-if (a <= 0 || a == Inf || b <= 0 || b == Inf)
-  rnd = NaN;
-else
-  x = gamrnd(a, 1);
-  rnd = x/(x+gamrnd(b, 1));
+[ma,na] = size(a);
+[mb,nb] = size(b);
+
+if ma~=mb || na~=nb
+    error('betarnd:: Input arguments must have the same size!');
 end
+
+if na~=1
+    error('betarnd:: Input arguments must be column vectors');
+end
+
+x = gamrnd(a,ones(ma,1));
+rnd = x./(x+gamrnd(b, ones(mb,1)));
\ No newline at end of file
diff --git a/matlab/distributions/toolbox/gamrnd.m b/matlab/distributions/toolbox/gamrnd.m
index 79cb37607..5260832bb 100644
--- a/matlab/distributions/toolbox/gamrnd.m
+++ b/matlab/distributions/toolbox/gamrnd.m
@@ -1,12 +1,19 @@
 function rnd = gamrnd(a,b,method)
-% GAMRND  Random samples from the Gamma distribution
-%  RND = gamrnd(A,B) returns a random sample from the
-%  Gamma distribution with parameters A and B (i.e. mean of
-%  the distribution is A*B and variance is A*B^2).
+% This function produces independent random variates from the Gamma distribution.
 %
-% Algorithm of Bauwens, Lubrano & Richard (page 316)
+%  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 variance 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).
 
-% Copyright (C) 2006-2008 Dynare Team
+% Copyright (C) 2006-2009 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -36,7 +43,10 @@ if nargin==2
     end
     if ~strcmpi(method,'BauwensLubranoRichard')
         Devroye.big = 'Best'; % 'Cheng' , 'Best' 
-        % REMARK: The first algorithm (Cheng) is still producing obviously wrong simulations. 
+        % REMARK 1: The first algorithm (Cheng) is still producing obviously wrong simulations.
+        % REMARK 2: The second algorithm seems slightly slower than the algorithm advocated by Bauwens, 
+        %           Lubrano and Richard, but the comparison depends on the value of a (this should be 
+        %           investigated further).
     end
 end