function r=gamrnd_marsg(a,b,n)
% analogous to gamrnd but uses the algorithm of Marsaglia and Tsang(2000)
% only use for a>1. 
%   R = GAMRND(A,B) returns a matrix of random numbers chosen   
%  
%   References:
%      [1]  Marsaglia, G. and Tsang, W.W.  (2000) "a simple method for generating gamma variates", 
%      ACM Transactions on Mth. Soft. 26, 3, 363-372

% Initialize r to zero.
r = zeros(1,n);

% Use marsaglia's rejection method for a > 1.
d=a-1/3;
x=randn(1,n);
u=rand(1,n);
v=(1+x/sqrt(9*d)).^3;
w=(log(u)<(x.^2/2+d-d*v+d*log(v)));
r=d*v(w);
r = b*r;