为什么 r 中有关 Gamma 随机变量生成的代码没有 return 预期输出？

Question

我有以下算法

生成U~(0,1)
X=-beta*ln(u_1*u_2*...*u_alpha)
重复步骤 1 和步骤 2 N 次

我觉得在r中是这样的

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

至少没有错误，但应该不是return相同的数字，

我在当前代码中做错了什么？

Answer 1

根据@Andrew Gustar 的评论，这是您想要的代码：

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

此外，此代码使用形状参数 alpha 和 scale 参数 beta 对 Gamma 分布进行采样，而 beta 在 [=16] =] 是 rate 参数。因此，要与 rgamma 进行比较，您必须这样做：

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)

或

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

（比率是比例的倒数）

为什么 r 中有关 Gamma 随机变量生成的代码没有 return 预期输出？

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

algorithm

simulation

statistics

r

gamma-distribution