Алгоритм в r о биномиальных случайных величинах

Question

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

1. Generate U
2. Определить x = 0, P = (1-p) ^n и F = P
3. Если U
4. Определить P = (nx) p P / (x + 1) (1-p), F = F + P и x = x + 1
5. Вернитесь к шагу 3.

Насколько мне известно, r будет

varBinom<-function(n,p)
{
  U<-runif(n)
  x<-0
  P<-(1-p)^n
  FF<-P

  for(i in 1:n)
  {
    if(U<FF)
    {
      X<-x
      break
    }
    P<-(n-x)*p*P/(x+1)*(1-p)
    FF<-FF+P
    x<-x+1
  }
  return(x)
}

Однако при компиляции кода я получаю десять предупреждающих сообщений, все они говорят:

Предупреждающие сообщения: 1: In if (U 1 и будет использоваться только первый элемент

Почему это происходит?Как я могу исправить код?

---

Answer 1

Я думаю, что вы сделали две маленькие ошибки.

1. В 4 это должно быть P = (nx) * p * P / ( (x + 1) * (1-p) ) поэтому (1-p) находится в знаменателе, а не в числителе.
2. U - одно случайное число, но вы создаете n различных случайных чисел, используя runif(n).Вот почему вы получаете предупреждение.

Итак, вот исправленный алгоритм:

varBinom<-function(n, p)
{
  U <- runif(1)
  x <- 0
  P <- (1-p)^n
  FF <- P

  for(i in 1:n)
  {
    if(U<FF) return(x)
    P <- (n-x) * p * P/((x+1)*(1-p))
    FF <- FF+P
    x <- x+1
  }
  return(x)
}

Вот результат, когда вы вызываете функцию 15 раз:

set.seed(1)
replicate(15, varBinom(10, 1/2))
[1] 4 4 5 7 4 7 7 6 6 3 4 4 6 5 6