在R中手动模拟泊松过程

Question

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

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

• 让τ0 = 0。 • 生成独立同分布的指数随机变量ρ1，ρ2等等。 • 对于n = 1, 2等，令τn = ρ1 + … + ρn。 • 对于每个k = 0, 1等，当τk ≤ t < τk+1时，令Nt = k。

使用该方法，在区间[0,20]上生成一个λ为0.5的泊松过程（Nt）t。

生成10000个λ = 0.5的泊松过程（Nt）t的实现，并使用您的结果估计E（Nt）和Var（Nt）。将估计值与理论值进行比较。

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

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

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

    return (vec)
}

然后，我通过逐步累加 ρ 来创建 τ。

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值的。比如说，当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，则会得到以下输出：

.. 这是一个阶梯函数，但所有步骤都相等。

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

- user366312

---

那不是 rexp 的正确参数。你正在编写 C-as-R 代码。R 具有矢量化思维方式，并希望您一次生成所有随机抽样。请参阅 ?rexp 获取一些示例。 - MichaelChirico

或者至少你需要 rexp(n=1, rate=0.5)；rexp(0.5) 返回一个长度为0的数值向量... - Ben Bolker

Answer 1

看起来你在rhos函数中使用max1来指示指数分布采样的次数。我建议像这样：

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

这将在间隔中的每个整数时间生成泊松过程的值。我怀疑你的问题部分在于尝试将tau值放在y轴上，而不是已发生事件的计数。以下代码对我有效，可以产生类似于你的示例的随机外观阶梯形状。

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

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

- mcz

---

Ntk() 函数没问题吧？ - user366312

@user366312 Ntk函数正在计算val，它看起来是时间k之前最后一个事件发生的时间值。它还在一次k时计算此值。但是您需要计算在时间k之前发生的事件数量作为计数，并且需要一次性计算所有这些事件。目前编写Ntk的方式是每次查看新的时间值k时都会创建一个新的泊松过程。ppGen函数返回所有时间值k的计数值，然后您问题的第二部分可以使用ppGen的最后一个元素，即maxTime时的值。 - mcz

Answer 2

这里是另一种可能的实现方式。思路是生成一个等待时间向量（tau），并将其与我们正在等待的事件列表（max1）绘制在一起。

poi.process <- function(lambda,n){
  
  # initialize vector of total wait time for the arrival of each event:
  s<-numeric(n+1)
  # set S_0 = 0
  s[1] <-0
  # generate vector of iid Exp random variables:
  x <-replicate(n,rexp(1,lambda))
  # assign wait time to vector s in for loop:
  for (k in 1:n){
    
    s[k+1] <-sum(x[1:k])
    
  }
  # return vector of wait time
  return(s)
  
}

使用stepfun绘制它，我们会得到如下图所示的结果：

n<-20

lambda <-3

# simulate list of wait time:
s_list <-poi.process(lambda,n)

# plot function:
plot(stepfun(0:(n-1), s_list), 
     do.points = TRUE,
     pch = 16,
     col.points = "red",
     verticals = FALSE,
     main = 'Realization of a Poisson process with lambda = 3',
     xlab = 'Time of arrival',
     ylab = 'Number of arrivals')

泊松过程示例：

- sami c