Моделирование случайных величин с использованием метода принятия-отклонения

Question

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

Шаг 1. Имитировать значение Y с помощью qj = P (Y = j)

Шаг 2. Сгенерируйте равномерную переменную

Шаг 3. Если U <= Pj / (c * qj), то X = j и остановка. В противном случае вернитесь к шагу 1. </p>

И конкретный пример:

X = 1,2,3,4 с P1 = .2, P2 = .15, P3 = .25, P4 = .4

Создание значений для X

1. Пусть Y ~ UD (1,4)

2. с = .4 / .25

Вот мой подход к реализации этого алгоритма в R:

f<-function(){
n<-4

probY<-c(.25,.25,.25,.25)
probX<-c(2,.15,.25,.4)

X<-rep(0,4)

U<-runif(n)

c<-.4/.25

for(j in 1:4)
{
if(U[j]<= probX[j]/(c*probY[j]))
   X<-j

}


return(X)

}

Вывод 4 Не думаю, что это правильно. Я не уверен, должен ли я написать Y<-runif(n,1,4) или эту строку probY<-c(.25,.25,.25,.25). Также эта строка «В противном случае вернитесь к шагу 1.» отсутствует в цикле, хотя всегда одинаково .25

Может кто-нибудь помочь, пожалуйста?

Answer 1

Я думаю, что проблема здесь в некоторой путанице с тем, как работает алгоритм.

Чтобы сгенерировать одиночное значение из вашего дистрибутива (X = 1,2,3,4, где P (X = 1) = .2, P (X = 2) = .15, P (X = 3) = .25, P (X = 4) = .4), нам нужно следовать шагам алгоритма. Предполагая, что мы выбираем c = .4 / .25 :

1.Generate y из Y ~ UD (1,4).

2. Генерировать u из U ~ U (0,1).

3. Проверьте, нет ли u≤f (y) / cg (y) . Если равно , определите x = y и все готово! Если не , вернитесь к шагу 1.

В коде, который вы дали, вы фактически никогда не генерируете переменную y . Вот функция, которая должна работать вместо вас! Надеюсь, мои комментарии достаточно хорошо это объясняют!

accRej <- function(){

  #The probabilities for generating a r.v. from X
  probX <- c(.2,.15,.25,.4)

  #The Value for c
  c <- .4/.25

  #x is a placeholder for our final value of x 
  #and i is a counter variable which will allow us to stop the loop when we complete step 3
  x <- numeric(1)
  i <- 1

  #Now, start the loop!
  while(i <= 1){
    #Step 1, get y
    y <- sample(1:4,1)
    #Step 2, get u
    u <- runif(1)
    #Step 3, check to see if the inequality holds
    #If it does, assign y to x and add 1 to i (making it 2) to stop the while loop
    #If not, do nothing and try again!
    if(u <= probX[y]/c*.25){
      x[i] <- y
      i <- i+1
    }
  }
  #Return our value of x
  return(x)
}

Обратите внимание, что в этом коде probX[i] равно f (y) в нашем алгоритме, и, поскольку Y ~ UD (1,4), .25 = g (y) всегда. Надеюсь, это поможет!

Кроме того, здесь приведен код для генерации n случайных величин с помощью этого метода. По сути, это то же самое, что и выше, только с возможностью изменить 1 на n.

accRej <- function(n){

  #The probabilities for generating a r.v. from X
  probX <- c(.2,.15,.25,.4)

  #The Value for c
  c <- .4/.25

  #x is a placeholder for our final value of x 
  #and i is a counter variable which will allow us to stop the loop when we complete step 3
  x <- numeric(n)
  i <- 1

  #Now, start the loop!
  while(i <= n){
    #Step 1, get y
    y <- sample(1:4,1)
    #Step 2, get u
    u <- runif(1)
    #Step 3, check to see if the inequality holds
    #If it does, assign y to x and add 1 to i
    #If not, do nothing and try again!
    if(u <= probX[y]/c*.25){
      x[i] <- y
      i <- i+1
    }
  }
  #Return our value of x
  return(x)
}