r - 在 R 中手动模拟泊松过程

Question

以下问题告诉我们从ρ（到达间隔时间）和τ（到达时间）逐步生成泊松过程。

讲座中提出的理论成果之一给出了以下模拟泊松过程的直接方法：

• 令τ ₀ = 0。
• 生成iid 指数随机变量ρ1, ρ2, 。. ..
• 令τ _n = ρ ₁ + 。. . + ρ _n对于 n = 1, 2, . . . .
• 对于每个 k = 0, 1, 。. ., 令 N _t = k 对于 τ _k ≤ t < τ _k+1。

1. 使用此方法，在区间 [0, 20] 上生成一个泊松过程 (N _t ) _{t的实现，其中 λ = 0.5。}
2. _{生成 λ = 0.5 的泊松过程 (N t} ) _t的 10000 个实现，并使用您的结果来估计 E(N _t ) 和 Var(N _t )。将估计值与理论值进行比较。

我尝试的解决方案：

首先，我使用R 中的函数生成了ρ的值。rexp()

rhos <-function(lambda, max1)
{
    vec <- vector()

    for (i in 1:max1) 
    {
        vec[i] <- rexp(0.5)
    }

    return (vec)
}

然后，我通过对ρ s进行累进求和来创建τ s。

taos <- function(lambda, max)
{
    rho_vec <- rhos(lambda, max)
    #print(rho_vec)

    vec <- vector()
    vec[1] <- 0
    sum <- 0
    for(i in 2:max)
    {
        sum <- sum + rho_vec[i]
        vec[i] <- sum
    }

    return (vec)
}

以下函数用于在给定 k 值时找到N _t =k的值。比如说，它是7等等。

Ntk <- function(lambda, max, k)
{
    tao_vec <- taos(lambda, max)
    val <- max(tao_vec[tao_vec < k])
}

y <- taos(0.5, 20)
x <- seq(0, 20-1, by=1)

plot(x,y, type="s")

输出：

如您所见，泊松过程的图是空白的，而不是楼梯。

如果我更改rexp为exp，我会得到以下输出：

..这是一个楼梯功能，但所有步骤都是平等的。

为什么我的源代码没有产生预期的输出？

标签： rmontecarlopoissonstochastic-process

Answer 1

看起来您正在使用 max1 来指示对rhos函数中的指数分布进行采样的次数。我会推荐这样的东西：

rhosGen <- function(lambda, maxTime){
  rhos <- NULL
  i <- 1
  while(sum(rhos) < maxTime){
    samp <- rexp(n = 1, rate = lambda)
    rhos[i] <- samp
    i <- i+1
  }
  return(head(rhos, -1))
}

这将继续从指数采样，直到这些保持时间的总和大于给定间隔的长度。head删除最后一个样本，以便我们跟踪的所有事件都肯定发生在我们感兴趣的时间间隔内。从这里你必须通过总结之前的持有时间 (rhos) 来生成 taos：

taosGen <- function(lambda, maxTime){
  rhos <- rhosGen(lambda, maxTime)
  taos <- NULL
  cumSum <- 0
  for(i in 1:length(rhos)){
    taos[i] <- sum(rhos[1:i])
  }
  return(taos)
}

现在您已经了解了我们知道时间间隔 (0,maxTime) 中每个事件发生的时间。这导致我们通过查找时间间隔中每个 t 的 Nt 值来生成相关的泊松过程：

ppGen <- function(lambda, maxTime){
  taos <- taosGen(lambda, maxTime)
  pp <- NULL
  for(i in 1:maxTime){
    pp[i] <- sum(taos <= i)
  }
  return(pp)
}

这会在间隔中的每个整数时间生成泊松过程的值。我怀疑您的部分问题是试图将 tao 值放在 y 轴上，而不是已经发生的事件计数。以下代码对我有用，可以生成一个随机的楼梯箱，类似于您的示例。

y <- ppGen(0.5, 20)
x <- seq(0, 20-1, by=1)

plot(x,y, type="s")