r - 在 R 中不使用内置函数生成几何变量

Question

我正在学习 R 并且刚刚开始从不同的概率分布进行模拟。我试图从 生成一个随机变量几何分布 .我不愿意使用内置功能rgeom在 R 中做同样的事情。在其他方法中，我决定使用以下结果:

The Geometric distribution is the probability distribution of the number of Bernoulli trials needed to get one success.

我将生成伯努利变量的函数编写为 in my previous post命名为 ber , 与 p是成功概率。
现在要生成一个几何变量，我可以运行一个循环(我在这里使用了 repeat 循环)并不断生成一个伯努利变量，检查它是否成功等等......
我是这样做的:

geom<-function(p)
{
    x<-1
    repeat
    {
        if(ber(p)==1) break
        x<-x+1
    }
    return(x)
}

由 R 生成的一个输出是:

geom(0.035)
[1] 9

我想做什么以及我被困在哪里:
我想使用矢量化方法而不是这个 repeat循环方法来生成几何变量，但我似乎对如何实现这一点没有任何想法。在这里非常欢迎任何帮助/建议和提示。
注 1:至于为什么我定义了一个新函数 geom(p) ，我最初计划生成不是一个，而是许多(比如 n )几何变量。我想像这样定义一个函数，然后使用 replicate函数生成更多具有相同成功概率的相似变量更容易。
注 2:如果我以前的一些帖子看起来很相似(在某种意义上说，要求用一些向量方法替换循环的建议)，我真的很抱歉，但我是 R 的新手。我仍在学习矢量化方法。

Answer 1

在不使用循环的情况下，您可以利用几何随机变量具有可逆累积分布函数的事实并使用 inverse transform sampling .您的几何变量版本的 CDF(首次成功次数，而不是您计算首次成功前失败次数的版本)是:

F(x) = 1 - (1-p)^x,  x = 1, 2, 3, ....

我们可以这样求解 x:

x = log(1-F(x))/(log(1-p))

要使用反向采样，我们需要取其上限以使其离散:

geom <- function(n,p) {
  ceiling(log(runif(n))/log(1-p))
}

上面的代码利用了 if U 的事实。在 [0,1] 中是统一的那么 1-U .
例如:

> set.seed(42)
> geom(10,1/6)
 [1]  1  1  7  2  3  4  2 11  3  2
> mean(geom(100000,1/6))
[1] 6.00359

在编辑 一种更符合您目前所拥有的方法的方法是创建一个完整的随机变量向量，全部初始化为 0 ，并在实验仍在等待成功的意义上保持变量索引列表仍然“活跃”。在每一步，给这些事件变量加 1，概率为 1-p在下一步中保持一些索引处于事件状态:

geom2 <- function(n,p) {
  nums <- numeric(n)
  active <- 1:n
  k <- n
  while(k > 0){
    nums[active] <- nums[active]+1
    active <- active[runif(k) <= 1-p]
    k = length(active)
  }
  nums
}

这仍然相当快，尽管曾经 n与第一个定义相比，它明显落后于数百万。
我没有看到任何方法可以避免使用循环或逆变换。