关于二项式随机变量的r算法

Question

我有以下算法（基于this algorithm）

1. 生成U.
2. 定义x = 0，P =（1-p）^ n和F = P.
3. 如果U <F，则X = x并停止。
4. 定义P =（n-x）pP /（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：如果（U <FF）{：条件长度> 1且仅使用第一个元素

为什么会这样？我该如何修复代码？

---

Answer 1

我认为你犯了两个小错误。

1. 在4中它应该是P =（n-x）* 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