Why code in r about generation of Gamma random variables does not return the expected output?

Question

I have the following algorithm

1. Generate U~(0,1)
2. X=-beta*ln(u_1*u_2*...*u_alpha)
3. repeat step 1 an step 2 N times

I think in r is as follows

Gam<-function(alpha,beta){
N<-1000
u<-runif(alpha)
x<-rep(0,alpha)
X<-rep(0,N)

for(j in 1:N){
for(i in 1:alpha){
x[i]<--beta*log(prod(u[i]))
X[j]<-x[i]
}
}
return(X)}

#Comparation with the predefined gamma function

y<-rgamma(alpha,beta)
par(mfrow=c(1,2))

hist(X,freq=F,nclass=10)
hist(y,freq=F,nclass=10)

The output is 1000 times the same number .2139196 when I call the function with alpha=3 and beta=4

At least there is no error but it should have not return the same number,

what am I doing wrong in the current code?

Answer 1

Following @Andrew Gustar's comment, here is the code you want:

Gam <- function(alpha, beta){
  N <- 1000
  X <- rep(0,N)
  for(j in 1:N){
    u <- runif(alpha)
    X[j] <- -beta * log(prod(u))
  }
  return(X)
}

Moreover, this code samples the Gamma distribution with shape parameter alpha and scale parameter beta, while beta in rgamma(N, alpha, beta) is the rate parameter. So, to compare with rgamma you have to do:

alpha <- 3; beta <- 4
y1 <- Gam(alpha, beta)
y2 <- rgamma(1000, alpha, scale=beta)

par(mfrow=c(1,2))
hist(y1, freq=F, nclass=10)
hist(y2, freq=F, nclass=10)

or

y2 <- rgamma(1000, alpha, 1/beta)

(the rate is the inverse of the scale).