73 lines
2.0 KiB
Matlab
73 lines
2.0 KiB
Matlab
function gb = gamm_rnd(nrow,m,k)
|
|
% PURPOSE: a matrix of random draws from the gamma distribution
|
|
%---------------------------------------------------
|
|
% USAGE: r = gamm_rnd(n,m,k)
|
|
% where: n = the size of the vector drawn
|
|
% m = a parameter such that the mean of the gamma = m/k
|
|
% k = a parameter such that the variance of the gamma = m/(k^2)
|
|
% note: m=r/2, k=2 equals chisq r random deviate
|
|
%---------------------------------------------------
|
|
% RETURNS:
|
|
% r = an n x 1 vector of random numbers from the gamma distribution
|
|
% --------------------------------------------------
|
|
% SEE ALSO: gamm_inv, gamm_pdf, gamm_cdf
|
|
%---------------------------------------------------
|
|
% NOTE: written by: Michael Gordy, 15 Sept 1993
|
|
% mbgordy@athena.mit.edu
|
|
%---------------------------------------------------
|
|
% REFERENCES: Luc Devroye, Non-Uniform Random Variate Generation,
|
|
% New York: Springer Verlag, 1986, ch 9.3-6.
|
|
|
|
if nargin ~= 3
|
|
error('Wrong # of arguments to gamm_rnd');
|
|
end;
|
|
|
|
ncol = 1;
|
|
gb=zeros(nrow,ncol);
|
|
if m<=1
|
|
% Use RGS algorithm by Best, p. 426
|
|
c=1/m;
|
|
t=0.07+0.75*sqrt(1-m);
|
|
b=1+exp(-t)*m/t;
|
|
for i1=1:nrow
|
|
for i2=1:ncol
|
|
accept=0;
|
|
while accept==0
|
|
u=rand; w=rand; v=b*u;
|
|
if v<=1
|
|
x=t*(v^c);
|
|
accept=((w<=((2-x)/(2+x))) | (w<=exp(-x)));
|
|
else
|
|
x=-log(c*t*(b-v));
|
|
y=x/t;
|
|
accept=(((w*(m+y-m*y))<=1) | (w<=(y^(m-1))));
|
|
end
|
|
end
|
|
gb(i1,i2)=x;
|
|
end
|
|
end
|
|
else
|
|
% Use Best's rejection algorithm XG, p. 410
|
|
b=m-1;
|
|
c=3*m-0.75;
|
|
for i1=1:nrow
|
|
for i2=1:ncol
|
|
accept=0;
|
|
while accept==0
|
|
u=rand; v=rand;
|
|
w=u*(1-u); y=sqrt(c/w)*(u-0.5);
|
|
x=b+y;
|
|
if x >= 0
|
|
z=64*(w^3)*v*v;
|
|
accept=(z<=(1-2*y*y/x)) ...
|
|
| (log(z)<=(2*(b*log(x/b)-y)));
|
|
end
|
|
end
|
|
gb(i1,i2)=x;
|
|
end
|
|
end
|
|
end
|
|
gb=gb/k;
|
|
|
|
|