CDF 很容易反转。 反演公式导致必须是最简单和最方便的可能解决方案之一。

首先观察结果的成正比。因此，如果我们在和 _{\ max }=\ sum_ =，我们只需要找到最大的$k$$0 \le k \le n$$e^{-b k}$$q$$0$$q_{\max}=\sum_{k=0}^{n} e^{-b k}$$(1 - e^{-b(n+1)})/(1 - e^{-b})$$k$

$$q \ge \sum_{i=0}^{k} e^{-bi} = \frac{1 - e^{-(k+1)b}}{1 -e^{-b}}.$$

简单代数给出了解决方案

$$k = - \text{ceiling}\left(\frac{\log(1 - q (1-e^{-b}))}{b}\right).$$

这是一个与所有其他随机数生成器一样构造的R实现：它的第一个参数指定要生成多少个iid值，其余参数命名参数（ as和 as ）：$b$b$n$n.max

rgeom.truncated <- function(n=1, b, n.max) {
  a <- 1 - exp(-b)
  q.max <- (1 - exp(-b*(n.max+1))) / a
  q <- runif(n, 0, q.max)
  return(-ceiling(log(1 - q*a) / b))
}

作为其使用示例，让我们根据此分布生成一百万个随机变量：

b <- 0.001
n.max <- 3500
n.sim <- 10^6
set.seed(17)
system.time(sim <- rgeom.truncated(n.sim, b,n.max))

（秒。）$0.10$

h <- hist(sim+1, probability=TRUE, breaks=50, xlab="Outcome+1")
pmf <- exp(-b * (0: n.max)); pmf <- pmf / sum(pmf)
lines(0:n.max, pmf, col="Red", lwd=2)

（添加到每个值以创建更好的直方图：的过程有一个特质（=错误），其中当左端点设置为零时，第一个条形太高。）红色曲线是参考分布这个模拟试图重现。让我们用卡方检验评估拟合优度：$1$Rhist

observed <- table(sim)
expected <- n.sim * pmf
chi.square <- (observed-expected)^2 / expected
pchisq(sum(chi.square), n.max, lower.tail=FALSE)

p 值为：非常合适。$0.84$

您正在处理的截断几何分布。有多种方法可以解决这个问题。$p = 1-e^{-b}$

我会在不同的情况下建议不同的选择；一些选项将涉及从几何模拟并在其超出范围时重新生成，采用适当的截断指数的整数部分（如此处），或使用针对有限范围内的离散分布量身定制的几种快速技术中的任何一种。鉴于很大，取截断指数的下限可能相对较快，但它是否是最佳选择还取决于。$n$$b$

这是关于 math.SE的相关问题

在我尝试具体建议之前，的典型值范围是多少？$b$

首先，请注意如果是连续的，它将与指数分布有关。然后，您可以做的是从截断的指数分布进行模拟并获取（整数部分）观测值。$P(x)\propto e^{-bx}$$x$floor()

截断指数的 cdf 是

$$F(x;n,b)= \dfrac{1-e^{-bx}}{1-e^{-bn}}.$$

然后，如果我们使，我们得到。如果很大，则建议近似。$F(x;n,b)=u$$x=-\dfrac{1}{b}\log[1-u(1-e^{-bn})]$$bn$$e^{-bn}\approx 0$$x\approx -\dfrac{1}{b}\log[1-u]$

rweirdp <- function(ns,n,b){
u <- runif(ns)
samp <- - log(1-u*(1-exp(-n*b)))/b
return(floor(samp))
}

rweirdp(1000,10,1)

从目标分布中采样的一种方法是$p(k)\propto \exp\{-bk\}$

1. 运行 Metropolis-Hastings 实验以确定分布的（有趣的）支持，即它集中在的哪个子集中；$\{0,1,\ldots,n\}$

   metro=function(N,b,n){
   x=sample(0:n,N,rep=TRUE)
   for (t in 2:N){
     x[t]=prop=x[t-1]+sample(c(-1,1),1)
   
     if ((prop<0)||(prop>n)||(log(runif(1))>b*(x[t]-prop)))
         x[t]=x[t-1]
     }
   return(x)
   }
   

2. 使用如此确定的支持，来计算精确概率为以避免溢出。$\{k_0,\ldots,k_1\}$$p(k)\propto \exp\{-bk+bk_0\}$

更新：考虑更多时，由于在 k 中减小，分布的有效支持将始终从开始。如果很大，这种支持将很快结束，在这种情况下，无关紧要，因为永远不会访问如果非常小，则 pdf 几乎是平坦的，这意味着可以使用上的均匀分布作为接受-拒绝提议。并在接受步骤中使用日志以避免溢出。$p(\cdot)$$k_0=0$$b$$n$$k$$b$$\{0,1,\ldots,n\}$