Почему код в r о генерации гамма-случайных переменных не возвращает ожидаемый результат? - PullRequest
Новая
       76

Почему код в r о генерации гамма-случайных переменных не возвращает ожидаемый результат?

1 голос
/ 12 апреля 2019

У меня есть следующий алгоритм

  1. Генерация U ~ (0,1)
  2. X = -бета * ln (u_1 * u_2 * ... * u_alpha)
  3. повторите шаг 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)

Выход в 1000 раз больше того же числа .2139196когда я вызываю функцию с альфа = 3 и бета = 4

По крайней мере, нет ошибки, но она не должна возвращать тот же номер,

что я делаю неправильнотекущий код?

1 Ответ

1 голос
/ 12 апреля 2019

После комментария @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 и масштаб параметр beta, в то время как beta в rgamma(N, alpha, beta) является параметром 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)

(коэффициент обратный к шкале).

...