手动模拟R中的泊松过程

Question

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

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

•设τ0= 0。 •生成i.i.d.指数随机变量ρ1，ρ2,. 。 .. •设τn=ρ1+。 。 。 +ρn表示n = 1,2，.... 。 。 。 •对于每个k = 0,1，.... 。对于τk≤t<τk+ 1，设Nt = k。

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

我试图解决的问题

首先，我在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的值时找到Nt = 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，我会得到以下输出：

..这是一个阶梯功能，但所有步骤都相同。

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

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)
}

现在您已经知道了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")