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

Question

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

生成 U
定义 x=0,P=(1-p)^n 和 F=P
如果 U < F，则 X=x 并停止。
定义 P=(n-x)pP/(x+1)(1-p), F=F+P 和 x=x+1
回到第 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)
}

但是，在编译代码时，我收到十条警告消息，它们都说：

Warning messages: 1: In if (U < FF) { : the condition has length > 1 and only the first element will be used

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

Answer 1

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

在4中应该是P=(n-x) * p * P / ((x+1)*(1-p)) 所以 (1-p) 是分母而不是分子。
U 是一个随机数，但是您使用 runif(n) 创建了 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

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

An algorithm in r about binomial random variables

random

simulation

r