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

Question

I have the following algorithm

Generate U~(0,1)
X=-beta*ln(u_1*u_2*...*u_alpha)
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?

A few obvious things... `prod(u[i])` just gives `u[i]` which is a single value. You probably want `prod(u)` - and don't need the `i` loop. You are also only generating one set of random `u` variables - you might want to do that inside the loop to get a different set each time. — Andrew Gustar, Apr 12 '19 at 18:05

score 1 · Accepted Answer · answered Apr 12 '19 at 18:38

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).

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

1 Answers1