为何生成泊松分布随机数的此R代码仅适用于小lambda值？

你这段代码用的是逆变换采样法的实现，靠累积概率和与均匀随机数对比来生成泊松变量，但lambda值大时会失效，核心原因有两个：

• 数值下溢导致循环卡死
  当lambda很大时，初始计算的p = exp(-lambda)会极小，直接跌到R浮点数精度的下限变成0。后续p = p * lambda / i的计算结果始终是0，suma就停留在初始的0（或接近0的极小值），永远无法超过0-1之间的U[j]，while循环会无限执行，程序直接卡死。

• 累积概率的数值误差失控
  就算初始p没完全下溢，大lambda对应的泊松分布，只有当i接近lambda时概率才会显著，前面小i对应的p都极小。累加这些极小值时，浮点运算的累加误差会被不断放大，导致suma无法准确反映真实的累积概率。最终要么循环停不下来，要么生成的i值完全不符合泊松分布。

如果要适配大lambda场景，要么用正态近似法（lambda足够大时，泊松分布近似N(lambda, lambda)，再做整数修正），要么直接用R内置的rpois()函数——它内部会根据lambda大小自动切换优化算法，比如小lambda用逆变换，大lambda用更稳定的方法。

原代码：

poisson_numbers <- function(lambda, n) {
  x <- numeric(n)
  U <- runif(n, 0, 1)
  print(U)
  for (j in 1:n) {
    i <- 0
    p <- exp(-lambda)
    suma <- p
    while (suma < U[j]) {
      i <- i + 1
      p <- p * lambda / i
      suma <- suma + p
    }
    x[j] <- i
  }
  return(x)
}

内容的提问来源于stack exchange，提问作者cosmic