I only read section 4, and translated it as follows:
<br />
rmvgamma=function(n,p,shape,rate=1/scale,rho=0,scale=1)<br />
{<br />
stopifnot(length(rho)==1,length(shape)==1,length(rate)==1,<br />
length(scale)==1,length(n)==1,length(p)==1,<br />
rate>0, shape>0, rho>=0, rho<=1, scale>0, n>=0, p>0<br />
)<br />
n=round(n);p=ceiling(p)<br />
theta=rate*rho/(1-rho)<br />
k=rnbinom(n,shape,rate/(rate+theta))<br />
matrix(rgamma(p*n,shape+k,rate+theta),n)<br />
}<br />
It _SEEMS_ to be working:
<br />
> x=rmvgamma(100000,2,1.9,4,.5)<br />
> hist(x,pr=T,br=100)<br />
> curve(dgamma(x,1.9,4),0,5,add=T)<br />
> ks.test(x,pgamma,shape=1.9,rate=4)<br />
<br />
One-sample Kolmogorov-Smirnov test<br />
<br />
data: x <br />
D = 0.0024, p-value = 0.2159<br />
alternative hypothesis: two-sided <br />
> cor(x)<br />
[,1] [,2]<br />
[1,] 1.0000000 0.5013881<br />
[2,] 0.5013881 1.0000000<br />
> qqplot(x[,1],x[,2]);abline(0,1)<br />
I was surprised at the size parameter of neg-binom not necessary to be an integer...